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))
# }
# })
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.