context("sum of llcont") test_that("lavaan object", { if (isTRUE(require("lavaan"))) { HS.model <- 'visual =~ x1 + x2 + x3 textual =~ x4 + x5 + x6 speed =~ x7 + x8 + x9 ' fit1 <- cfa(HS.model, data=HolzingerSwineford1939) fit2 <- cfa(HS.model, data=HolzingerSwineford1939, group="school") fit3 <- cfa(HS.model, data = HolzingerSwineford1939, group = "school", group.equal = c("loadings")) fit4 <- cfa(HS.model, data = HolzingerSwineford1939, group = "school", group.equal = c("loadings"), group.partial = c("visual=~x2", "x7~1")) expect_equal(sum(llcont(fit1)), as.numeric(logLik(fit1))) expect_equal(sum(llcont(fit2)), as.numeric(logLik(fit2))) expect_equal(round(sum(llcont(fit3)) - as.numeric(logLik(fit3)), 8), 0L) expect_equal(sum(llcont(fit4)), as.numeric(logLik(fit4))) HS.model2 <- 'visual =~ x1 + 0.5*x2 + c(0.6, 0.8)*x3 textual =~ x4 + start(c(1.2, 0.6))*x5 + a*x6 speed =~ x7 + x8 + x9' fit5 <- cfa(HS.model2, data=HolzingerSwineford1939, group="school") expect_equal(round(sum(llcont(fit5)) - as.numeric(logLik(fit5)), 8), 0L) HS.model3 <- 'visual =~ x1 + x2 + c(v3,v3)*x3 textual =~ x4 + x5 + x6 speed =~ x7 + x8 + x9' fit6 <- cfa(HS.model3, data=HolzingerSwineford1939, group="school") expect_equal(sum(llcont(fit6)), as.numeric(logLik(fit6))) ## a more stringent test with missing, fixed.x, empty case: set.seed(1234) pop.model <- ' f =~ 0.7*y1 + 0.7*y2 + 0.7*y3 + 0.7*y4 + 0.7*y5 f ~ (-2.3)*x1 + 0.8*x2 y1 ~ 0.2*x2 y3 ~ 0.7*x1 ' Data <- simulateData(pop.model, sample.nobs=100) model <- ' f =~ y1 + y2 + y3 + y4 + y5 f ~ x1 + x2 y1 ~ x2 y3 ~ x1 ' obs <- rbinom(prod(dim(Data)), 1, .9) Data <- Data*obs Data[Data==0] <- NA Data[95,] <- NA fit7 <- sem(model, data=Data, fixed.x=TRUE, meanstructure=TRUE, missing='ml') expect_equal(sum(llcont(fit7)), as.numeric(logLik(fit7))) } }) test_that("glm object", { if (isTRUE(require("faraway")) && isTRUE(require("MASS"))) { ## binomial bin1 <- glm(formula=am ~ hp + wt, data=mtcars, family=binomial) bin2 <- glm(cbind(Menarche, Total-Menarche) ~ Age, family=binomial(logit), data=menarche) expect_equal(sum(llcont(bin1)), as.numeric(logLik(bin1))) expect_equal(sum(llcont(bin2)), as.numeric(logLik(bin2))) ## quasibinomial qbin1 <- glm(formula=am ~ hp + wt, data=mtcars, family=quasibinomial) qbin2 <- glm(cbind(Menarche, Total-Menarche) ~ Age, family=quasibinomial, data=menarche) expect_equal(sum(llcont(qbin1)), as.numeric(logLik(qbin1))) expect_equal(sum(llcont(qbin2)), as.numeric(logLik(qbin2))) ## gaussian gau1 <- glm(Species ~ Area + Elevation + Nearest + Scruz + Adjacent, data=gala, family=gaussian) gau2 <- glm(Species ~ Area + Elevation + Nearest, data=gala, family=gaussian) expect_equal(sum(llcont(gau1)), as.numeric(logLik(gau1))) expect_equal(sum(llcont(gau2)), as.numeric(logLik(gau2))) ## inverse.gaussian invGau1 <- glm(actual ~ projected-1, family=inverse.gaussian(link="identity"), cpd) expect_equal(sum(llcont(invGau1)), as.numeric(logLik(invGau1))) ## Gamma clotting <- data.frame(u = c(5,10,15,20,30,40,60,80,100), lot1 = c(118,58,42,35,27,25,21,19,18), lot2 = c(69,35,26,21,18,16,13,12,12)) gam1 <- glm(lot1 ~ log(u), data = clotting, family = Gamma) expect_equal(sum(llcont(gam1)), as.numeric(logLik(gam1))) ## poisson counts <- c(18,17,15,20,10,20,25,13,12) outcome <- gl(3,1,9) treatment <- gl(3,3) d.AD <- data.frame(treatment, outcome, counts) pois1 <- glm(counts ~ outcome + treatment, family = poisson) pois2 <- glm(counts ~ outcome, family = poisson) expect_equal(sum(llcont(pois1)), as.numeric(logLik(pois1))) expect_equal(sum(llcont(pois2)), as.numeric(logLik(pois2))) ## quasipoisson qpois1 <- glm(counts ~ outcome + treatment, family = quasipoisson) qpois2 <- glm(counts ~ outcome, family = quasipoisson) expect_equal(sum(llcont(qpois1)), as.numeric(logLik(qpois1))) expect_equal(sum(llcont(qpois2)), as.numeric(logLik(qpois2))) ## negative-binomial: MASS nb1 <- glm.nb(Days ~ Sex/(Age + Eth*Lrn), data = quine) expect_equal(sum(llcont(nb1)), as.numeric(logLik(nb1))) } }) test_that("clm object", { if (isTRUE(require("ordinal")) && isTRUE(require("MASS"))) { clm1 <- clm(rating ~ temp * contact, data = wine) clm2 <- update(clm1, ~.-temp:contact) clm3 <- update(clm1, link = "logit") clm4 <- update(clm1, link = "probit") clm5 <- update(clm1, link = "loglog") clm6 <- update(clm1, link = "cloglog") clm7 <- update(clm1, link = "cauchit") clm8 <- update(clm1, threshold = "symmetric") clm9 <- update(clm1, threshold = "equidistant") clm10 <- clm(Sat ~ Infl + Type + Cont, weights = Freq, data = housing) expect_equal(sum(llcont(clm1)), as.numeric(logLik(clm1))) expect_equal(sum(llcont(clm2)), as.numeric(logLik(clm2))) expect_equal(sum(llcont(clm3)), as.numeric(logLik(clm3))) expect_equal(sum(llcont(clm4)), as.numeric(logLik(clm4))) expect_equal(sum(llcont(clm5)), as.numeric(logLik(clm5))) expect_equal(sum(llcont(clm6)), as.numeric(logLik(clm6))) expect_equal(sum(llcont(clm7)), as.numeric(logLik(clm7))) expect_equal(sum(llcont(clm8)), as.numeric(logLik(clm8))) expect_equal(sum(llcont(clm9)), as.numeric(logLik(clm9))) expect_equal(sum(llcont(clm10)), as.numeric(logLik(clm10))) } }) test_that("hurdle object", { if (isTRUE(require("pscl"))) { hurdle1 <- hurdle(formula = art ~ ., data = bioChemists) hurdle2 <- hurdle(formula = art ~ ., data = bioChemists, separate=FALSE) hurdle3 <- hurdle(art ~ ., data = bioChemists, zero = "geometric") hurdle4 <- hurdle(art ~ fem + ment, data = bioChemists, dist = "negbin", zero = "negbin") hurdle5 <- hurdle(art ~ ., data = bioChemists, dist = "negbin") expect_equal(sum(llcont(hurdle1)), as.numeric(logLik(hurdle1))) expect_equal(sum(llcont(hurdle2)), as.numeric(logLik(hurdle2))) expect_equal(sum(llcont(hurdle3)), as.numeric(logLik(hurdle3))) expect_equal(sum(llcont(hurdle4)), as.numeric(logLik(hurdle4))) expect_equal(sum(llcont(hurdle5)), as.numeric(logLik(hurdle5))) } }) test_that("zeroinfl object", { if (isTRUE(require("pscl"))) { zi1 <- zeroinfl(art ~ . | 1, data = bioChemists) zi2 <- zeroinfl(art ~ . | 1, data = bioChemists, dist = "negbin") zi3 <- zeroinfl(art ~ . | ., data = bioChemists) zi4 <- zeroinfl(art ~ . | ., data = bioChemists, dist = "negbin") expect_equal(sum(llcont(zi1)), as.numeric(logLik(zi1))) expect_equal(sum(llcont(zi2)), as.numeric(logLik(zi2))) expect_equal(sum(llcont(zi3)), as.numeric(logLik(zi3))) expect_equal(sum(llcont(zi4)), as.numeric(logLik(zi4))) } }) test_that("lm object", { ctl <- c(4.17,5.58,5.18,6.11,4.50,4.61,5.17,4.53,5.33,5.14) trt <- c(4.81,4.17,4.41,3.59,5.87,3.83,6.03,4.89,4.32,4.69) group <- gl(2, 10, 20, labels = c("Ctl","Trt")) weight <- c(ctl, trt) lm1 <- lm(weight ~ group) expect_equal(sum(llcont(lm1)), as.numeric(logLik(lm1))) }) test_that("mlogit object", { if (isTRUE(require("mlogit")) & isTRUE(require("AER"))) { data("Fishing", package = "mlogit") Fish <- mlogit.data(Fishing, varying = c(2:9), shape = "wide", choice = "mode") mlog1 <- mlogit(mode ~ price + catch, data = Fish) mlog2 <- mlogit(mode ~ 0 | income, data = Fish) mlog3 <- mlogit(mode ~ price+ catch | income, data = Fish) mlog4 <- mlogit(mode ~ price+ catch | income, data = Fish, reflevel = "charter") mlog5 <- mlogit(mode ~ price+ catch | income, data = Fish, alt.subset = c("charter", "pier", "beach")) Fishing2 <- Fishing Fishing2[1, "price.pier"] <- Fishing2[3, "price.beach"] <- NA mlog6 <- mlogit(mode~price+catch|income, Fishing2, shape="wide", choice="mode", varying = 2:9) data("TravelMode", package = "AER") Tr2 <- TravelMode[-c(2, 7, 9),] mlog7 <- mlogit(choice~wait+gcost|income+size, Tr2, shape = "long", chid.var = "individual", alt.var="mode", choice = "choice") data("TravelMode", package = "AER") mlog8 <- mlogit(choice ~ wait + travel + vcost, TravelMode, shape = "long", chid.var = "individual", alt.var = "mode", method = "bfgs", heterosc = TRUE, tol = 10) ## fails with "TravelModel object cannot be found"; no idea why ## TravelMode$avincome <- with(TravelMode, income * (mode == "air")) ## TravelMode$time <- with(TravelMode, travel + wait)/60 ## TravelMode$timeair <- with(TravelMode, time * I(mode == "air")) ## TravelMode$income <- with(TravelMode, income / 10) ## TravelMode$incomeother <- with(TravelMode, ## ifelse(mode %in% c('air', 'car'), ## income, 0)) ## mlog9 <- mlogit(choice~gcost+wait+incomeother, TravelMode, ## shape='long', alt.var='mode', ## nests=list(public=c('train', 'bus'), ## other=c('car','air'))) data("Game", package = "mlogit") mlog10 <- mlogit(ch~own|hours, Game, choice='ch', varying = 1:12, ranked=TRUE, shape="wide", reflevel="PC") expect_equal(sum(llcont(mlog1)), as.numeric(logLik(mlog1))) expect_equal(sum(llcont(mlog2)), as.numeric(logLik(mlog2))) expect_equal(sum(llcont(mlog3)), as.numeric(logLik(mlog3))) expect_equal(sum(llcont(mlog4)), as.numeric(logLik(mlog4))) expect_equal(sum(llcont(mlog5)), as.numeric(logLik(mlog5))) expect_equal(sum(llcont(mlog6)), as.numeric(logLik(mlog6))) expect_equal(sum(llcont(mlog7)), as.numeric(logLik(mlog7))) expect_equal(sum(llcont(mlog8)), as.numeric(logLik(mlog8))) ##expect_equal(sum(llcont(mlog9)), as.numeric(logLik(mlog9))) expect_equal(sum(llcont(mlog10)), as.numeric(logLik(mlog10))) } }) test_that("nls object", { DNase1 <- subset(DNase, Run == 1) nls1 <- nls(density ~ SSlogis(log(conc), Asym, xmid, scal), DNase1) nls2 <- nls(density ~ 1/(1 + exp((xmid - log(conc))/scal)), data = DNase1, start = list(xmid = 0, scal = 1), algorithm = "plinear") nls3 <- nls(density ~ Asym/(1 + exp((xmid - log(conc))/scal)), data = DNase1, start = list(Asym = 3, xmid = 0, scal = 1)) nls4 <- nls(density ~ Asym/(1 + exp((xmid - log(conc))/scal)), data = DNase1, start = list(Asym = 3, xmid = 0, scal = 1), algorithm = "port") Treated <- Puromycin[Puromycin$state == "treated", ] weighted.MM <- function(resp, conc, Vm, K) { pred <- (Vm * conc)/(K + conc) (resp - pred) / sqrt(pred) } nls5 <- nls( ~ weighted.MM(rate, conc, Vm, K), data = Treated, start = list(Vm = 200, K = 0.1)) lisTreat <- with(Treated, list(conc1 = conc[1], conc.1 = conc[-1], rate = rate)) weighted.MM1 <- function(resp, conc1, conc.1, Vm, K) { conc <- c(conc1, conc.1) pred <- (Vm * conc)/(K + conc) (resp - pred) / sqrt(pred) } nls6 <- nls( ~ weighted.MM1(rate, conc1, conc.1, Vm, K), data = lisTreat, start = list(Vm = 200, K = 0.1)) weighted.MM.grad <- function(resp, conc1, conc.1, Vm, K) { conc <- c(conc1, conc.1) K.conc <- K+conc dy.dV <- conc/K.conc dy.dK <- -Vm*dy.dV/K.conc pred <- Vm*dy.dV pred.5 <- sqrt(pred) dev <- (resp - pred) / pred.5 Ddev <- -0.5*(resp+pred)/(pred.5*pred) attr(dev, "gradient") <- Ddev * cbind(Vm = dy.dV, K = dy.dK) dev } nls7 <- nls( ~ weighted.MM.grad(rate, conc1, conc.1, Vm, K), data = lisTreat, start = list(Vm = 200, K = 0.1)) if(isTRUE(require("MASS"))){ utils::data(muscle, package = "MASS") nls9 <- nls(Length ~ cbind(1, exp(-Conc/th)), muscle, start = list(th = 1), algorithm = "plinear") b <- coef(nls9) nls10 <- nls(Length ~ a[Strip] + b[Strip]*exp(-Conc/th), muscle, start = list(a = rep(b[2], 21), b = rep(b[3], 21), th = b[1])) expect_equal(sum(llcont(nls9)), as.numeric(logLik(nls9))) expect_equal(sum(llcont(nls10)), as.numeric(logLik(nls10))) } expect_equal(sum(llcont(nls1)), as.numeric(logLik(nls1))) expect_equal(sum(llcont(nls2)), as.numeric(logLik(nls2))) expect_equal(sum(llcont(nls3)), as.numeric(logLik(nls3))) expect_equal(sum(llcont(nls4)), as.numeric(logLik(nls4))) expect_equal(sum(llcont(nls5)), as.numeric(logLik(nls5))) expect_equal(sum(llcont(nls6)), as.numeric(logLik(nls6))) expect_equal(sum(llcont(nls7)), as.numeric(logLik(nls7))) }) test_that("polr object", { if (isTRUE(require("MASS"))) { options(contrasts = c("contr.treatment", "contr.poly")) polr1 <- polr(Sat ~ Infl + Type + Cont, weights = Freq, data = housing) polr2 <- update(polr1, method = "probit", Hess = TRUE) polr3 <- update(polr1, method = "loglog", Hess = TRUE) polr4 <- update(polr1, method = "cloglog", Hess = TRUE) expect_equal(sum(llcont(polr1)), as.numeric(logLik(polr1))) expect_equal(sum(llcont(polr2)), as.numeric(logLik(polr2))) expect_equal(sum(llcont(polr3)), as.numeric(logLik(polr3))) expect_equal(sum(llcont(polr4)), as.numeric(logLik(polr4))) } }) test_that("rlm object", { if (isTRUE(require("MASS"))) { rlm1 <- rlm(stack.loss ~ ., stackloss) rlm2 <- rlm(stack.loss ~ ., stackloss, psi = psi.hampel, init = "lts") rlm3 <- rlm(stack.loss ~ ., stackloss, psi = psi.bisquare) expect_equal(sum(llcont(rlm1)), as.numeric(logLik(rlm1))) expect_equal(sum(llcont(rlm2)), as.numeric(logLik(rlm2))) expect_equal(sum(llcont(rlm3)), as.numeric(logLik(rlm3))) } }) test_that("OpenMx object", { if (isTRUE(require("OpenMx")) & isTRUE(require("tidySEM"))) { res <- mx_lca(data = data_mix_ordinal, classes = 1:2, run = FALSE) res[[1]]$fitfunction$rowDiagnostics <- TRUE res[[2]]$class1$fitfunction$rowDiagnostics <- TRUE res[[2]]$class2$fitfunction$rowDiagnostics <- TRUE res[[1]] <- mxTryHardOrdinal(res[[1]]) res[[2]] <- mxTryHardOrdinal(res[[2]]) expect_equal(sum(llcont(res[[1]])), as.numeric(logLik(res[[1]]))) expect_equal(sum(llcont(res[[2]])), as.numeric(logLik(res[[2]]))) } }) # test_that("SingleGroupClass", { # if (isTRUE(require("mirt"))) { # data <- expand.table(LSAT7) # mirts1 <- mirt(data, 1) # mirts2 <- mirt(data, 1, SE = TRUE) # mirts3 <- mirt(data, 1, SE = TRUE, SE.type = 'SEM') # mirts4 <- mirt(data, 1, SE = TRUE, SE.type = 'BL') # mirts5 <- mirt(data, 1, itemtype = c('2PL', '2PL', '2PL', '2PL', '3PL')) # model <- 'F = 1-5 # PRIOR = (5, g, norm, -1.5, 3)' # mirts6 <- mirt(data, model, itemtype = c('2PL', '2PL', '2PL', '2PL', '3PL')) # mirts7 <- mirt(data, 2) # # expect_equal(sum(llcont(mirts1)), as.numeric(mirts1@logLik)) # expect_equal(sum(llcont(mirts2)), as.numeric(mirts2@logLik)) # expect_equal(sum(llcont(mirts3)), as.numeric(mirts3@logLik)) # expect_equal(sum(llcont(mirts4)), as.numeric(mirts4@logLik)) # expect_equal(sum(llcont(mirts5)), as.numeric(mirts5@logLik)) # expect_equal(sum(llcont(mirts6)), as.numeric(mirts6@logLik)) # expect_equal(sum(llcont(mirts7)), as.numeric(mirts7@logLik)) # } # })