Nothing
tests/testthat/test-1-main.R
In WeMix: Weighted Mixed-Effects Models Using Multilevel Pseudo Maximum Likelihood Estimation
require(testthat)
# for sleepstudy, calls in here to lmer
library(lme4)
# for grad function
options(width = 500)
options(useFancyQuotes = FALSE)
### Data Used: sleepstudy
tolerance <- 1E-3
data("sleepstudy")
### Unweigted =====================================
sleepstudyU <- sleepstudy
sleepstudyU$weight1L1 <- 1
sleepstudyU$weight1L2 <- 1
context("The model runs")
test_that("The model runs", {
wm0 <- mix(Reaction ~ Days + (1|Subject), data=sleepstudyU, weights=c("weight1L1", "weight1L2"))
lme1 <- lmer(Reaction ~ Days + (1|Subject), data=sleepstudyU, REML=FALSE)
expect_equal(wm0$lnl, -897.039321502613, tolerance=tolerance*897) # value from lmer(Reaction ~ Days + (1 | Subject), data=sleepstudy, REML=FALSE)
expect_equal(unname(c(wm0$coef)), unname(fixef(lme1)), tolerance=tolerance)
expect_equal(unname(wm0$theta), unname(lme1@theta), tolerance=tolerance)
expect_equal(coef(wm0),
expected = getME(lme1, "fixef"),
tolerance = tolerance)
lmevars1 <- data.frame(summary(lme1)$varcor)$sdcor
expect_equal(unname(wm0$vars),
expected = unname(lmevars1)^2,
tolerance = tolerance)
# agrees with GLLAMM
#source: logit_command.do
gllamm_model1 <- c("(Intercept)" = 251.4051,
"Days" = 10.46729,
"Residual" = 1296.8768,
"Subject" = 954.52789,
"lnl" = -897.03932)
expect_equal(unname(coef(wm0)),
expected = unname(gllamm_model1[1:2]),
tolerance = abs(tolerance))
expect_equal(unname(wm0$vars),
expected = unname(gllamm_model1[3:4]),
tolerance = tolerance)
expect_equal(unname(wm0$lnl),
expected=unname(gllamm_model1[5]),
tolerance=abs(tolerance))
sleepstudy2 <- sleepstudy
sleepstudy2$block <- rep(1:6,each=30)
sleepstudy2$weight1L1 <- sleepstudy2$weight1L2 <- sleepstudy2$weight1L3 <- 1
# change Reaction by the artifical block we just added
sleepstudy2$Reaction <- sleepstudy2$Reaction + sleepstudy2$block*20
names(sleepstudy2) <- gsub("weight1L","pwt",names(sleepstudy2))
lmr <- lmer(Reaction ~ Days + (1|Subject) + (1|block), data=sleepstudy2, REML=FALSE)
mm <- mix(Reaction ~ Days + (1|Subject) + (1|block), data=sleepstudy2, weights=c("pwt1", "pwt2", "pwt3"))
expect_equal(unname(mm$lnl),
expected=unname(logLik(lmr)[[1]]),
tolerance=tolerance)
})
context("Factor binomial")
test_that("Factor binomial", {
sleepstudyM <- sleepstudyU
sleepstudyM$highR <- factor(ifelse(sleepstudyM$Reaction > 340, 2, 1), levels=c(1,2), labels=c("L","H"))
sleepstudyM$Sub <- factor(sleepstudyM$Subject, levels=300:400)
m1 <- mix(highR ~ Days + (1|Subject), data=sleepstudyM, weights=c("weight1L1", "weight1L2"),
family="binomial")
expect_equal(coef(m1), c(`(Intercept)` = -8.69336570860094, Days = 1.17236933934874), tol=1e-3)
summaryREF <- c("Call:",
"mix(formula = highR ~ Days + (1 | Subject), data = sleepstudyM, ",
" weights = c(\"weight1L1\", \"weight1L2\"), family = \"binomial\")",
"",
"Variance terms:",
" Level Group Name Variance Std. Error Std.Dev.",
" 2 Subject (Intercept) 7.99 5.14 2.83", "Groups:",
" Level Group n size mean wgt sum wgt",
" 2 Subject 18 1 18",
" 1 Obs 180 1 180",
"",
"Fixed Effects:",
" Estimate Std. Error t value",
"(Intercept) -8.694 2.073 -4.19",
"Days 1.172 0.278 4.21",
"",
"lnl= -53.90 ")
withr::with_options(list(digits=2),
co <- capture.output(summary(m1))
)
expect_equal(co, summaryREF)
})
context("Unweighted model with 2 random effects")
test_that("Agrees with lme4 3,handles missing data", {
sleepstudyM <- sleepstudyU
sleepstudyM <- sleepstudyM[order(runif(nrow(sleepstudyM))), ]
#introduce a missing value
sleepstudyM$Days[3] <- NA
sleepstudyM$Subject <- ifelse(sleepstudyM$Subject %in% c("310"),"2",as.character(sleepstudyM$Subject))
sleepstudyM$Subject <- ifelse(sleepstudyM$Subject %in% c("309"),"1",as.character(sleepstudyM$Subject))
sleepstudyM$Subject <- ifelse(sleepstudyM$Subject %in% c("330"),"3",as.character(sleepstudyM$Subject))
lme2 <- lmer(Reaction ~ Days + (1 | Subject) + (0 + Days | Subject), data=sleepstudyM, REML=FALSE)
# the dropped row should cause a warning
expect_warning(wm2 <- mix(Reaction ~ Days + (1 | Subject) + (0 + Days | Subject), data=sleepstudyM, weights=c("weight1L1", "weight1L2")),"with missing data")
expect_equal(wm2$lnl, as.numeric(logLik(lme2)), tol=2E-7*abs(as.numeric(logLik(lme2))))
# check coef
expect_equal(coef(wm2),
expected = getME(lme2, "fixef"),
tolerance = tolerance)
# check vars
lmewm2vars <- data.frame(summary(lme2)$varcor)$sdcor
expect_equal(unname(wm2$vars),
expected = unname(lmewm2vars)^2,
tolerance = tolerance)
lme4ranef <- ranef(lme2)$Subject
# not in WeMix output
attr(lme4ranef, "postVar") <- NULL
expect_equal(lme4ranef, wm2$ranefMat$Subject, 0.01)
expect_equal(predict(lme2), predict(wm2), tol=1e-4)
sleepstudyM2 <- sleepstudyM[order(runif(nrow(sleepstudyM))), ]
sleepstudyM2$Subject <- ifelse(sleepstudyM2$Subject %in% c("332"),"4",as.character(sleepstudyM$Subject))
sleepstudyM2$Subject <- ifelse(sleepstudyM2$Subject %in% c("308"),"1",as.character(sleepstudyM$Subject))
sleepstudyM2$Subject <- ifelse(sleepstudyM2$Subject %in% c("351"),"5",as.character(sleepstudyM$Subject))
expect_equal(predict(lme2,sleepstudyM2, allow.new.levels=TRUE), predict(wm2,sleepstudyM2, allow.new.levels=TRUE), tol=1e-4)
})
context("Mean Centering Matches HLM results")
test_that("Mean Centering Matches HLM results", {
wm1 <- mix(Reaction ~ Days + (1|Subject), data=sleepstudyU, weights=c("weight1L1", "weight1L2"), nQuad=13,center_group=list("Subject"= as.formula(~Days)), verbose=FALSE)
expect_equal(unname(wm1$coef), c(298.507892, 10.467286), tolerance=1E-1)
expect_equal(wm1$lnl, -897.0393, tolerance=tolerance*897) # value from lmer(Reaction ~ Days + (1 | Subject), data=sleepstudy, REML=FALSE)
})
ss <- sleepstudy
ss1 <- ss
ss2 <- ss
doubles <- c(308, 309, 310) # subject with double obs
ss2 <- rbind(ss, subset(ss, Subject %in% doubles))
ss$w <- ifelse(ss$Subject %in% doubles, 2, 1)
contrasts(ss2$Subject) <- "contr.sum"
ss1$W1 <- ifelse(ss1$Subject %in% doubles, 2, 1)
ss1$W2 <- 1
ss1$bin <- ifelse(sleepstudy$Reaction<300,0,1) #for the binomial test
ss2$W2 <- ss2$W1 <- 1
doubles <- c(308, 309, 310) # subject with double obs
ss30 <- subset(ss, Subject %in% doubles)
ss30$Subject <- as.numeric(as.character(ss30$Subject)) + 1000
ss0 <- ss
ss0$Subject <- as.numeric(as.character(ss$Subject))
ss3 <- rbind(ss0, ss30)
ss3$Subject <- as.factor(ss3$Subject)
ss3$W2 <- 1
ss3$W1 <- 1
ss4 <- ss
ss4$W2 <- ifelse(ss4$Subject %in% doubles, 2, 1)
ss4$W1 <- ifelse(ss4$Subject %in% doubles, 2, 1) # make unconditional weight
context("GLM works: Binomial")
test_that("GLM works: Binomial", {
skip_on_cran()
#full test for binomial
bi_1 <- mix(bin~Days + (1|Subject), data=ss1, family=binomial(link="logit"), verbose=FALSE,
weights=c("W1", "W2"), nQuad=13)
expect_equal(unname(bi_1$coef), c(-3.3448,.5928), tolerance=1E-3)
expect_equal(as.numeric(bi_1$lnl), -93.751679, tolerance=1E-5)
sum_bi <- summary(bi_1)
expect_is(summary(bi_1), "summaryWeMixResults")
})
context("Repeating is the same as weighting: L1 replicate vs weighting")
test_that("Repeating is the same as weighting: L1 replicate vs weighting", {
# mix for L1, weighted
wmeL1W <- mix(formula=Reaction ~ Days + (1 | Subject), data=ss1,
weights=c("W1", "W2"))
# mix for L1, duplicated
system.time(wmeL1D <- mix(formula=Reaction ~ Days + (1 | Subject), data=ss2,
weights=c("W1", "W2")))
# check weighted agrees with duplicated lme4 results
expect_equal(wmeL1W$lnl, -1048.34318418762, tolerance=1050*tolerance)
# check duplicated agrees with duplicated lme4 results
expect_equal(wmeL1D$lnl, -1048.34318418762, tolerance=1050*tolerance)
})
context("grouping factor not sorted")
test_that("grouping factor not sorted", {
skip_on_cran()
ss1_mixed <- ss1[c(125:180,1,100,2:99,101:124),]
row.names(ss1_mixed) <- NULL
# mix for L1, weighted
wmeL1W <- mix(formula=Reaction ~ Days + (1 | Subject), data=ss1_mixed,
weights=c("W1", "W2"), nQuad=13, run=TRUE, verbose=FALSE)
# mix for L1, duplicated
system.time(wmeL1D <- mix(formula=Reaction ~ Days + (1 | Subject), data=ss2,
weights=c("W1", "W2"), nQuad=13, run=FALSE, verbose=FALSE))
#statares <- c(251.4619, 10.40726, 1000.7466, 1338.0865) # not used
# check weighted agrees with duplicated lme4 results
expect_equal(wmeL1W$lnl, -1048.34318418762, tolerance=tolerance)
expect_equal(wmeL1D$lnl, -1048.34318418762, tolerance=tolerance)
# check final results
suppressWarnings(mix1 <- mix(formula=Reaction ~ Days + (1 | Subject), data=ss1_mixed,
weights=c("W1", "W2")))
suppressWarnings(mix1REF <- mix(formula=Reaction ~ Days + (1 | Subject), data=ss1,
weights=c("W1", "W2")))
expect_equal(mix1$coef, mix1REF$coef, tolerance=1e3)
expect_equal(mix1$vars, mix1REF$vars, tolerance=1e-3)
expect_equal(mix1$lnl, mix1REF$lnl, tolerance=1e-3)
#check weighted fixed effects variances
expect_equal(unname(sqrt(diag(mix1$cov_mat))), unname(sqrt(diag(mix1REF$cov_mat))),tolerance = tolerance)
})
context("Weighted three level model unsorted")
test_that("Weighted three level model unsorted", {
skip_on_cran()
sleepstudy2 <- sleepstudy
sleepstudy2$Group <- 1
sleepstudy2$Group <- ifelse(sleepstudy2$Subject %in% c("310", "309", "349", "335"), 2, sleepstudy2$Group)
sleepstudy2$Group <- ifelse(sleepstudy2$Subject %in% c("330", "352", "337", "369"), 3, sleepstudy2$Group)
sleepstudy2$Group <- ifelse(sleepstudy2$Subject %in% c("331", "332", "334", "372"), 4, sleepstudy2$Group)
sleepstudy2$Group <- factor(sleepstudy2$Group)
ss2 <- sleepstudy2
w1c <- w1 <- rep(1,180)
w2c <- w2 <- rep(1,18)
w3c <- w3 <- rep(1,4)
# unbalanced (non-identical within a Subject), level-1 (obs level) weights
w1c[1:4] <- w1[1:4] <- 2
uR <- sleepstudy2[1:4,]
sleepstudy2 <- rbind(sleepstudy2, uR)
# level-2 weights
w2c[1] <- w2[1] <- 2
w1[ss2$Subject=="308"] <- 2*w1[ss2$Subject=="308"]
sR <- subset(sleepstudy2, Subject == "308")
sR$Subject <- "S308"
sleepstudy2 <- rbind(sleepstudy2, sR)
# level-3 weights
w3c[1] <- w3[1] <- 2
w2[c(1,7,12,13,16,17)] <- 2*w2[c(1,7,12,13,16,17)]
w1[ss2$Group==1] <- 2*w1[ss2$Group==1]
gR <- subset(sleepstudy2, Group %in% 1)
gR$Subject <- paste0("G", gR$Subject)
gR$Group <- "11"
sleepstudy2 <- rbind(sleepstudy2, gR)
# lmr for reference
lmr <- lmer(Reaction ~ Days + (1|Subject) + (1|Group),
data=sleepstudy2, control=lmerControl(optimizer="bobyqa"),
REML=FALSE)
ss2$w1 <- w1
ss2$w2 <- rep(w2,each=10)
ss2$w3 <- ifelse(ss2$Subject %in% c("308", "333", "350", "351", "370", "371"),2,1)
ss3 <- ss2[sample(row.names(ss2), size=nrow(ss2)), ]
wm0 <- mix(Reaction ~ Days + (1|Subject) + (1|Group), data=ss3,
weights=c("w1", "w2", "w3"))
# check coef
expect_equal(coef(wm0),
expected = getME(lmr, "fixef"),
tolerance = tolerance)
lmeRanef <- ranef(lmr)$Subject
# not in WeMix output
attr(lmeRanef, "postVar") <- NULL
# only first 18 of lmeRanef are unique/comparable
expect_equal(lmeRanef[1:18,,drop=FALSE], wm0$ranefMat$Subject, 20*sqrt(.Machine$double.eps))
# check vars
lmewm2vars <- data.frame(summary(lmr)$varcor)$sdcor
expect_equal(unname(wm0$vars),
expected = unname(lmewm2vars)^2,
tolerance = tolerance)
# var beta not expected to be equal
})
context("repeating is the same as weighting: L2 replicate vs weighting")
test_that("L2 replicate vs weighting", {
# mix for L2, duplicated
system.time(wmeL2D <- mix(formula=Reaction ~ Days + (1 | Subject),
data=ss3, weights=c("W1", "W2")))
# mix for L2, weighted
system.time(wmeL2W <- mix(formula=Reaction ~ Days + (1 | Subject), data=ss4,
weights=c("W1", "W2")))
expect_equal(wmeL2W$lnl, -1055.34690957995, tolerance=2E-7)
expect_equal(wmeL2D$lnl, -1055.34690957995, tolerance=2E-7)
})
context("Repeating is the same as weighting: L1 replicate vs weighting, 2 REs")
test_that("Repeating is the same as weighting: L1 replicate vs weighting, 2 REs", {
# mix for L1, weighted, 2 REs
wmeL1WRE2 <- mix(formula=Reaction ~ Days + (1 | Subject) + (0+Days|Subject),
data=ss1, weights=c("W1", "W2"), nQuad=13, run=FALSE, verbose=FALSE)
# mix for L1, duplicated, 2 REs
wmeL1DRE2 <- mix(formula=Reaction ~ Days + (1 | Subject) + (0+Days|Subject),
data=ss2, weights=c("W1", "W2"),nQuad=13, run=FALSE, verbose=FALSE)
expect_equal(wmeL1WRE2$lnl, -1018.29298875158, tolerance=1050*2E-7)
expect_equal(wmeL1DRE2$lnl, -1018.29298875158, tolerance=1050*2E-7)
})
ssB <- sleepstudy
set.seed(2)
ssB$Reaction <- ssB$Days * 3.141 + rnorm(nrow(ssB))
ssB$W2 <- 1
ssB$W1 <- 1
context("Zero variance estimate")
test_that("Zero variance estimate", {
skip_on_cran()
# this has 0 variance estimate in lmer
lmeB <- lmer(Reaction ~ Days + (1|Subject), data=ssB, REML=FALSE)
mixB <- mix(formula=Reaction ~ Days + (1 | Subject), data=ssB,
weights=c("W1", "W2"), nQuad=13, run=TRUE, verbose=FALSE)
expect_equal(mixB$lnl, as.numeric(logLik(lmeB)), tol=1e-5)
expect_equal(coef(mixB), fixef(lmeB), tol=1e-5)
expect_equal(unname(mixB$vars[length(mixB$vars)]), unname(lmeB@devcomp$cmp["sigmaML"]^2), tol=1e-5)
ss1 <- sleepstudy
#add group variables for 3 level model
ss1$Group <- 1
ss1$Group <- ifelse(ss1$Subject %in% c(349,335,330, 352, 337, 369), 2, ss1$Group)
# Create weights
ss1$W1 <- ifelse(ss1$Subject %in% c(308, 309, 310), 2, 1)
ss1$W2 <- 1
ss1$W3 <- ifelse(ss1$Group == 2, 2, 1)
#Run three level model with random slope and intercept.
three_level <- mix(Reaction ~ Days + (1|Subject) + (1+Days|Group), data=ss1, weights = c("W1","W2","W3"))
expect_is(three_level, "WeMixResults")
})
context("Unweighted three level model")
test_that("Unweighted three level model", {
sleepstudy2 <- sleepstudy
sleepstudy2$Group <- 1
sleepstudy2$Group <- ifelse(sleepstudy2$Subject %in% c("310", "309", "349", "335"), 2, sleepstudy2$Group)
sleepstudy2$Group <- ifelse(sleepstudy2$Subject %in% c("330", "352", "337", "369"), 3, sleepstudy2$Group)
sleepstudy2$Group <- ifelse(sleepstudy2$Subject %in% c("331", "332", "334", "372"), 4, sleepstudy2$Group)
sleepstudy2$Group <- factor(sleepstudy2$Group)
sleepstudy2$w1 <- 1
sleepstudy2$w2 <- 1
sleepstudy2$w3 <- 1
wm0 <- mix(Reaction ~ Days + (1|Subject) + (0+Days|Subject) + (1 | Group), data=sleepstudy2, weights=c("w1", "w2","w3"), verbose=FALSE, run=TRUE)
lm0 <- lmer(Reaction ~ Days + (1|Subject) + (0+Days|Subject) + (1 | Group), data=sleepstudy2,REML=FALSE)
wm0RE <- wm0$ranefMat
lm0RE <- ranef(lm0)
expect_equal(names(wm0RE), names(lm0RE))
attr(lm0RE$Subject, "postVar") <- NULL
expect_equal(wm0RE$Subject, lm0RE$Subject, tol=.Machine$double.eps^0.25)
attr(lm0RE$Group, "postVar") <- NULL
expect_equal(wm0RE$Group, lm0RE$Group, tol=.Machine$double.eps^0.25)
wm0RE <- wm0$ranefMat
lm0RE <- ranef(lm0)
attr(lm0RE$Subject, "postVar") <- NULL
expect_equal(wm0RE$Subject, lm0RE$Subject, tol=.Machine$double.eps^0.25)
# check vars
lmevars1 <- data.frame(summary(lm0)$varcor)$sdcor
expect_equal(unname(wm0$vars),
expected = unname(lmevars1)^2,
tolerance = 1e-3)
expect_equal(wm0$lnl, as.numeric(logLik(lm0)), tol=1e-3)
expect_equal(coef(wm0), fixef(lm0), tol=1e-4)
expect_equal(unname(wm0$vars[length(wm0$vars)]), unname(lm0@devcomp$cmp["sigmaML"]^2), tol=1e-4)
ranef <- ranef(lm0)
# WeMix does not have postVar attribute
attr(ranef$Subject, "postVar") <- NULL
attr(ranef$Group, "postVar") <- NULL
expect_equal(ranef$Subject, wm0$ranefMat$Subject, (.Machine$double.eps)^0.25)
expect_equal(ranef$Group, wm0$ranefMat$Group, (.Machine$double.eps)^0.25)
})
context("Unweighted three level model, binomial")
test_that("Unweighted three level model, binomial", {
sleepstudy2 <- sleepstudy
set.seed(564)
sleepstudy2$rand_slope <- rnorm(180,mean=3.25,sd=1.79)
sleepstudy2$highR <- ifelse(sleepstudy2$Reaction > 340, 1, 0)
sleepstudy2$Group <- 1
sleepstudy2$Group <- ifelse(sleepstudy2$Subject %in% c("310", "309", "349", "335"), 2, sleepstudy2$Group)
sleepstudy2$Group <- ifelse(sleepstudy2$Subject %in% c("330", "352", "337", "369"), 3, sleepstudy2$Group)
sleepstudy2$Group <- ifelse(sleepstudy2$Subject %in% c("331", "332", "334", "372"), 4, sleepstudy2$Group)
sleepstudy2$Group <- factor(sleepstudy2$Group)
sleepstudy2$w1 <- 1
sleepstudy2$w2 <- 1
sleepstudy2$w3 <- 1
sleepstudy2 <- sleepstudy2[sample(1:nrow(sleepstudy2), nrow(sleepstudy2)),]
# up the max_iterations for more complex model structure
suppressWarnings(wm0 <- mix(highR ~ Days + (rand_slope|Subject) + (1 | Group), data=sleepstudy2,
family=binomial(),weights=c("w1", "w2","w3"), verbose=FALSE, run=TRUE,
nQuad=1,max_iteration = 20))
lm0 <- glmer(highR ~ Days + (rand_slope|Subject) + (1 | Group), data=sleepstudy2,
family=binomial())
# check vars
lmevars1 <- data.frame(summary(lm0)$varcor)$vcov
expect_equal(unname(wm0$vars),
expected = unname(lmevars1),
tolerance = 1e-3)
expect_equal(as.numeric(wm0$lnl), as.numeric(logLik(lm0)), tol=1e-3)
expect_equal(coef(wm0), fixef(lm0), tol=1e-4)
expect_equal(predict(wm0), predict(lm0), tol=1e-4)
expect_equal(predict(wm0, sleepstudy2), predict(lm0, sleepstudy2), tol=1e-4)
sleepstudy3 <- sleepstudy2[sample(1:nrow(sleepstudy2), floor(nrow(sleepstudy2)/4)),]
expect_equal(a1 <- predict(wm0, sleepstudy3), b1 <- predict(lm0, sleepstudy3), tol=1e-4)
sleepstudy4 <- sleepstudy2
sleepstudy4$Subject <- ifelse(sleepstudy4$Subject == "332", "1", as.character(sleepstudy4$Subject))
expect_equal(a2 <- predict(wm0, sleepstudy4, allow.new.levels=TRUE), b2 <- predict(lm0, sleepstudy4, allow.new.levels=TRUE), tol=1e-4)
})
context("Weighted v unweighted replicated two level model")
test_that("Weighted v unweighted replicated two level model", {
sleepstudy2 <- sleepstudy
sleepstudy2$w1 <- 1
sleepstudy2$w2 <- 1
sleepstudy2$w3 <- 1
sleepstudy2$Subject <- as.character(sleepstudy2$Subject)
g1 <- subset(sleepstudy2, sleepstudy2$Subject %in% c("310", "309", "349", "335"))
g1$Subject <- paste0("R", g1$Subject)
sleepstudyrep <- rbind(sleepstudy2, g1)
sleepstudy2$Subject <- factor(sleepstudy2$Subject)
sleepstudy2$w2 <- ifelse(sleepstudy2$Subject %in% c("310", "309", "349", "335"), 2, 1)
# two level
wm0 <- mix(Reaction ~ Days + (1|Subject), data=sleepstudyrep, weights=c("w1", "w2"), verbose=FALSE, run=TRUE)
wmw <- mix(Reaction ~ Days + (1|Subject), data=sleepstudy2, weights=c("w1", "w2"), verbose=FALSE, run=TRUE, cWeights=TRUE)
lm0 <- lmer(Reaction ~ Days + (1|Subject), data=sleepstudyrep, REML=FALSE)
# check vars
lmevars1 <- data.frame(summary(lm0)$varcor)$sdcor
expect_equal(unname(wm0$vars),
expected = unname(lmevars1)^2,
tolerance = 1e-3)
expect_equal(wm0$lnl, as.numeric(logLik(lm0)), tol=1e-3)
expect_equal(coef(wm0), fixef(lm0), tol=1e-4)
expect_equal(unname(wm0$vars[length(wm0$vars)]), unname(lm0@devcomp$cmp["sigmaML"]^2), tol=1e-4)
lm0ranef <- ranef(lm0)$Subject
# not in WeMix output
attr(lm0ranef, "postVar") <- NULL
expect_equal(lm0ranef, wm0$ranefMat$Subject, 20*sqrt(.Machine$double.eps))
expect_equal(wmw$lnl, wm0$lnl, tol=1e-3)
expect_equal(coef(wmw), coef(wm0), tol=1e-4)
expect_equal(wmw$vars, wm0$vars, tol=1e-4)
# should be the same if not conditional as well
sleepstudy2$w1 <- sleepstudy2$w2
wmw <- mix(Reaction ~ Days + (1|Subject), data=sleepstudy2, weights=c("w1", "w2"), verbose=FALSE, run=TRUE, cWeights=FALSE)
expect_equal(wmw$lnl, wm0$lnl, tol=1e-3)
expect_equal(coef(wmw), coef(wm0), tol=1e-4)
expect_equal(wmw$vars, wm0$vars, tol=1e-4)
# sub 1 conditional weights work
sleepstudy2 <- sleepstudy
sleepstudy2$w1 <- 1
sleepstudy2$w2 <- 1
sleepstudy2$w3 <- 1
sleepstudy2$Subject <- as.character(sleepstudy2$Subject)
g1 <- subset(sleepstudy2, sleepstudy2$Subject %in% c("310", "309", "349", "335"))
# do not renaem subject
sleepstudyrep <- rbind(sleepstudy2, g1)
sleepstudyrep$w1 <- ifelse(sleepstudyrep$Subject %in% c("310", "309", "349", "335"), 0.5, 1)
# two level
wm0 <- mix(Reaction ~ Days + (1|Subject), data=sleepstudyrep, weights=c("w1", "w2"), verbose=FALSE, run=TRUE, cWeights=TRUE)
wmw <- mix(Reaction ~ Days + (1|Subject), data=sleepstudy2, weights=c("w1", "w2"), verbose=FALSE, run=TRUE)
lm0 <- lmer(Reaction ~ Days + (1|Subject), data=sleepstudy2, REML=FALSE)
# check vars
lmevars1 <- data.frame(summary(lm0)$varcor)$sdcor
expect_equal(unname(wm0$vars),
expected = unname(lmevars1)^2,
tolerance = 1e-3)
expect_equal(wm0$lnl, as.numeric(logLik(lm0)), tol=1e-3)
expect_equal(coef(wm0), fixef(lm0), tol=1e-4)
expect_equal(unname(wm0$vars[length(wm0$vars)]), unname(lm0@devcomp$cmp["sigmaML"]^2), tol=1e-4)
expect_equal(wmw$lnl, wm0$lnl, tol=1e-3)
expect_equal(coef(wmw), coef(wm0), tol=1e-4)
expect_equal(wmw$vars, wm0$vars, tol=1e-4)
# should be the same if not conditional as well
wm0 <- mix(Reaction ~ Days + (1|Subject), data=sleepstudyrep, weights=c("w1", "w2"), verbose=FALSE, run=TRUE)
expect_equal(wmw$lnl, wm0$lnl, tol=1e-3)
expect_equal(coef(wmw), coef(wm0), tol=1e-4)
expect_equal(wmw$vars, wm0$vars, tol=1e-4)
})
context("Weighted v unweighted replicated three level model")
test_that("Weighted v unweighted replicated three level model", {
sleepstudy2 <- sleepstudy
sleepstudy2$Group <- 1
sleepstudy2$Group <- ifelse(sleepstudy2$Subject %in% c("310", "309", "349", "335"), 2, sleepstudy2$Group)
sleepstudy2$Group <- ifelse(sleepstudy2$Subject %in% c("330", "352", "337", "369"), 3, sleepstudy2$Group)
sleepstudy2$Group <- ifelse(sleepstudy2$Subject %in% c("331", "332", "334", "372"), 4, sleepstudy2$Group)
#sleepstudy2$Group <- factor(sleepstudy2$Group)
sleepstudy2$w1 <- 1
sleepstudy2$w2 <- 1
sleepstudy2$w3 <- 1
sleepstudy2$Subject <- as.character(sleepstudy2$Subject)
g1 <- subset(sleepstudy2, Group == 1)
g1$Group <- 5
g1$Subject <- paste0("5", g1$Subject)
sleepstudyrep <- rbind(sleepstudy2, g1)
s308 <- subset(sleepstudy2, Subject %in% "308")
s308$Subject <- "R308"
sleepstudyrep <- rbind(sleepstudyrep, s308)
s308$Group <- 5
s308$Subject <- "5R308"
sleepstudyrep <- rbind(sleepstudyrep, s308)
sleepstudy2$w1u <- ifelse(sleepstudy2$Group == 1, 2, 1)
s308unit <- subset(sleepstudyrep, Subject %in% c("308", "5308", "R308", "5R308") & Days %in% c(0,3))
sleepstudyrep <- rbind(sleepstudyrep, s308unit)
sleepstudy2$Group <- factor(sleepstudy2$Group)
sleepstudy2$Subject <- factor(sleepstudy2$Subject)
sleepstudy2$w1 <- ifelse(sleepstudy2$Subject %in% c("308", "R308", "5R308") & sleepstudy2$Days %in% c(0,3), 2, 1)
sleepstudy2$w2 <- ifelse(sleepstudy2$Subject == "308", 2, 1)
sleepstudy2$w3 <- ifelse(sleepstudy2$Group == 1, 2, 1)
sleepstudy2$w2u <- ifelse(sleepstudy2$Group == 1, 2, 1)
sleepstudy2$w2u <- ifelse(sleepstudy2$Subject == "308", 4, sleepstudy2$w2u)
sleepstudy2$w1u <- ifelse(sleepstudy2$Group == 1, 2, 1)
sleepstudy2$w1u <- ifelse(sleepstudy2$Subject == "308", 4, sleepstudy2$w1u)
sleepstudy2$w1u <- ifelse(sleepstudy2$Subject == "308" & sleepstudy2$Days %in% c(0,3), 8, sleepstudy2$w1u)
# replication
wm0 <- mix(Reaction ~ Days + (1|Subject) + (1 | Group), data=sleepstudyrep, weights=c("w1", "w2","w3"), verbose=FALSE, run=TRUE)
# conditional weights
wmw <- mix(Reaction ~ Days + (1|Subject) + (1 | Group), data=sleepstudy2, weights=c("w1", "w2","w3"), verbose=FALSE, run=TRUE, cWeights=TRUE)
expect_equal(wm0$lnl, wmw$lnl, tol=20*sqrt(.Machine$double.eps))
expect_equal(coef(wmw), coef(wm0), tol=20*sqrt(.Machine$double.eps))
# unconditional weights
wmwc <- mix(Reaction ~ Days + (1|Subject) + (1 | Group), data=sleepstudy2, weights=c("w1u", "w2u","w3"), verbose=FALSE, run=TRUE, cWeights=FALSE)
expect_equal(wm0$lnl, wmwc$lnl, tol=20*sqrt(.Machine$double.eps))
expect_equal(coef(wmwc), coef(wm0), tol=20*sqrt(.Machine$double.eps))
lm0 <- lmer(Reaction ~ Days + (1|Subject) + (1 | Group), data=sleepstudyrep,REML=FALSE)
# check vars
lmevars1 <- data.frame(summary(lm0)$varcor)$sdcor
expect_equal(unname(wm0$vars),
expected = unname(lmevars1)^2,
tolerance = 1e-3)
expect_equal(wm0$lnl, as.numeric(logLik(lm0)), tol=1e-3)
expect_equal(coef(wm0), fixef(lm0), tol=1e-4)
expect_equal(unname(wm0$vars[length(wm0$vars)]), unname(lm0@devcomp$cmp["sigmaML"]^2), tol=1e-4)
expect_equal(wmw$lnl, wm0$lnl, tol=1e-3)
expect_equal(coef(wmw), coef(wm0), tol=1e-4)
expect_equal(wmw$vars, wm0$vars, tol=1e-4)
expect_equal(wm0$ranefMat$Subject[1:18,,drop=FALSE], wmw$ranefMat$Subject, 20*sqrt(.Machine$double.eps))
expect_equal(wm0$ranefMat$Group[1:4,,drop=FALSE], wmw$ranefMat$Group, 20*sqrt(.Machine$double.eps))
})
context("Weighted v unweighted replicated two level model, binomial")
test_that("Weighted v unweighted replicated two level model, binomial", {
sleepstudy2 <- sleepstudy
sleepstudy2$w1 <- 1
sleepstudy2$w2 <- 1
sleepstudy2$w3 <- 1
sleepstudy2$Subject <- as.character(sleepstudy2$Subject)
sleepstudy2$highR <- ifelse(sleepstudy2$Reaction > 340, 1, 0)
g1 <- subset(sleepstudy2, sleepstudy2$Subject %in% c("310", "309", "349", "335"))
g1$Subject <- paste0("R", g1$Subject)
sleepstudyrep <- rbind(sleepstudy2, g1)
sleepstudy2$Subject <- factor(sleepstudy2$Subject)
sleepstudy2$w2 <- ifelse(sleepstudy2$Subject %in% c("310", "309", "349", "335"), 2, 1)
# two level
wm0 <- mix(highR ~ Days + (1|Subject), data=sleepstudyrep, weights=c("w1", "w2"),
family=binomial(), verbose=FALSE)
wmw <- mix(highR ~ Days + (1|Subject), data=sleepstudy2, weights=c("w1", "w2"),
family=binomial(), verbose=FALSE, cWeights=TRUE,
start=c(wm0$coef,wm0$vars))
lm0 <- glmer(highR ~ Days + (1|Subject), data=sleepstudyrep, family=binomial(),
nAGQ = 13)
# check vars
lmevars1 <- data.frame(summary(lm0)$varcor)$sdcor
expect_equal(unname(wm0$vars),
expected = unname(lmevars1)^2,
tolerance = 1e-3)
expect_equal(as.numeric(wm0$lnl), as.numeric(logLik(lm0)), tol=1e-3)
expect_equal(coef(wm0), fixef(lm0), tol=1e-4)
expect_equal(wmw$lnl, wm0$lnl, tol=1e-3)
expect_equal(coef(wmw), coef(wm0), tol=1e-4)
expect_equal(wmw$vars, wm0$vars, tol=1e-4)
# should be the same if not conditional as well
sleepstudy2$w1 <- sleepstudy2$w2
wmw <- mix(highR ~ Days + (1|Subject), data=sleepstudy2, weights=c("w1", "w2"),
family=binomial(),
cWeights=FALSE)
expect_equal(wmw$lnl, wm0$lnl, tol=1e-3)
expect_equal(coef(wmw), coef(wm0), tol=1e-4)
expect_equal(wmw$vars, wm0$vars, tol=1e-4)
# sub 1 conditional weights work
sleepstudy2 <- sleepstudy
sleepstudy2$highR <- ifelse(sleepstudy2$Reaction > 340, 1, 0)
sleepstudy2$w1 <- 1
sleepstudy2$w2 <- 1
sleepstudy2$w3 <- 1
sleepstudy2$Subject <- as.character(sleepstudy2$Subject)
g1 <- subset(sleepstudy2, sleepstudy2$Subject %in% c("310", "309", "349", "335"))
# do not renaem subject
sleepstudyrep <- rbind(sleepstudy2, g1)
sleepstudyrep$w1 <- ifelse(sleepstudyrep$Subject %in% c("310", "309", "349", "335"), 0.5, 1)
# two level
wm0 <- mix(highR ~ Days + (1|Subject), data=sleepstudyrep, weights=c("w1", "w2"),
family=binomial(),verbose=FALSE, run=TRUE, cWeights=TRUE)
wmw <- mix(highR ~ Days + (1|Subject), data=sleepstudy2, weights=c("w1", "w2"),
family=binomial(),verbose=FALSE, run=TRUE)
lm0 <- glmer(highR ~ Days + (1|Subject), data=sleepstudy2, family=binomial(),
nAGQ=13)
# check vars
lmevars1 <- data.frame(summary(lm0)$varcor)$sdcor
expect_equal(unname(wm0$vars),
expected = unname(lmevars1)^2,
tolerance = 1e-3)
expect_equal(as.numeric(wm0$lnl), as.numeric(logLik(lm0)), tol=1e-3)
expect_equal(coef(wm0), fixef(lm0), tol=1e-4)
expect_equal(wmw$lnl, wm0$lnl, tol=1e-3)
expect_equal(coef(wmw), coef(wm0), tol=1e-4)
expect_equal(wmw$vars, wm0$vars, tol=1e-4)
# should be the same if not conditional as well
wm0 <- mix(highR ~ Days + (1|Subject), data=sleepstudyrep, weights=c("w1", "w2"),
family=binomial(),verbose=FALSE, run=TRUE)
expect_equal(wmw$lnl, wm0$lnl, tol=1e-3)
expect_equal(coef(wmw), coef(wm0), tol=1e-4)
expect_equal(wmw$vars, wm0$vars, tol=1e-4)
})
context("Weighted v unweighted replicated two level model, binomial with probit link")
test_that("Weighted v unweighted replicated two level model, binomial with probit link", {
sleepstudy2 <- sleepstudy
sleepstudy2$w1 <- 1
sleepstudy2$w2 <- 1
sleepstudy2$w3 <- 1
sleepstudy2$Subject <- as.character(sleepstudy2$Subject)
sleepstudy2$highR <- ifelse(sleepstudy2$Reaction > 340, 1, 0)
g1 <- subset(sleepstudy2, sleepstudy2$Subject %in% c("310", "309", "349", "335"))
g1$Subject <- paste0("R", g1$Subject)
sleepstudyrep <- rbind(sleepstudy2, g1)
sleepstudy2$Subject <- factor(sleepstudy2$Subject)
sleepstudy2$w2 <- ifelse(sleepstudy2$Subject %in% c("310", "309", "349", "335"), 2, 1)
# two level
suppressWarnings(wm0 <- mix(highR ~ Days + (1|Subject), data=sleepstudyrep, weights=c("w1", "w2"),
family=binomial(link="probit"), verbose=FALSE))
suppressWarnings(wmw <- mix(highR ~ Days + (1|Subject), data=sleepstudy2, weights=c("w1", "w2"),
family=binomial(link="probit"), verbose=FALSE, cWeights=TRUE,
start=c(wm0$coef,wm0$vars)))
lm0 <- glmer(highR ~ Days + (1|Subject), data=sleepstudyrep, family=binomial(link="probit"),
nAGQ = 13)
# check vars
lmevars1 <- data.frame(summary(lm0)$varcor)$sdcor
expect_equal(unname(wm0$vars),
expected = unname(lmevars1)^2,
tolerance = 1e-3)
expect_equal(as.numeric(wm0$lnl), as.numeric(logLik(lm0)), tol=1e-3)
expect_equal(coef(wm0), fixef(lm0), tol=1e-4)
expect_equal(wmw$lnl, wm0$lnl, tol=1e-3)
expect_equal(coef(wmw), coef(wm0), tol=1e-4)
expect_equal(wmw$vars, wm0$vars, tol=1e-4)
# should be the same if not conditional as well
sleepstudy2$w1 <- sleepstudy2$w2
suppressWarnings(wmw <- mix(highR ~ Days + (1|Subject), data=sleepstudy2, weights=c("w1", "w2"),
family=binomial(link="probit"),
cWeights=FALSE))
expect_equal(wmw$lnl, wm0$lnl, tol=1e-3)
expect_equal(coef(wmw), coef(wm0), tol=1e-4)
expect_equal(wmw$vars, wm0$vars, tol=1e-4)
# sub 1 conditional weights work
sleepstudy2 <- sleepstudy
sleepstudy2$highR <- ifelse(sleepstudy2$Reaction > 340, 1, 0)
sleepstudy2$w1 <- 1
sleepstudy2$w2 <- 1
sleepstudy2$w3 <- 1
sleepstudy2$Subject <- as.character(sleepstudy2$Subject)
g1 <- subset(sleepstudy2, sleepstudy2$Subject %in% c("310", "309", "349", "335"))
# do not renaem subject
sleepstudyrep <- rbind(sleepstudy2, g1)
sleepstudyrep$w1 <- ifelse(sleepstudyrep$Subject %in% c("310", "309", "349", "335"), 0.5, 1)
# two level
suppressWarnings(wm0 <- mix(highR ~ Days + (1|Subject), data=sleepstudyrep, weights=c("w1", "w2"),
family=binomial(link="probit"),verbose=FALSE, run=TRUE, cWeights=TRUE))
suppressWarnings(wmw <- mix(highR ~ Days + (1|Subject), data=sleepstudy2, weights=c("w1", "w2"),
family=binomial(link="probit"),verbose=FALSE, run=TRUE))
lm0 <- glmer(highR ~ Days + (1|Subject), data=sleepstudy2, family=binomial(link="probit"),
nAGQ=13)
# check vars
lmevars1 <- data.frame(summary(lm0)$varcor)$sdcor
expect_equal(unname(wm0$vars),
expected = unname(lmevars1)^2,
tolerance = 1e-3)
expect_equal(as.numeric(wm0$lnl), as.numeric(logLik(lm0)), tol=1e-3)
expect_equal(coef(wm0), fixef(lm0), tol=1e-4)
expect_equal(wmw$lnl, wm0$lnl, tol=1e-3)
expect_equal(coef(wmw), coef(wm0), tol=1e-4)
expect_equal(wmw$vars, wm0$vars, tol=1e-4)
# should be the same if not conditional as well
suppressWarnings(wm0 <- mix(highR ~ Days + (1|Subject), data=sleepstudyrep, weights=c("w1", "w2"),
family=binomial(link="probit"),verbose=FALSE, run=TRUE))
expect_equal(wmw$lnl, wm0$lnl, tol=1e-3)
expect_equal(coef(wmw), coef(wm0), tol=1e-4)
expect_equal(wmw$vars, wm0$vars, tol=1e-4)
})
context("Weighted v unweighted replicated two level model, poisson")
test_that("Weighted v unweighted replicated two level model, poisson", {
skip_on_cran()
skip_if_not_installed("glmmTMB")
require(glmmTMB)
owls2 <- Owls
owls2$w1 <- 1
owls2$w2 <- 1
owls2$w3 <- 1
g1 <- subset(owls2, owls2$Nest %in% c("Jeuss", "Lully", "Seiry", "Oleyes"))
g1$Nest <- paste0("R", g1$Nest)
owlsrep <- rbind(owls2, g1)
owls2$w2 <- ifelse(owls2$Nest %in% c("Jeuss", "Lully", "Seiry", "Oleyes"), 2, 1)
# glmer seems to have a bug with nAGQ > 1 for poisson models
# (see https://stat.ethz.ch/pipermail/r-sig-mixed-models/2020q2/028606.html),
# so can only compare with nQuad = 1
suppressWarnings(wm0 <- mix(SiblingNegotiation ~ ArrivalTime + (1|Nest), data=owlsrep, weights=c("w1", "w2"),
family=poisson(), verbose=FALSE,nQuad=1))
suppressWarnings(wmw <- mix(SiblingNegotiation ~ ArrivalTime + (1|Nest), data=owls2, weights=c("w1", "w2"), nQuad = 1,
family=poisson(), verbose=FALSE, cWeights=TRUE))
lm0 <- glmer(SiblingNegotiation ~ ArrivalTime + (1|Nest), data=owlsrep,
family="poisson")
# check vars
lmevars1 <- data.frame(summary(lm0)$varcor)$sdcor
expect_equal(unname(wm0$vars),
expected = unname(lmevars1)^2,
tolerance = 1e-3)
expect_equal(as.numeric(wm0$lnl), as.numeric(logLik(lm0)), tol=1e-3)
expect_equal(coef(wm0), fixef(lm0), tol=1e-4)
expect_equal(wmw$lnl, wm0$lnl, tol=1e-3)
expect_equal(coef(wmw), coef(wm0), tol=1e-4)
expect_equal(wmw$vars, wm0$vars, tol=1e-4)
# should be the same if not conditional as well
owls2$w1 <- owls2$w2
suppressWarnings(wmw <- mix(SiblingNegotiation ~ ArrivalTime + (1|Nest), data=owls2, weights=c("w1", "w2"),
nQuad = 1,family=poisson(), verbose=FALSE, cWeights=FALSE))
expect_equal(wmw$lnl, wm0$lnl, tol=1e-3)
expect_equal(coef(wmw), coef(wm0), tol=1e-4)
expect_equal(wmw$vars, wm0$vars, tol=1e-4)
# sub 1 conditional weights work
owls2 <- Owls
owls2$w1 <- 1
owls2$w2 <- 1
owls2$w3 <- 1
owls2$Nest <- as.character(owls2$Nest)
g1 <- subset(owls2, owls2$Nest %in% c("Jeuss", "Lully", "Seiry", "Oleyes"))
# do not renaem subject
owlsrep <- rbind(owls2, g1)
owlsrep$w1 <- ifelse(owlsrep$Nest %in% c("Jeuss", "Lully", "Seiry", "Oleyes"), 0.5, 1)
# two level
suppressWarnings(wm0 <- mix(SiblingNegotiation ~ ArrivalTime + (1|Nest), data=owlsrep, weights=c("w1", "w2"),
family=poisson(), verbose=FALSE,cWeights=TRUE,nQuad=1))
suppressWarnings(wmw <- mix(SiblingNegotiation ~ ArrivalTime + (1|Nest), data=owls2, weights=c("w1", "w2"), nQuad = 1,
family=poisson(), verbose=FALSE))
lm0 <- glmer(SiblingNegotiation ~ ArrivalTime + (1|Nest), data=owls2,
family="poisson")
# check vars
lmevars1 <- data.frame(summary(lm0)$varcor)$sdcor
expect_equal(unname(wm0$vars),
expected = unname(lmevars1)^2,
tolerance = 1e-3)
expect_equal(as.numeric(wm0$lnl), as.numeric(logLik(lm0)), tol=1e-3)
expect_equal(coef(wm0), fixef(lm0), tol=1e-4)
expect_equal(wmw$lnl, wm0$lnl, tol=1e-3)
expect_equal(coef(wmw), coef(wm0), tol=1e-4)
expect_equal(wmw$vars, wm0$vars, tol=1e-4)
# should be the same if not conditional as well
suppressWarnings(wm0 <- mix(SiblingNegotiation ~ ArrivalTime + (1|Nest), data=owlsrep, weights=c("w1", "w2"),
family=poisson(), verbose=FALSE,nQuad=1))
expect_equal(wmw$lnl, wm0$lnl, tol=1e-3)
expect_equal(coef(wmw), coef(wm0), tol=1e-4)
expect_equal(wmw$vars, wm0$vars, tol=1e-4)
})
context("Quadrature works, 3 RE at level 2")
test_that("Quadrature works, 3 RE at level 2", {
skip_on_cran()
skip_if_not_installed("glmmTMB")
require(glmmTMB)
owls2 <- Owls
owls2$w1 <- 1
owls2$w2 <- 1
owls2$w3 <- 1
suppressWarnings(wm0 <- mix(SiblingNegotiation ~ ArrivalTime + (ArrivalTime + NegPerChick|Nest),
data=owls2, weights=c("w1", "w2"),
family=poisson(),verbose=FALSE,nQuad=7))
expect_is(wm0, "WeMixResults")
})
context("Three level model slash and colon")
test_that("Three level model slash and colon", {
sleepstudy2 <- sleepstudy
sleepstudy2$Group <- 1
sleepstudy2$Group <- ifelse(sleepstudy2$Subject %in% c("310", "309", "349", "335"), 2, sleepstudy2$Group)
sleepstudy2$Group <- ifelse(sleepstudy2$Subject %in% c("330", "352", "337", "369"), 3, sleepstudy2$Group)
sleepstudy2$Group <- ifelse(sleepstudy2$Subject %in% c("331", "332", "334", "372"), 4, sleepstudy2$Group)
sleepstudy2$Group <- factor(sleepstudy2$Group)
# make subj that restarts at 1 per group, so there are four Subjects with subj=1, one per group.
sleepstudy2 <- sleepstudy2[order(sleepstudy2$Group, sleepstudy2$Subject),]
sleepstudy2$subj <- factor(rep(c(1:6,1:4,1:4,1:4),each=10))
# table(sleepstudy2$subj, sleepstudy2$Group) # shows each group has a subject 1:4 and group 1 as a 5 and 6 too
ss2 <- sleepstudy2
ss2$w1 <- w1c <- w1 <- rep(1,180)
w2c <- w2 <- rep(1,18)
w3c <- w3 <- rep(1,4)
ss2$w2 <- rep(1, 180)
ss2$w3 <- rep(1, 180)
lmr <- lmer(Reaction ~ Days + (1|Group) + (1|Subject), data=sleepstudy2, REML=FALSE)
# next line is a bad specification: confounds group=1, subj=1 person with group=2, subj=2 person
# lmr <- lmer(Reaction ~ Days + (1|Group) + (1|subj), data=sleepstudy2, REML=FALSE)
wm0 <- mix(Reaction ~ Days + (1|Group:subj) + (1|Group) , data=ss2, weights=c("w1", "w2","w3"))
expect_equal(wm0$lnl, as.numeric(logLik(lmr)), tol=1e-6)
expect_equal(coef(wm0), fixef(lmr), tol=1e-6)
expect_equal(unname(wm0$vars[length(wm0$vars)]), unname(lmr@devcomp$cmp["sigmaML"]^2), tol=1e-4)
#group mean center Days to test group mean centering as well
# we know the mean is 4.5 for all groups because all individuals were tested on all days
sleepstudy2$gmc_days <- sleepstudy2$Days-4.5
lmr <- lmer(Reaction ~ gmc_days + (1|Group/subj), data=sleepstudy2, REML=FALSE)
# above is same as:
# lmr <- lmer(Reaction ~ Days + (1|Group) + (1|Group:subj), data=sleepstudy2, REML=FALSE)
wm0 <- mix(Reaction ~ Days + (1|Group/subj), data=ss2, weights=c("w1", "w2","w3"),center_group = list("Group/subj" = ~Days))
expect_equal(wm0$lnl, as.numeric(logLik(lmr)), tol=1e-6)
expect_equal(unname(coef(wm0)), unname(fixef(lmr)), tol=1e-6)
expect_equal(unname(wm0$vars[length(wm0$vars)]), unname(lmr@devcomp$cmp["sigmaML"]^2), tol=1e-4)
})
# check the format of summary output
context("Summary output format")
test_that("summary output format", {
skip_on_cran()
sleepstudy2 <- sleepstudy
sleepstudy2$Group <- 1
sleepstudy2$Group <- ifelse(sleepstudy2$Subject %in% c("310", "309", "349", "335"), 2, sleepstudy2$Group)
sleepstudy2$Group <- ifelse(sleepstudy2$Subject %in% c("330", "352", "337", "369"), 3, sleepstudy2$Group)
sleepstudy2$Group <- ifelse(sleepstudy2$Subject %in% c("331", "332", "334", "372"), 4, sleepstudy2$Group)
sleepstudy2$Group <- factor(sleepstudy2$Group)
# make subj that restarts at 1 per group, so there are four Subjects with subj=1, one per group.
sleepstudy2 <- sleepstudy2[order(sleepstudy2$Group, sleepstudy2$Subject),]
sleepstudy2$subj <- factor(rep(c(1:6,1:4,1:4,1:4),each=10))
sleepstudy2$carPr <- pnorm(sleepstudy2$Reaction -285,sd=50)
set.seed(2)
sleepstudy2$car <- runif(180) < sleepstudy2$carPr
# table(sleepstudy2$subj, sleepstudy2$Group) # shows each group has a subject 1:4 and group 1 as a 5 and 6 too
ss2 <- sleepstudy2
ss2$w1 <- w1c <- w1 <- rep(1,180)
w2c <- w2 <- rep(1,18)
w3c <- w3 <- rep(1,4)
ss2$w2 <- rep(1, 180)
ss2$w3 <- rep(1, 180)
co0 <- c("Call:",
"mix(formula = Reaction ~ Days + (Days + car | Subject), data = ss2, ",
" weights = c(\"w1\", \"w2\"))", "",
"Variance terms:",
" Level Group Name Variance Std. Error Std.Dev. Corr1 Corr2",
" 2 Subject (Intercept) 629.4 254.1 25.09 ",
" 2 Subject Days 36.7 12.2 6.06 0.17 ",
" 2 Subject carTRUE 648.5 309.3 25.46 -0.42 -0.39",
" 1 Residual 528.8 145.7 23.00 ",
"Groups:",
" Level Group n size mean wgt sum wgt",
" 2 Subject 18 1 18",
" 1 Obs 180 1 180",
"",
"Fixed Effects:",
" Estimate Std. Error t value",
"(Intercept) 252.72 6.53 38.70",
"Days 11.19 1.50 7.44",
"",
"lnl= -868.13 ",
"Intraclass Correlation= 0.713 ")
wm0 <- mix(Reaction ~ Days + (Days+car|Subject), data=ss2, weights=c("w1", "w2"))
withr::with_options(list(digits=2),
co <- capture.output(summary(wm0))
)
expect_equal(co, co0)
# this model is not accuractly estimated so is not a good test of output format.
# nevertheless, it uses data from this block, so appears here
expect_warning(wm0 <- mix(Reaction ~ Days + (car||Subject), data=ss2, weights=c("w1", "w2")), "standard error of variance terms")
varsref <- structure(list(Variance = c(698.125, 213.50 , 1618.92, 808.25 ),
`Std. Error` = c(385.8 , 232 , 850.5 , 162.0972),
Std.Dev. = c(26.4220, 14.6116, 40.2358, 28.42977)),
class = "data.frame",
row.names = c("Subject.(Intercept)", "Subject.carFALSE", "Subject.carTRUE", "Residual"))
coefref <- structure(c(236.842700778804, 9.59410909072588, 6.7061583949811,
1.56079000196884, 35.317194558968, 6.14695704010372), .Dim = 2:3,
.Dimnames = list(c("(Intercept)", "Days"),
c("Estimate", "Std. Error", "t value")))
sum0 <- summary(wm0)
expect_equal(sum0$coef, coefref, tolerance=1e-4)
expect_equal(sum0$vars, varsref, tolerance=4e-3)
})
context("Weighted three level model")
test_that("Weighted three level model", {
sleepstudy2 <- sleepstudy
sleepstudy2$Group <- 1
sleepstudy2$Group <- ifelse(sleepstudy2$Subject %in% c("310", "309", "349", "335"), 2, sleepstudy2$Group)
sleepstudy2$Group <- ifelse(sleepstudy2$Subject %in% c("330", "352", "337", "369"), 3, sleepstudy2$Group)
sleepstudy2$Group <- ifelse(sleepstudy2$Subject %in% c("331", "332", "334", "372"), 4, sleepstudy2$Group)
sleepstudy2$Group <- factor(sleepstudy2$Group)
#sleepstudy2$NewVar <- sleepstudy2$Reaction/100 + rnorm(180,10,5)
sleepstudy2$NewVar <- c(9.43457771946133, 12.7183652886321, 11.242490364576, 8.32049769489214, 21.4769136458094, 12.801691804759, 11.2944480319612, 23.2877967612343, 17.6554807720438, 15.6839576150421,
12.0415732799454, 12.6326929342708, 10.7497191412678, 8.67408880103278, 16.7872963953428, 14.6248001125594, 14.3553882674867, 11.8275049579491, 18.7525279773655, 5.57650843995549,
19.7623973348143, 19.1858339464439, 16.9845555399241, 5.79214586016569, 13.6760775160985, 8.20870365747105, 10.469622698121, 9.5294364080553, 0.474272552642729, 9.90998530696626,
7.85598926755497, 13.6908324775899, 8.363578606421, 12.4585008913942, 20.8371544197279, 6.95572306644593, 4.49104017816579, 13.529348727477, 8.31255622061886, 8.78452997885101,
17.216378291436, 13.6427335980228, 19.498592420836, 7.77906722601152, 15.1266474871889, 12.8030124499139, 6.62054433117513, 8.66974559938289, 11.9540398510289, 14.4271543736958,
10.0496588773239, 15.5177485124749, 9.12750122498998, 10.1801282848309, 14.2562515194416, 19.3249322019648, 10.7021090784807, 19.5336266451282, 7.58133671725254, 19.889529248361,
5.32630146259616, 13.6044477039579, 4.98570236262749, 23.8801945918399, 14.447634009947, 15.6043025989833, 22.411260279495, 16.66695588071, 13.2723030644539, 17.6752877679585,
13.4050771535649, 16.8476633007641, 14.4754165647234, 13.2597805356783, 12.2720489082021, 16.3515751846435, 18.2589186962967, 13.5608781488381, 22.128982372665, 14.0340540905538,
2.08454544532827, 10.8863647522786, 21.0103204678449, 18.4595691535641, 11.9914499297614, 9.96791953740012, 20.0728999616675, 10.7237225063556, 18.6117287544885, 15.6290844263721,
13.584186380461, 12.8188973592142, 8.73588740878364, 22.1049445716224, 8.82190889109436, 16.7501795316447, 11.7894955061621, 8.67751729426691, 8.28560204932531, 25.3725616103908,
14.0620375346182, 14.9882805901214, 17.181272042065, 19.6581190406994, 11.6757977643058, 16.638959090019, 9.68436466184703, 23.1282879805077, 20.7690736926973, 14.3401954001056,
16.2747349705209, 8.61220295948687, 19.869302293393, 5.78478905807136, 6.59988762289215, 20.489844013367, 27.239468499143, 20.5797203990527, 19.4544929129352, 12.5320006698243,
18.2664525619468, 16.6022500953422, 13.0216995414837, 12.4573456587016, 18.6922622366627, 7.37338081469625, 18.2016743123929, 12.9216600044986, 17.2543801572655, 20.1850819959469,
9.46342758905614, 16.1505649461445, 11.644647125078, 17.6992809329738, 13.5726089392145, 11.9006738563255, 9.11529530502925, 11.9473408244809, 13.7037965361481, 21.052192466607,
8.94718178506257, 15.165337112518, 8.23261412105177, 12.2188157585207, 14.1455895553682, 4.94479286361234, 16.9113105275149, 12.7165809453375, 5.58774136619933, 8.31670629685991,
6.67728085523989, 11.7741431392146, 7.91590113917014, 9.13428547474213, 18.5689544076891, 18.1995195536489, 13.3510327364336, 21.4047919488063, 17.7821335119331, 9.19515401988151,
10.7918240635605, 6.63952336871624, 15.3733302621749, 15.5262403960572, 13.4421915705244, 11.5980148647886, 17.3920393426159, 6.6883267354037, 10.6865661488425, 15.1104585129107,
9.99190979270527, 16.54240123493, 13.5687018760124, 13.5592897858269, 19.4019579232211, 15.0561754718535, 9.36036662501835, 4.27173388144395, 19.2040511963462, 6.36268714321505)
ss2 <- sleepstudy2
w1c <- w1 <- rep(1,180)
w2c <- w2 <- rep(1,18)
w3c <- w3 <- rep(1,4)
# unbalanced (non-identical within a Subject), level-1 (obs level) weights
w1c[1:4] <- w1[1:4] <- 2
uR <- sleepstudy2[1:4,]
sleepstudy2 <- rbind(sleepstudy2, uR)
# level-2 weights
w2c[1] <- w2[1] <- 2
w1[ss2$Subject=="308"] <- 2*w1[ss2$Subject=="308"]
sR <- subset(sleepstudy2, Subject == "308")
sR$Subject <- "S308"
sleepstudy2 <- rbind(sleepstudy2, sR)
# level-3 weights
w3c[1] <- w3[1] <- 2
w2[c(1,7,12,13,16,17)] <- 2*w2[c(1,7,12,13,16,17)]
w1[ss2$Group==1] <- 2*w1[ss2$Group==1]
gR <- subset(sleepstudy2, Group %in% 1)
gR$Subject <- paste0("G", gR$Subject)
gR$Group <- "11"
sleepstudy2 <- rbind(sleepstudy2, gR)
ss2$w1 <- w1
ss2$w2 <- rep(w2,each=10)
ss2$w3 <- ifelse(ss2$Subject %in% c("308", "333", "350", "351", "370", "371"),2,1)
# lmr for reference
lmr <- lmer(Reaction ~ Days + (1|Subject) + (1|Group), data=sleepstudy2, REML=FALSE)
#these are the conditional weights, used in the stata comparison
ss2$c1 <- w1c
ss2$c2 <- rep(w2c,each=10)
ss2$c3 <- ifelse(ss2$Subject %in% c("308", "333", "350", "351", "370", "371"),2,1)
#compare against mixed stata resutls
wm0 <- mix(Reaction ~ Days + (1|Subject) + (1|Group), data=ss2, weights=c("w1", "w2","w3"))
expect_equal(wm0$lnl, -1389.0983, tolerance=1e-3)
expect_equal(unname(wm0$coef), c(243.6831,12.67954),tolerance = 1e-3)
expect_equal(unname(wm0$vars), c(360.97 ,756.5857, 1153.704 ),tolerance = 1e-3)
expect_equal(unname( wm0$SE), c( 11.68881 , 2.468474 ),tolerance = 1e-3)
wm1 <- mix(Reaction ~ Days + (1 |Subject) + (1|Group)+ (0+Days|Group), data=ss2, weights=c("w1", "w2","w3"))
expect_equal(wm1$lnl, -1377.6876 , tolerance=1e-3)
expect_equal(unname(wm1$coef), c(247.6234 , 11.71211 ),tolerance = 1e-3)
expect_equal(unname(wm1$vars), c(397.6301 , 220.4182 , 16.86835 , 1022.88 ),tolerance = 1e-3)
expect_equal(unname( wm1$SE), c(9.97378 , 2.60348 ),tolerance = 1e-3)
})
context("Wald Test")
test_that("Test for Wald Tests using ", {
wm1<- mix(Reaction ~ Days + (0+Days|Subject), data=sleepstudyU, weights=c("weight1L1", "weight1L2"))
beta_1 <- waldTest(wm1, type="beta", coefs="Days", hypothesis=10.4073- 1.5458 *1.96)
#wald test and T test are the same for individual betas
expect_equal(beta_1$p,.05,.01)
beta_2<- waldTest(wm1, type="beta", coefs="(Intercept)", hypothesis=251.4051- 6.8246 *1.96)
#wald test and T test are the same for individual betas
expect_equal(beta_2$p,.05,.01)
# Ensure that tests for Lambda run
wm2<- mix(Reaction ~ Days + (1+Days|Subject), data=sleepstudyU, weights=c("weight1L1", "weight1L2"))
waldTest(wm2,type="Lambda")
waldTest(wm2,type="Lambda",coefs = c("Subject.Days.(Intercept)","Subject.Days" ))
})
context("Complex weighted three level model")
test_that("Complex weighted three level model", {
skip_on_cran()
sleepstudy2 <- sleepstudy
sleepstudy2$Group <- 1
sleepstudy2$Group <- ifelse(sleepstudy2$Subject %in% c("310", "309", "349", "335"), 2, sleepstudy2$Group)
sleepstudy2$Group <- ifelse(sleepstudy2$Subject %in% c("330", "352", "337", "369"), 3, sleepstudy2$Group)
sleepstudy2$Group <- ifelse(sleepstudy2$Subject %in% c("331", "332", "334", "372"), 4, sleepstudy2$Group)
sleepstudy2$Group <- factor(sleepstudy2$Group)
sleepstudy2$NewVar <- c(9.43457771946133, 12.7183652886321, 11.242490364576, 8.32049769489214, 21.4769136458094, 12.801691804759, 11.2944480319612, 23.2877967612343, 17.6554807720438, 15.6839576150421,
12.0415732799454, 12.6326929342708, 10.7497191412678, 8.67408880103278, 16.7872963953428, 14.6248001125594, 14.3553882674867, 11.8275049579491, 18.7525279773655, 5.57650843995549,
19.7623973348143, 19.1858339464439, 16.9845555399241, 5.79214586016569, 13.6760775160985, 8.20870365747105, 10.469622698121, 9.5294364080553, 0.474272552642729, 9.90998530696626,
7.85598926755497, 13.6908324775899, 8.363578606421, 12.4585008913942, 20.8371544197279, 6.95572306644593, 4.49104017816579, 13.529348727477, 8.31255622061886, 8.78452997885101,
17.216378291436, 13.6427335980228, 19.498592420836, 7.77906722601152, 15.1266474871889, 12.8030124499139, 6.62054433117513, 8.66974559938289, 11.9540398510289, 14.4271543736958,
10.0496588773239, 15.5177485124749, 9.12750122498998, 10.1801282848309, 14.2562515194416, 19.3249322019648, 10.7021090784807, 19.5336266451282, 7.58133671725254, 19.889529248361,
5.32630146259616, 13.6044477039579, 4.98570236262749, 23.8801945918399, 14.447634009947, 15.6043025989833, 22.411260279495, 16.66695588071, 13.2723030644539, 17.6752877679585,
13.4050771535649, 16.8476633007641, 14.4754165647234, 13.2597805356783, 12.2720489082021, 16.3515751846435, 18.2589186962967, 13.5608781488381, 22.128982372665, 14.0340540905538,
2.08454544532827, 10.8863647522786, 21.0103204678449, 18.4595691535641, 11.9914499297614, 9.96791953740012, 20.0728999616675, 10.7237225063556, 18.6117287544885, 15.6290844263721,
13.584186380461, 12.8188973592142, 8.73588740878364, 22.1049445716224, 8.82190889109436, 16.7501795316447, 11.7894955061621, 8.67751729426691, 8.28560204932531, 25.3725616103908,
14.0620375346182, 14.9882805901214, 17.181272042065, 19.6581190406994, 11.6757977643058, 16.638959090019, 9.68436466184703, 23.1282879805077, 20.7690736926973, 14.3401954001056,
16.2747349705209, 8.61220295948687, 19.869302293393, 5.78478905807136, 6.59988762289215, 20.489844013367, 27.239468499143, 20.5797203990527, 19.4544929129352, 12.5320006698243,
18.2664525619468, 16.6022500953422, 13.0216995414837, 12.4573456587016, 18.6922622366627, 7.37338081469625, 18.2016743123929, 12.9216600044986, 17.2543801572655, 20.1850819959469,
9.46342758905614, 16.1505649461445, 11.644647125078, 17.6992809329738, 13.5726089392145, 11.9006738563255, 9.11529530502925, 11.9473408244809, 13.7037965361481, 21.052192466607,
8.94718178506257, 15.165337112518, 8.23261412105177, 12.2188157585207, 14.1455895553682, 4.94479286361234, 16.9113105275149, 12.7165809453375, 5.58774136619933, 8.31670629685991,
6.67728085523989, 11.7741431392146, 7.91590113917014, 9.13428547474213, 18.5689544076891, 18.1995195536489, 13.3510327364336, 21.4047919488063, 17.7821335119331, 9.19515401988151,
10.7918240635605, 6.63952336871624, 15.3733302621749, 15.5262403960572, 13.4421915705244, 11.5980148647886, 17.3920393426159, 6.6883267354037, 10.6865661488425, 15.1104585129107,
9.99190979270527, 16.54240123493, 13.5687018760124, 13.5592897858269, 19.4019579232211, 15.0561754718535, 9.36036662501835, 4.27173388144395, 19.2040511963462, 6.36268714321505)
ss2 <- sleepstudy2
w1c <- w1 <- rep(1,180)
w2c <- w2 <- rep(1,18)
w3c <- w3 <- rep(1,4)
# unbalanced (non-identical within a Subject), level-1 (obs level) weights
w1c[1:4] <- w1[1:4] <- 2
uR <- sleepstudy2[1:4,]
sleepstudy2 <- rbind(sleepstudy2, uR)
# level-2 weights
w2c[1] <- w2[1] <- 2
w1[ss2$Subject=="308"] <- 2*w1[ss2$Subject=="308"]
sR <- subset(sleepstudy2, Subject == "308")
sR$Subject <- "S308"
sleepstudy2 <- rbind(sleepstudy2, sR)
# level-3 weights
w3c[1] <- w3[1] <- 2
w2[c(1,7,12,13,16,17)] <- 2*w2[c(1,7,12,13,16,17)]
w1[ss2$Group==1] <- 2*w1[ss2$Group==1]
gR <- subset(sleepstudy2, Group %in% 1)
gR$Subject <- paste0("G", gR$Subject)
gR$Group <- "11"
sleepstudy2 <- rbind(sleepstudy2, gR)
ss2$w1 <- w1
ss2$w2 <- rep(w2,each=10)
ss2$w3 <- ifelse(ss2$Subject %in% c("308", "333", "350", "351", "370", "371"),2,1)
ss2$n1 <- ss2$n2 <- ss2$n3 <- 1
#lme reference with duplicated data
lmr <- lmer(Reaction ~ Days + (1|Subject) + (1+Days|Group), data=sleepstudy2, REML=FALSE)
#Does WeMix match to duplicated lmr with one+two random effects
wmr <- mix(Reaction ~ Days + (1|Subject) + (1+Days|Group), data=ss2,weights=c("w1","w2","w3"))
expect_equal(wmr$lnl, as.numeric(logLik(lmr)), tol=1e-3)
expect_equal(coef(wmr), fixef(lmr), tol=1e-4)
lmevars1 <- data.frame(summary(lmr)$varcor)
vars <- lmevars1[is.na(lmevars1$var2),"vcov"]
expect_equal(unname(wmr$vars),
expected = unname(vars),
tolerance = 1e-3)
expect_equal(coef(wmr), fixef(lmr), tol=1e-4)
expect_equal(unname(wmr$vars[length(wmr$vars)]), unname(lmr@devcomp$cmp["sigmaML"]^2), tol=1e-4)
#lme reference with duplicated data
lmr2 <- lmer(Reaction ~ Days + (1 + NewVar |Subject) + (1+Days|Group),
data=sleepstudy2, control=lmerControl(optimizer="bobyqa"), REML=FALSE)
#wemix with weights
suppressWarnings(wmr2 <- mix(Reaction ~ Days + (1+ NewVar |Subject) + (1+Days|Group), data=ss2,weights=c("w1","w2","w3")))
#Does WeMix match to duplicated lmr with two correlated random effects at two levesl
expect_equal(wmr2$lnl, as.numeric(logLik(lmr2)), tol=1e-3)
expect_equal(coef(wmr2), fixef(lmr2), tol=1e-4)
lmevars2 <- data.frame(summary(lmr2)$varcor)
vars <- lmevars2[is.na(lmevars2$var2),"vcov"]
expect_equal(unname(wmr2$vars),
expected = unname(vars),
tolerance = 1e-3)
expect_equal(coef(wmr2), fixef(lmr2), tol=1e-4)
expect_equal(unname(wmr2$vars[length(wmr2$vars)]), unname(lmr2@devcomp$cmp["sigmaML"]^2), tol=1e-4)
})
context("Complex weighted three level model, binomial")
test_that("Complex weighted three level model, binomial", {
skip_on_cran()
sleepstudy2 <- sleepstudy
sleepstudy2$highR <- ifelse(sleepstudy2$Reaction > 340, 1, 0)
sleepstudy2$Group <- 1
sleepstudy2$Group <- ifelse(sleepstudy2$Subject %in% c("310", "309", "349", "335"), 2, sleepstudy2$Group)
sleepstudy2$Group <- ifelse(sleepstudy2$Subject %in% c("330", "352", "337", "369"), 3, sleepstudy2$Group)
sleepstudy2$Group <- ifelse(sleepstudy2$Subject %in% c("331", "332", "334", "372"), 4, sleepstudy2$Group)
sleepstudy2$Group <- factor(sleepstudy2$Group)
sleepstudy2$NewVar <- c(9.43457771946133, 12.7183652886321, 11.242490364576, 8.32049769489214, 21.4769136458094, 12.801691804759, 11.2944480319612, 23.2877967612343, 17.6554807720438, 15.6839576150421,
12.0415732799454, 12.6326929342708, 10.7497191412678, 8.67408880103278, 16.7872963953428, 14.6248001125594, 14.3553882674867, 11.8275049579491, 18.7525279773655, 5.57650843995549,
19.7623973348143, 19.1858339464439, 16.9845555399241, 5.79214586016569, 13.6760775160985, 8.20870365747105, 10.469622698121, 9.5294364080553, 0.474272552642729, 9.90998530696626,
7.85598926755497, 13.6908324775899, 8.363578606421, 12.4585008913942, 20.8371544197279, 6.95572306644593, 4.49104017816579, 13.529348727477, 8.31255622061886, 8.78452997885101,
17.216378291436, 13.6427335980228, 19.498592420836, 7.77906722601152, 15.1266474871889, 12.8030124499139, 6.62054433117513, 8.66974559938289, 11.9540398510289, 14.4271543736958,
10.0496588773239, 15.5177485124749, 9.12750122498998, 10.1801282848309, 14.2562515194416, 19.3249322019648, 10.7021090784807, 19.5336266451282, 7.58133671725254, 19.889529248361,
5.32630146259616, 13.6044477039579, 4.98570236262749, 23.8801945918399, 14.447634009947, 15.6043025989833, 22.411260279495, 16.66695588071, 13.2723030644539, 17.6752877679585,
13.4050771535649, 16.8476633007641, 14.4754165647234, 13.2597805356783, 12.2720489082021, 16.3515751846435, 18.2589186962967, 13.5608781488381, 22.128982372665, 14.0340540905538,
2.08454544532827, 10.8863647522786, 21.0103204678449, 18.4595691535641, 11.9914499297614, 9.96791953740012, 20.0728999616675, 10.7237225063556, 18.6117287544885, 15.6290844263721,
13.584186380461, 12.8188973592142, 8.73588740878364, 22.1049445716224, 8.82190889109436, 16.7501795316447, 11.7894955061621, 8.67751729426691, 8.28560204932531, 25.3725616103908,
14.0620375346182, 14.9882805901214, 17.181272042065, 19.6581190406994, 11.6757977643058, 16.638959090019, 9.68436466184703, 23.1282879805077, 20.7690736926973, 14.3401954001056,
16.2747349705209, 8.61220295948687, 19.869302293393, 5.78478905807136, 6.59988762289215, 20.489844013367, 27.239468499143, 20.5797203990527, 19.4544929129352, 12.5320006698243,
18.2664525619468, 16.6022500953422, 13.0216995414837, 12.4573456587016, 18.6922622366627, 7.37338081469625, 18.2016743123929, 12.9216600044986, 17.2543801572655, 20.1850819959469,
9.46342758905614, 16.1505649461445, 11.644647125078, 17.6992809329738, 13.5726089392145, 11.9006738563255, 9.11529530502925, 11.9473408244809, 13.7037965361481, 21.052192466607,
8.94718178506257, 15.165337112518, 8.23261412105177, 12.2188157585207, 14.1455895553682, 4.94479286361234, 16.9113105275149, 12.7165809453375, 5.58774136619933, 8.31670629685991,
6.67728085523989, 11.7741431392146, 7.91590113917014, 9.13428547474213, 18.5689544076891, 18.1995195536489, 13.3510327364336, 21.4047919488063, 17.7821335119331, 9.19515401988151,
10.7918240635605, 6.63952336871624, 15.3733302621749, 15.5262403960572, 13.4421915705244, 11.5980148647886, 17.3920393426159, 6.6883267354037, 10.6865661488425, 15.1104585129107,
9.99190979270527, 16.54240123493, 13.5687018760124, 13.5592897858269, 19.4019579232211, 15.0561754718535, 9.36036662501835, 4.27173388144395, 19.2040511963462, 6.36268714321505)
ss2 <- sleepstudy2
w1c <- w1 <- rep(1,180)
w2c <- w2 <- rep(1,18)
w3c <- w3 <- rep(1,4)
# unbalanced (non-identical within a Subject), level-1 (obs level) weights
w1c[1:4] <- w1[1:4] <- 2
uR <- sleepstudy2[1:4,]
sleepstudy2 <- rbind(sleepstudy2, uR)
# level-2 weights
w2c[1] <- w2[1] <- 2
w1[ss2$Subject=="308"] <- 2*w1[ss2$Subject=="308"]
sR <- subset(sleepstudy2, Subject == "308")
sR$Subject <- "S308"
sleepstudy2 <- rbind(sleepstudy2, sR)
# level-3 weights
w3c[1] <- w3[1] <- 2
w2[c(1,7,12,13,16,17)] <- 2*w2[c(1,7,12,13,16,17)]
w1[ss2$Group==1] <- 2*w1[ss2$Group==1]
gR <- subset(sleepstudy2, Group %in% 1)
gR$Subject <- paste0("G", gR$Subject)
gR$Group <- "11"
sleepstudy2 <- rbind(sleepstudy2, gR)
ss2$w1 <- w1
ss2$w2 <- rep(w2,each=10)
ss2$w3 <- ifelse(ss2$Subject %in% c("308", "333", "350", "351", "370", "371"),2,1)
ss2$n1 <- ss2$n2 <- ss2$n3 <- 1
#lme reference with duplicated data
glmr <- glmer(highR ~ Days + (1|Subject) + (1+Days|Group), data=sleepstudy2,
family=binomial())
#Does WeMix match to duplicated lmr with one+two random effects
suppressWarnings(wmr <- mix(highR ~ Days + (1|Subject) + (Days|Group), data=ss2,
weights=c("w1","w2","w3"),
family=binomial(),nQuad=1,verbose=FALSE))
expect_equal(as.numeric(wmr$lnl), as.numeric(logLik(glmr)), tol=1e-3)
expect_equal(coef(wmr), fixef(glmr), tol=1e-4)
lmevars1 <- data.frame(summary(glmr)$varcor)
vars <- lmevars1$vcov
expect_equal(unname(wmr$vars),
expected = unname(vars),
tolerance = 1e-3)
#lme reference with duplicated data
glmr2 <- glmer(highR ~ Days + (1 + NewVar|Subject) + (1+Days|Group),
data=sleepstudy2,
family=binomial(),control=glmerControl(optimizer="bobyqa"))
#wemix with weights
suppressWarnings(wmr2 <- mix(highR ~ Days + (1+ NewVar |Subject) + (1+Days|Group), data=ss2,
weights=c("w1","w2","w3"), family=binomial(),nQuad=1,
max_iteration = 20,verbose=FALSE))
#Does WeMix match to duplicated lmr with two correlated random effects at two levesl
expect_equal(as.numeric(wmr2$lnl), as.numeric(logLik(glmr2)), tol=1e-3)
expect_equal(coef(wmr2), fixef(glmr2), tol=1e-4)
lmevars2 <- data.frame(summary(glmr2)$varcor)
vars <- lmevars2$vcov
expect_equal(unname(wmr2$vars),
expected = unname(vars),
tolerance = 1e-3)
})
context("Model Matrix has a hard time with")
test_that("Model Matrix has a hard time with", {
skip_on_cran()
require(EdSurvey)
sdf <- readNAEP(system.file("extdata/data", "M36NT2PM.dat", package = "NAEPprimer"))
gg <- EdSurvey::getData(varnames=c("composite", "dsex", "b017451", "scrpsu", "origwt", "smsrswt"), data=sdf, returnJKreplicates=FALSE)
gg2 <- gg[gg$origwt > 0 & gg$smsrswt > 0,]
suppressMessages(m4 <- mix(mrpcm2 ~ dsex + b017451 + (1|scrpsu), data=gg2, weights=c("origwt", "smsrswt")))
expect_equal(m4$lnl, -81882.3634148408, tol=1e-5)
expect_equal(m4$coef, c(`(Intercept)` = 270.571817625202, dsexFemale = -2.14907551600309, `b017451Once every few weeks` = 3.84433128452533, `b017451About once a week` = 9.19954009166631, `b0174512 or 3 times a week` = 12.8701977366809, `b017451Every day` = 6.29635843233831))
expect_equal(m4$SE, c(`(Intercept)` = 1.14859092700564, dsexFemale = 0.636207386283148, `b017451Once every few weeks` = 1.05190060207538, `b017451About once a week` = 0.986062663493941, `b0174512 or 3 times a week` = 1.00366545196335, `b017451Every day` = 1.10858750669702))
m4varDF <- structure(list(grp = c("scrpsu", "Residual"),
var1 = c("(Intercept)", NA),
var2 = c(NA_character_, NA_character_),
vcov = c(300.695102536901, 969.501918646696),
ngrp = c(672, 16331),
level = c(2, 1),
SEvcov = c(25.0004604746421, 16.2937862202021),
fullGroup = c("scrpsu.(Intercept)", "Residual")),
row.names = c(NA, -2L),
class = "data.frame")
expect_equal(m4$varDF, m4varDF, tol=1e-5)
})
context("examples run")
test_that("examples run", {
skip_on_cran()
ss1 <- sleepstudy
# Create weights
ss1$W1 <- ifelse(ss1$Subject %in% c(308, 309, 310), 2, 1)
ss1$W2 <- 1
# Run random-intercept 2-level model
two_level <- mix(Reaction ~ Days + (1|Subject), data=ss1, weights=c("W1", "W2"))
expect_is(two_level, "WeMixResults")
#Run random-intercept 2-level model with group-mean centering
grp_centered <- mix(Reaction ~ Days + (1|Subject), data=ss1,
weights = c("W1", "W2"),
center_group = list("Subject" = ~Days))
expect_is(grp_centered, "WeMixResults")
#Run three level model with random slope and intercept.
#add group variables for 3 level model
ss1$Group <- 3
ss1$Group <- ifelse(as.numeric(ss1$Subject) %% 10 < 7, 2, ss1$Group)
ss1$Group <- ifelse(as.numeric(ss1$Subject) %% 10 < 4, 1, ss1$Group)
# level-3 weights
ss1$W3 <- ifelse(ss1$Group == 2, 2, 1)
three_level <- mix(Reaction ~ Days + (1|Subject) + (1+Days|Group), data=ss1,
weights=c("W1", "W2", "W3"))
expect_is(three_level, "WeMixResults")
})
Try the WeMix package in your browser
Any scripts or data that you put into this service are public.
WeMix documentation built on Nov. 3, 2023, 9:06 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
require(testthat)
# for sleepstudy, calls in here to lmer
library(lme4)
# for grad function
options(width = 500)
options(useFancyQuotes = FALSE)
### Data Used: sleepstudy
tolerance <- 1E-3
data("sleepstudy")
### Unweigted =====================================
sleepstudyU <- sleepstudy
sleepstudyU$weight1L1 <- 1
sleepstudyU$weight1L2 <- 1
context("The model runs")
test_that("The model runs", {
wm0 <- mix(Reaction ~ Days + (1|Subject), data=sleepstudyU, weights=c("weight1L1", "weight1L2"))
lme1 <- lmer(Reaction ~ Days + (1|Subject), data=sleepstudyU, REML=FALSE)
expect_equal(wm0$lnl, -897.039321502613, tolerance=tolerance*897) # value from lmer(Reaction ~ Days + (1 | Subject), data=sleepstudy, REML=FALSE)
expect_equal(unname(c(wm0$coef)), unname(fixef(lme1)), tolerance=tolerance)
expect_equal(unname(wm0$theta), unname(lme1@theta), tolerance=tolerance)
expect_equal(coef(wm0),
expected = getME(lme1, "fixef"),
tolerance = tolerance)
lmevars1 <- data.frame(summary(lme1)$varcor)$sdcor
expect_equal(unname(wm0$vars),
expected = unname(lmevars1)^2,
tolerance = tolerance)
# agrees with GLLAMM
#source: logit_command.do
gllamm_model1 <- c("(Intercept)" = 251.4051,
"Days" = 10.46729,
"Residual" = 1296.8768,
"Subject" = 954.52789,
"lnl" = -897.03932)
expect_equal(unname(coef(wm0)),
expected = unname(gllamm_model1[1:2]),
tolerance = abs(tolerance))
expect_equal(unname(wm0$vars),
expected = unname(gllamm_model1[3:4]),
tolerance = tolerance)
expect_equal(unname(wm0$lnl),
expected=unname(gllamm_model1[5]),
tolerance=abs(tolerance))
sleepstudy2 <- sleepstudy
sleepstudy2$block <- rep(1:6,each=30)
sleepstudy2$weight1L1 <- sleepstudy2$weight1L2 <- sleepstudy2$weight1L3 <- 1
# change Reaction by the artifical block we just added
sleepstudy2$Reaction <- sleepstudy2$Reaction + sleepstudy2$block*20
names(sleepstudy2) <- gsub("weight1L","pwt",names(sleepstudy2))
lmr <- lmer(Reaction ~ Days + (1|Subject) + (1|block), data=sleepstudy2, REML=FALSE)
mm <- mix(Reaction ~ Days + (1|Subject) + (1|block), data=sleepstudy2, weights=c("pwt1", "pwt2", "pwt3"))
expect_equal(unname(mm$lnl),
expected=unname(logLik(lmr)[[1]]),
tolerance=tolerance)
})
context("Factor binomial")
test_that("Factor binomial", {
sleepstudyM <- sleepstudyU
sleepstudyM$highR <- factor(ifelse(sleepstudyM$Reaction > 340, 2, 1), levels=c(1,2), labels=c("L","H"))
sleepstudyM$Sub <- factor(sleepstudyM$Subject, levels=300:400)
m1 <- mix(highR ~ Days + (1|Subject), data=sleepstudyM, weights=c("weight1L1", "weight1L2"),
family="binomial")
expect_equal(coef(m1), c(`(Intercept)` = -8.69336570860094, Days = 1.17236933934874), tol=1e-3)
summaryREF <- c("Call:",
"mix(formula = highR ~ Days + (1 | Subject), data = sleepstudyM, ",
" weights = c(\"weight1L1\", \"weight1L2\"), family = \"binomial\")",
"",
"Variance terms:",
" Level Group Name Variance Std. Error Std.Dev.",
" 2 Subject (Intercept) 7.99 5.14 2.83", "Groups:",
" Level Group n size mean wgt sum wgt",
" 2 Subject 18 1 18",
" 1 Obs 180 1 180",
"",
"Fixed Effects:",
" Estimate Std. Error t value",
"(Intercept) -8.694 2.073 -4.19",
"Days 1.172 0.278 4.21",
"",
"lnl= -53.90 ")
withr::with_options(list(digits=2),
co <- capture.output(summary(m1))
)
expect_equal(co, summaryREF)
})
context("Unweighted model with 2 random effects")
test_that("Agrees with lme4 3,handles missing data", {
sleepstudyM <- sleepstudyU
sleepstudyM <- sleepstudyM[order(runif(nrow(sleepstudyM))), ]
#introduce a missing value
sleepstudyM$Days[3] <- NA
sleepstudyM$Subject <- ifelse(sleepstudyM$Subject %in% c("310"),"2",as.character(sleepstudyM$Subject))
sleepstudyM$Subject <- ifelse(sleepstudyM$Subject %in% c("309"),"1",as.character(sleepstudyM$Subject))
sleepstudyM$Subject <- ifelse(sleepstudyM$Subject %in% c("330"),"3",as.character(sleepstudyM$Subject))
lme2 <- lmer(Reaction ~ Days + (1 | Subject) + (0 + Days | Subject), data=sleepstudyM, REML=FALSE)
# the dropped row should cause a warning
expect_warning(wm2 <- mix(Reaction ~ Days + (1 | Subject) + (0 + Days | Subject), data=sleepstudyM, weights=c("weight1L1", "weight1L2")),"with missing data")
expect_equal(wm2$lnl, as.numeric(logLik(lme2)), tol=2E-7*abs(as.numeric(logLik(lme2))))
# check coef
expect_equal(coef(wm2),
expected = getME(lme2, "fixef"),
tolerance = tolerance)
# check vars
lmewm2vars <- data.frame(summary(lme2)$varcor)$sdcor
expect_equal(unname(wm2$vars),
expected = unname(lmewm2vars)^2,
tolerance = tolerance)
lme4ranef <- ranef(lme2)$Subject
# not in WeMix output
attr(lme4ranef, "postVar") <- NULL
expect_equal(lme4ranef, wm2$ranefMat$Subject, 0.01)
expect_equal(predict(lme2), predict(wm2), tol=1e-4)
sleepstudyM2 <- sleepstudyM[order(runif(nrow(sleepstudyM))), ]
sleepstudyM2$Subject <- ifelse(sleepstudyM2$Subject %in% c("332"),"4",as.character(sleepstudyM$Subject))
sleepstudyM2$Subject <- ifelse(sleepstudyM2$Subject %in% c("308"),"1",as.character(sleepstudyM$Subject))
sleepstudyM2$Subject <- ifelse(sleepstudyM2$Subject %in% c("351"),"5",as.character(sleepstudyM$Subject))
expect_equal(predict(lme2,sleepstudyM2, allow.new.levels=TRUE), predict(wm2,sleepstudyM2, allow.new.levels=TRUE), tol=1e-4)
})
context("Mean Centering Matches HLM results")
test_that("Mean Centering Matches HLM results", {
wm1 <- mix(Reaction ~ Days + (1|Subject), data=sleepstudyU, weights=c("weight1L1", "weight1L2"), nQuad=13,center_group=list("Subject"= as.formula(~Days)), verbose=FALSE)
expect_equal(unname(wm1$coef), c(298.507892, 10.467286), tolerance=1E-1)
expect_equal(wm1$lnl, -897.0393, tolerance=tolerance*897) # value from lmer(Reaction ~ Days + (1 | Subject), data=sleepstudy, REML=FALSE)
})
ss <- sleepstudy
ss1 <- ss
ss2 <- ss
doubles <- c(308, 309, 310) # subject with double obs
ss2 <- rbind(ss, subset(ss, Subject %in% doubles))
ss$w <- ifelse(ss$Subject %in% doubles, 2, 1)
contrasts(ss2$Subject) <- "contr.sum"
ss1$W1 <- ifelse(ss1$Subject %in% doubles, 2, 1)
ss1$W2 <- 1
ss1$bin <- ifelse(sleepstudy$Reaction<300,0,1) #for the binomial test
ss2$W2 <- ss2$W1 <- 1
doubles <- c(308, 309, 310) # subject with double obs
ss30 <- subset(ss, Subject %in% doubles)
ss30$Subject <- as.numeric(as.character(ss30$Subject)) + 1000
ss0 <- ss
ss0$Subject <- as.numeric(as.character(ss$Subject))
ss3 <- rbind(ss0, ss30)
ss3$Subject <- as.factor(ss3$Subject)
ss3$W2 <- 1
ss3$W1 <- 1
ss4 <- ss
ss4$W2 <- ifelse(ss4$Subject %in% doubles, 2, 1)
ss4$W1 <- ifelse(ss4$Subject %in% doubles, 2, 1) # make unconditional weight
context("GLM works: Binomial")
test_that("GLM works: Binomial", {
skip_on_cran()
#full test for binomial
bi_1 <- mix(bin~Days + (1|Subject), data=ss1, family=binomial(link="logit"), verbose=FALSE,
weights=c("W1", "W2"), nQuad=13)
expect_equal(unname(bi_1$coef), c(-3.3448,.5928), tolerance=1E-3)
expect_equal(as.numeric(bi_1$lnl), -93.751679, tolerance=1E-5)
sum_bi <- summary(bi_1)
expect_is(summary(bi_1), "summaryWeMixResults")
})
context("Repeating is the same as weighting: L1 replicate vs weighting")
test_that("Repeating is the same as weighting: L1 replicate vs weighting", {
# mix for L1, weighted
wmeL1W <- mix(formula=Reaction ~ Days + (1 | Subject), data=ss1,
weights=c("W1", "W2"))
# mix for L1, duplicated
system.time(wmeL1D <- mix(formula=Reaction ~ Days + (1 | Subject), data=ss2,
weights=c("W1", "W2")))
# check weighted agrees with duplicated lme4 results
expect_equal(wmeL1W$lnl, -1048.34318418762, tolerance=1050*tolerance)
# check duplicated agrees with duplicated lme4 results
expect_equal(wmeL1D$lnl, -1048.34318418762, tolerance=1050*tolerance)
})
context("grouping factor not sorted")
test_that("grouping factor not sorted", {
skip_on_cran()
ss1_mixed <- ss1[c(125:180,1,100,2:99,101:124),]
row.names(ss1_mixed) <- NULL
# mix for L1, weighted
wmeL1W <- mix(formula=Reaction ~ Days + (1 | Subject), data=ss1_mixed,
weights=c("W1", "W2"), nQuad=13, run=TRUE, verbose=FALSE)
# mix for L1, duplicated
system.time(wmeL1D <- mix(formula=Reaction ~ Days + (1 | Subject), data=ss2,
weights=c("W1", "W2"), nQuad=13, run=FALSE, verbose=FALSE))
#statares <- c(251.4619, 10.40726, 1000.7466, 1338.0865) # not used
# check weighted agrees with duplicated lme4 results
expect_equal(wmeL1W$lnl, -1048.34318418762, tolerance=tolerance)
expect_equal(wmeL1D$lnl, -1048.34318418762, tolerance=tolerance)
# check final results
suppressWarnings(mix1 <- mix(formula=Reaction ~ Days + (1 | Subject), data=ss1_mixed,
weights=c("W1", "W2")))
suppressWarnings(mix1REF <- mix(formula=Reaction ~ Days + (1 | Subject), data=ss1,
weights=c("W1", "W2")))
expect_equal(mix1$coef, mix1REF$coef, tolerance=1e3)
expect_equal(mix1$vars, mix1REF$vars, tolerance=1e-3)
expect_equal(mix1$lnl, mix1REF$lnl, tolerance=1e-3)
#check weighted fixed effects variances
expect_equal(unname(sqrt(diag(mix1$cov_mat))), unname(sqrt(diag(mix1REF$cov_mat))),tolerance = tolerance)
})
context("Weighted three level model unsorted")
test_that("Weighted three level model unsorted", {
skip_on_cran()
sleepstudy2 <- sleepstudy
sleepstudy2$Group <- 1
sleepstudy2$Group <- ifelse(sleepstudy2$Subject %in% c("310", "309", "349", "335"), 2, sleepstudy2$Group)
sleepstudy2$Group <- ifelse(sleepstudy2$Subject %in% c("330", "352", "337", "369"), 3, sleepstudy2$Group)
sleepstudy2$Group <- ifelse(sleepstudy2$Subject %in% c("331", "332", "334", "372"), 4, sleepstudy2$Group)
sleepstudy2$Group <- factor(sleepstudy2$Group)
ss2 <- sleepstudy2
w1c <- w1 <- rep(1,180)
w2c <- w2 <- rep(1,18)
w3c <- w3 <- rep(1,4)
# unbalanced (non-identical within a Subject), level-1 (obs level) weights
w1c[1:4] <- w1[1:4] <- 2
uR <- sleepstudy2[1:4,]
sleepstudy2 <- rbind(sleepstudy2, uR)
# level-2 weights
w2c[1] <- w2[1] <- 2
w1[ss2$Subject=="308"] <- 2*w1[ss2$Subject=="308"]
sR <- subset(sleepstudy2, Subject == "308")
sR$Subject <- "S308"
sleepstudy2 <- rbind(sleepstudy2, sR)
# level-3 weights
w3c[1] <- w3[1] <- 2
w2[c(1,7,12,13,16,17)] <- 2*w2[c(1,7,12,13,16,17)]
w1[ss2$Group==1] <- 2*w1[ss2$Group==1]
gR <- subset(sleepstudy2, Group %in% 1)
gR$Subject <- paste0("G", gR$Subject)
gR$Group <- "11"
sleepstudy2 <- rbind(sleepstudy2, gR)
# lmr for reference
lmr <- lmer(Reaction ~ Days + (1|Subject) + (1|Group),
data=sleepstudy2, control=lmerControl(optimizer="bobyqa"),
REML=FALSE)
ss2$w1 <- w1
ss2$w2 <- rep(w2,each=10)
ss2$w3 <- ifelse(ss2$Subject %in% c("308", "333", "350", "351", "370", "371"),2,1)
ss3 <- ss2[sample(row.names(ss2), size=nrow(ss2)), ]
wm0 <- mix(Reaction ~ Days + (1|Subject) + (1|Group), data=ss3,
weights=c("w1", "w2", "w3"))
# check coef
expect_equal(coef(wm0),
expected = getME(lmr, "fixef"),
tolerance = tolerance)
lmeRanef <- ranef(lmr)$Subject
# not in WeMix output
attr(lmeRanef, "postVar") <- NULL
# only first 18 of lmeRanef are unique/comparable
expect_equal(lmeRanef[1:18,,drop=FALSE], wm0$ranefMat$Subject, 20*sqrt(.Machine$double.eps))
# check vars
lmewm2vars <- data.frame(summary(lmr)$varcor)$sdcor
expect_equal(unname(wm0$vars),
expected = unname(lmewm2vars)^2,
tolerance = tolerance)
# var beta not expected to be equal
})
context("repeating is the same as weighting: L2 replicate vs weighting")
test_that("L2 replicate vs weighting", {
# mix for L2, duplicated
system.time(wmeL2D <- mix(formula=Reaction ~ Days + (1 | Subject),
data=ss3, weights=c("W1", "W2")))
# mix for L2, weighted
system.time(wmeL2W <- mix(formula=Reaction ~ Days + (1 | Subject), data=ss4,
weights=c("W1", "W2")))
expect_equal(wmeL2W$lnl, -1055.34690957995, tolerance=2E-7)
expect_equal(wmeL2D$lnl, -1055.34690957995, tolerance=2E-7)
})
context("Repeating is the same as weighting: L1 replicate vs weighting, 2 REs")
test_that("Repeating is the same as weighting: L1 replicate vs weighting, 2 REs", {
# mix for L1, weighted, 2 REs
wmeL1WRE2 <- mix(formula=Reaction ~ Days + (1 | Subject) + (0+Days|Subject),
data=ss1, weights=c("W1", "W2"), nQuad=13, run=FALSE, verbose=FALSE)
# mix for L1, duplicated, 2 REs
wmeL1DRE2 <- mix(formula=Reaction ~ Days + (1 | Subject) + (0+Days|Subject),
data=ss2, weights=c("W1", "W2"),nQuad=13, run=FALSE, verbose=FALSE)
expect_equal(wmeL1WRE2$lnl, -1018.29298875158, tolerance=1050*2E-7)
expect_equal(wmeL1DRE2$lnl, -1018.29298875158, tolerance=1050*2E-7)
})
ssB <- sleepstudy
set.seed(2)
ssB$Reaction <- ssB$Days * 3.141 + rnorm(nrow(ssB))
ssB$W2 <- 1
ssB$W1 <- 1
context("Zero variance estimate")
test_that("Zero variance estimate", {
skip_on_cran()
# this has 0 variance estimate in lmer
lmeB <- lmer(Reaction ~ Days + (1|Subject), data=ssB, REML=FALSE)
mixB <- mix(formula=Reaction ~ Days + (1 | Subject), data=ssB,
weights=c("W1", "W2"), nQuad=13, run=TRUE, verbose=FALSE)
expect_equal(mixB$lnl, as.numeric(logLik(lmeB)), tol=1e-5)
expect_equal(coef(mixB), fixef(lmeB), tol=1e-5)
expect_equal(unname(mixB$vars[length(mixB$vars)]), unname(lmeB@devcomp$cmp["sigmaML"]^2), tol=1e-5)
ss1 <- sleepstudy
#add group variables for 3 level model
ss1$Group <- 1
ss1$Group <- ifelse(ss1$Subject %in% c(349,335,330, 352, 337, 369), 2, ss1$Group)
# Create weights
ss1$W1 <- ifelse(ss1$Subject %in% c(308, 309, 310), 2, 1)
ss1$W2 <- 1
ss1$W3 <- ifelse(ss1$Group == 2, 2, 1)
#Run three level model with random slope and intercept.
three_level <- mix(Reaction ~ Days + (1|Subject) + (1+Days|Group), data=ss1, weights = c("W1","W2","W3"))
expect_is(three_level, "WeMixResults")
})
context("Unweighted three level model")
test_that("Unweighted three level model", {
sleepstudy2 <- sleepstudy
sleepstudy2$Group <- 1
sleepstudy2$Group <- ifelse(sleepstudy2$Subject %in% c("310", "309", "349", "335"), 2, sleepstudy2$Group)
sleepstudy2$Group <- ifelse(sleepstudy2$Subject %in% c("330", "352", "337", "369"), 3, sleepstudy2$Group)
sleepstudy2$Group <- ifelse(sleepstudy2$Subject %in% c("331", "332", "334", "372"), 4, sleepstudy2$Group)
sleepstudy2$Group <- factor(sleepstudy2$Group)
sleepstudy2$w1 <- 1
sleepstudy2$w2 <- 1
sleepstudy2$w3 <- 1
wm0 <- mix(Reaction ~ Days + (1|Subject) + (0+Days|Subject) + (1 | Group), data=sleepstudy2, weights=c("w1", "w2","w3"), verbose=FALSE, run=TRUE)
lm0 <- lmer(Reaction ~ Days + (1|Subject) + (0+Days|Subject) + (1 | Group), data=sleepstudy2,REML=FALSE)
wm0RE <- wm0$ranefMat
lm0RE <- ranef(lm0)
expect_equal(names(wm0RE), names(lm0RE))
attr(lm0RE$Subject, "postVar") <- NULL
expect_equal(wm0RE$Subject, lm0RE$Subject, tol=.Machine$double.eps^0.25)
attr(lm0RE$Group, "postVar") <- NULL
expect_equal(wm0RE$Group, lm0RE$Group, tol=.Machine$double.eps^0.25)
wm0RE <- wm0$ranefMat
lm0RE <- ranef(lm0)
attr(lm0RE$Subject, "postVar") <- NULL
expect_equal(wm0RE$Subject, lm0RE$Subject, tol=.Machine$double.eps^0.25)
# check vars
lmevars1 <- data.frame(summary(lm0)$varcor)$sdcor
expect_equal(unname(wm0$vars),
expected = unname(lmevars1)^2,
tolerance = 1e-3)
expect_equal(wm0$lnl, as.numeric(logLik(lm0)), tol=1e-3)
expect_equal(coef(wm0), fixef(lm0), tol=1e-4)
expect_equal(unname(wm0$vars[length(wm0$vars)]), unname(lm0@devcomp$cmp["sigmaML"]^2), tol=1e-4)
ranef <- ranef(lm0)
# WeMix does not have postVar attribute
attr(ranef$Subject, "postVar") <- NULL
attr(ranef$Group, "postVar") <- NULL
expect_equal(ranef$Subject, wm0$ranefMat$Subject, (.Machine$double.eps)^0.25)
expect_equal(ranef$Group, wm0$ranefMat$Group, (.Machine$double.eps)^0.25)
})
context("Unweighted three level model, binomial")
test_that("Unweighted three level model, binomial", {
sleepstudy2 <- sleepstudy
set.seed(564)
sleepstudy2$rand_slope <- rnorm(180,mean=3.25,sd=1.79)
sleepstudy2$highR <- ifelse(sleepstudy2$Reaction > 340, 1, 0)
sleepstudy2$Group <- 1
sleepstudy2$Group <- ifelse(sleepstudy2$Subject %in% c("310", "309", "349", "335"), 2, sleepstudy2$Group)
sleepstudy2$Group <- ifelse(sleepstudy2$Subject %in% c("330", "352", "337", "369"), 3, sleepstudy2$Group)
sleepstudy2$Group <- ifelse(sleepstudy2$Subject %in% c("331", "332", "334", "372"), 4, sleepstudy2$Group)
sleepstudy2$Group <- factor(sleepstudy2$Group)
sleepstudy2$w1 <- 1
sleepstudy2$w2 <- 1
sleepstudy2$w3 <- 1
sleepstudy2 <- sleepstudy2[sample(1:nrow(sleepstudy2), nrow(sleepstudy2)),]
# up the max_iterations for more complex model structure
suppressWarnings(wm0 <- mix(highR ~ Days + (rand_slope|Subject) + (1 | Group), data=sleepstudy2,
family=binomial(),weights=c("w1", "w2","w3"), verbose=FALSE, run=TRUE,
nQuad=1,max_iteration = 20))
lm0 <- glmer(highR ~ Days + (rand_slope|Subject) + (1 | Group), data=sleepstudy2,
family=binomial())
# check vars
lmevars1 <- data.frame(summary(lm0)$varcor)$vcov
expect_equal(unname(wm0$vars),
expected = unname(lmevars1),
tolerance = 1e-3)
expect_equal(as.numeric(wm0$lnl), as.numeric(logLik(lm0)), tol=1e-3)
expect_equal(coef(wm0), fixef(lm0), tol=1e-4)
expect_equal(predict(wm0), predict(lm0), tol=1e-4)
expect_equal(predict(wm0, sleepstudy2), predict(lm0, sleepstudy2), tol=1e-4)
sleepstudy3 <- sleepstudy2[sample(1:nrow(sleepstudy2), floor(nrow(sleepstudy2)/4)),]
expect_equal(a1 <- predict(wm0, sleepstudy3), b1 <- predict(lm0, sleepstudy3), tol=1e-4)
sleepstudy4 <- sleepstudy2
sleepstudy4$Subject <- ifelse(sleepstudy4$Subject == "332", "1", as.character(sleepstudy4$Subject))
expect_equal(a2 <- predict(wm0, sleepstudy4, allow.new.levels=TRUE), b2 <- predict(lm0, sleepstudy4, allow.new.levels=TRUE), tol=1e-4)
})
context("Weighted v unweighted replicated two level model")
test_that("Weighted v unweighted replicated two level model", {
sleepstudy2 <- sleepstudy
sleepstudy2$w1 <- 1
sleepstudy2$w2 <- 1
sleepstudy2$w3 <- 1
sleepstudy2$Subject <- as.character(sleepstudy2$Subject)
g1 <- subset(sleepstudy2, sleepstudy2$Subject %in% c("310", "309", "349", "335"))
g1$Subject <- paste0("R", g1$Subject)
sleepstudyrep <- rbind(sleepstudy2, g1)
sleepstudy2$Subject <- factor(sleepstudy2$Subject)
sleepstudy2$w2 <- ifelse(sleepstudy2$Subject %in% c("310", "309", "349", "335"), 2, 1)
# two level
wm0 <- mix(Reaction ~ Days + (1|Subject), data=sleepstudyrep, weights=c("w1", "w2"), verbose=FALSE, run=TRUE)
wmw <- mix(Reaction ~ Days + (1|Subject), data=sleepstudy2, weights=c("w1", "w2"), verbose=FALSE, run=TRUE, cWeights=TRUE)
lm0 <- lmer(Reaction ~ Days + (1|Subject), data=sleepstudyrep, REML=FALSE)
# check vars
lmevars1 <- data.frame(summary(lm0)$varcor)$sdcor
expect_equal(unname(wm0$vars),
expected = unname(lmevars1)^2,
tolerance = 1e-3)
expect_equal(wm0$lnl, as.numeric(logLik(lm0)), tol=1e-3)
expect_equal(coef(wm0), fixef(lm0), tol=1e-4)
expect_equal(unname(wm0$vars[length(wm0$vars)]), unname(lm0@devcomp$cmp["sigmaML"]^2), tol=1e-4)
lm0ranef <- ranef(lm0)$Subject
# not in WeMix output
attr(lm0ranef, "postVar") <- NULL
expect_equal(lm0ranef, wm0$ranefMat$Subject, 20*sqrt(.Machine$double.eps))
expect_equal(wmw$lnl, wm0$lnl, tol=1e-3)
expect_equal(coef(wmw), coef(wm0), tol=1e-4)
expect_equal(wmw$vars, wm0$vars, tol=1e-4)
# should be the same if not conditional as well
sleepstudy2$w1 <- sleepstudy2$w2
wmw <- mix(Reaction ~ Days + (1|Subject), data=sleepstudy2, weights=c("w1", "w2"), verbose=FALSE, run=TRUE, cWeights=FALSE)
expect_equal(wmw$lnl, wm0$lnl, tol=1e-3)
expect_equal(coef(wmw), coef(wm0), tol=1e-4)
expect_equal(wmw$vars, wm0$vars, tol=1e-4)
# sub 1 conditional weights work
sleepstudy2 <- sleepstudy
sleepstudy2$w1 <- 1
sleepstudy2$w2 <- 1
sleepstudy2$w3 <- 1
sleepstudy2$Subject <- as.character(sleepstudy2$Subject)
g1 <- subset(sleepstudy2, sleepstudy2$Subject %in% c("310", "309", "349", "335"))
# do not renaem subject
sleepstudyrep <- rbind(sleepstudy2, g1)
sleepstudyrep$w1 <- ifelse(sleepstudyrep$Subject %in% c("310", "309", "349", "335"), 0.5, 1)
# two level
wm0 <- mix(Reaction ~ Days + (1|Subject), data=sleepstudyrep, weights=c("w1", "w2"), verbose=FALSE, run=TRUE, cWeights=TRUE)
wmw <- mix(Reaction ~ Days + (1|Subject), data=sleepstudy2, weights=c("w1", "w2"), verbose=FALSE, run=TRUE)
lm0 <- lmer(Reaction ~ Days + (1|Subject), data=sleepstudy2, REML=FALSE)
# check vars
lmevars1 <- data.frame(summary(lm0)$varcor)$sdcor
expect_equal(unname(wm0$vars),
expected = unname(lmevars1)^2,
tolerance = 1e-3)
expect_equal(wm0$lnl, as.numeric(logLik(lm0)), tol=1e-3)
expect_equal(coef(wm0), fixef(lm0), tol=1e-4)
expect_equal(unname(wm0$vars[length(wm0$vars)]), unname(lm0@devcomp$cmp["sigmaML"]^2), tol=1e-4)
expect_equal(wmw$lnl, wm0$lnl, tol=1e-3)
expect_equal(coef(wmw), coef(wm0), tol=1e-4)
expect_equal(wmw$vars, wm0$vars, tol=1e-4)
# should be the same if not conditional as well
wm0 <- mix(Reaction ~ Days + (1|Subject), data=sleepstudyrep, weights=c("w1", "w2"), verbose=FALSE, run=TRUE)
expect_equal(wmw$lnl, wm0$lnl, tol=1e-3)
expect_equal(coef(wmw), coef(wm0), tol=1e-4)
expect_equal(wmw$vars, wm0$vars, tol=1e-4)
})
context("Weighted v unweighted replicated three level model")
test_that("Weighted v unweighted replicated three level model", {
sleepstudy2 <- sleepstudy
sleepstudy2$Group <- 1
sleepstudy2$Group <- ifelse(sleepstudy2$Subject %in% c("310", "309", "349", "335"), 2, sleepstudy2$Group)
sleepstudy2$Group <- ifelse(sleepstudy2$Subject %in% c("330", "352", "337", "369"), 3, sleepstudy2$Group)
sleepstudy2$Group <- ifelse(sleepstudy2$Subject %in% c("331", "332", "334", "372"), 4, sleepstudy2$Group)
#sleepstudy2$Group <- factor(sleepstudy2$Group)
sleepstudy2$w1 <- 1
sleepstudy2$w2 <- 1
sleepstudy2$w3 <- 1
sleepstudy2$Subject <- as.character(sleepstudy2$Subject)
g1 <- subset(sleepstudy2, Group == 1)
g1$Group <- 5
g1$Subject <- paste0("5", g1$Subject)
sleepstudyrep <- rbind(sleepstudy2, g1)
s308 <- subset(sleepstudy2, Subject %in% "308")
s308$Subject <- "R308"
sleepstudyrep <- rbind(sleepstudyrep, s308)
s308$Group <- 5
s308$Subject <- "5R308"
sleepstudyrep <- rbind(sleepstudyrep, s308)
sleepstudy2$w1u <- ifelse(sleepstudy2$Group == 1, 2, 1)
s308unit <- subset(sleepstudyrep, Subject %in% c("308", "5308", "R308", "5R308") & Days %in% c(0,3))
sleepstudyrep <- rbind(sleepstudyrep, s308unit)
sleepstudy2$Group <- factor(sleepstudy2$Group)
sleepstudy2$Subject <- factor(sleepstudy2$Subject)
sleepstudy2$w1 <- ifelse(sleepstudy2$Subject %in% c("308", "R308", "5R308") & sleepstudy2$Days %in% c(0,3), 2, 1)
sleepstudy2$w2 <- ifelse(sleepstudy2$Subject == "308", 2, 1)
sleepstudy2$w3 <- ifelse(sleepstudy2$Group == 1, 2, 1)
sleepstudy2$w2u <- ifelse(sleepstudy2$Group == 1, 2, 1)
sleepstudy2$w2u <- ifelse(sleepstudy2$Subject == "308", 4, sleepstudy2$w2u)
sleepstudy2$w1u <- ifelse(sleepstudy2$Group == 1, 2, 1)
sleepstudy2$w1u <- ifelse(sleepstudy2$Subject == "308", 4, sleepstudy2$w1u)
sleepstudy2$w1u <- ifelse(sleepstudy2$Subject == "308" & sleepstudy2$Days %in% c(0,3), 8, sleepstudy2$w1u)
# replication
wm0 <- mix(Reaction ~ Days + (1|Subject) + (1 | Group), data=sleepstudyrep, weights=c("w1", "w2","w3"), verbose=FALSE, run=TRUE)
# conditional weights
wmw <- mix(Reaction ~ Days + (1|Subject) + (1 | Group), data=sleepstudy2, weights=c("w1", "w2","w3"), verbose=FALSE, run=TRUE, cWeights=TRUE)
expect_equal(wm0$lnl, wmw$lnl, tol=20*sqrt(.Machine$double.eps))
expect_equal(coef(wmw), coef(wm0), tol=20*sqrt(.Machine$double.eps))
# unconditional weights
wmwc <- mix(Reaction ~ Days + (1|Subject) + (1 | Group), data=sleepstudy2, weights=c("w1u", "w2u","w3"), verbose=FALSE, run=TRUE, cWeights=FALSE)
expect_equal(wm0$lnl, wmwc$lnl, tol=20*sqrt(.Machine$double.eps))
expect_equal(coef(wmwc), coef(wm0), tol=20*sqrt(.Machine$double.eps))
lm0 <- lmer(Reaction ~ Days + (1|Subject) + (1 | Group), data=sleepstudyrep,REML=FALSE)
# check vars
lmevars1 <- data.frame(summary(lm0)$varcor)$sdcor
expect_equal(unname(wm0$vars),
expected = unname(lmevars1)^2,
tolerance = 1e-3)
expect_equal(wm0$lnl, as.numeric(logLik(lm0)), tol=1e-3)
expect_equal(coef(wm0), fixef(lm0), tol=1e-4)
expect_equal(unname(wm0$vars[length(wm0$vars)]), unname(lm0@devcomp$cmp["sigmaML"]^2), tol=1e-4)
expect_equal(wmw$lnl, wm0$lnl, tol=1e-3)
expect_equal(coef(wmw), coef(wm0), tol=1e-4)
expect_equal(wmw$vars, wm0$vars, tol=1e-4)
expect_equal(wm0$ranefMat$Subject[1:18,,drop=FALSE], wmw$ranefMat$Subject, 20*sqrt(.Machine$double.eps))
expect_equal(wm0$ranefMat$Group[1:4,,drop=FALSE], wmw$ranefMat$Group, 20*sqrt(.Machine$double.eps))
})
context("Weighted v unweighted replicated two level model, binomial")
test_that("Weighted v unweighted replicated two level model, binomial", {
sleepstudy2 <- sleepstudy
sleepstudy2$w1 <- 1
sleepstudy2$w2 <- 1
sleepstudy2$w3 <- 1
sleepstudy2$Subject <- as.character(sleepstudy2$Subject)
sleepstudy2$highR <- ifelse(sleepstudy2$Reaction > 340, 1, 0)
g1 <- subset(sleepstudy2, sleepstudy2$Subject %in% c("310", "309", "349", "335"))
g1$Subject <- paste0("R", g1$Subject)
sleepstudyrep <- rbind(sleepstudy2, g1)
sleepstudy2$Subject <- factor(sleepstudy2$Subject)
sleepstudy2$w2 <- ifelse(sleepstudy2$Subject %in% c("310", "309", "349", "335"), 2, 1)
# two level
wm0 <- mix(highR ~ Days + (1|Subject), data=sleepstudyrep, weights=c("w1", "w2"),
family=binomial(), verbose=FALSE)
wmw <- mix(highR ~ Days + (1|Subject), data=sleepstudy2, weights=c("w1", "w2"),
family=binomial(), verbose=FALSE, cWeights=TRUE,
start=c(wm0$coef,wm0$vars))
lm0 <- glmer(highR ~ Days + (1|Subject), data=sleepstudyrep, family=binomial(),
nAGQ = 13)
# check vars
lmevars1 <- data.frame(summary(lm0)$varcor)$sdcor
expect_equal(unname(wm0$vars),
expected = unname(lmevars1)^2,
tolerance = 1e-3)
expect_equal(as.numeric(wm0$lnl), as.numeric(logLik(lm0)), tol=1e-3)
expect_equal(coef(wm0), fixef(lm0), tol=1e-4)
expect_equal(wmw$lnl, wm0$lnl, tol=1e-3)
expect_equal(coef(wmw), coef(wm0), tol=1e-4)
expect_equal(wmw$vars, wm0$vars, tol=1e-4)
# should be the same if not conditional as well
sleepstudy2$w1 <- sleepstudy2$w2
wmw <- mix(highR ~ Days + (1|Subject), data=sleepstudy2, weights=c("w1", "w2"),
family=binomial(),
cWeights=FALSE)
expect_equal(wmw$lnl, wm0$lnl, tol=1e-3)
expect_equal(coef(wmw), coef(wm0), tol=1e-4)
expect_equal(wmw$vars, wm0$vars, tol=1e-4)
# sub 1 conditional weights work
sleepstudy2 <- sleepstudy
sleepstudy2$highR <- ifelse(sleepstudy2$Reaction > 340, 1, 0)
sleepstudy2$w1 <- 1
sleepstudy2$w2 <- 1
sleepstudy2$w3 <- 1
sleepstudy2$Subject <- as.character(sleepstudy2$Subject)
g1 <- subset(sleepstudy2, sleepstudy2$Subject %in% c("310", "309", "349", "335"))
# do not renaem subject
sleepstudyrep <- rbind(sleepstudy2, g1)
sleepstudyrep$w1 <- ifelse(sleepstudyrep$Subject %in% c("310", "309", "349", "335"), 0.5, 1)
# two level
wm0 <- mix(highR ~ Days + (1|Subject), data=sleepstudyrep, weights=c("w1", "w2"),
family=binomial(),verbose=FALSE, run=TRUE, cWeights=TRUE)
wmw <- mix(highR ~ Days + (1|Subject), data=sleepstudy2, weights=c("w1", "w2"),
family=binomial(),verbose=FALSE, run=TRUE)
lm0 <- glmer(highR ~ Days + (1|Subject), data=sleepstudy2, family=binomial(),
nAGQ=13)
# check vars
lmevars1 <- data.frame(summary(lm0)$varcor)$sdcor
expect_equal(unname(wm0$vars),
expected = unname(lmevars1)^2,
tolerance = 1e-3)
expect_equal(as.numeric(wm0$lnl), as.numeric(logLik(lm0)), tol=1e-3)
expect_equal(coef(wm0), fixef(lm0), tol=1e-4)
expect_equal(wmw$lnl, wm0$lnl, tol=1e-3)
expect_equal(coef(wmw), coef(wm0), tol=1e-4)
expect_equal(wmw$vars, wm0$vars, tol=1e-4)
# should be the same if not conditional as well
wm0 <- mix(highR ~ Days + (1|Subject), data=sleepstudyrep, weights=c("w1", "w2"),
family=binomial(),verbose=FALSE, run=TRUE)
expect_equal(wmw$lnl, wm0$lnl, tol=1e-3)
expect_equal(coef(wmw), coef(wm0), tol=1e-4)
expect_equal(wmw$vars, wm0$vars, tol=1e-4)
})
context("Weighted v unweighted replicated two level model, binomial with probit link")
test_that("Weighted v unweighted replicated two level model, binomial with probit link", {
sleepstudy2 <- sleepstudy
sleepstudy2$w1 <- 1
sleepstudy2$w2 <- 1
sleepstudy2$w3 <- 1
sleepstudy2$Subject <- as.character(sleepstudy2$Subject)
sleepstudy2$highR <- ifelse(sleepstudy2$Reaction > 340, 1, 0)
g1 <- subset(sleepstudy2, sleepstudy2$Subject %in% c("310", "309", "349", "335"))
g1$Subject <- paste0("R", g1$Subject)
sleepstudyrep <- rbind(sleepstudy2, g1)
sleepstudy2$Subject <- factor(sleepstudy2$Subject)
sleepstudy2$w2 <- ifelse(sleepstudy2$Subject %in% c("310", "309", "349", "335"), 2, 1)
# two level
suppressWarnings(wm0 <- mix(highR ~ Days + (1|Subject), data=sleepstudyrep, weights=c("w1", "w2"),
family=binomial(link="probit"), verbose=FALSE))
suppressWarnings(wmw <- mix(highR ~ Days + (1|Subject), data=sleepstudy2, weights=c("w1", "w2"),
family=binomial(link="probit"), verbose=FALSE, cWeights=TRUE,
start=c(wm0$coef,wm0$vars)))
lm0 <- glmer(highR ~ Days + (1|Subject), data=sleepstudyrep, family=binomial(link="probit"),
nAGQ = 13)
# check vars
lmevars1 <- data.frame(summary(lm0)$varcor)$sdcor
expect_equal(unname(wm0$vars),
expected = unname(lmevars1)^2,
tolerance = 1e-3)
expect_equal(as.numeric(wm0$lnl), as.numeric(logLik(lm0)), tol=1e-3)
expect_equal(coef(wm0), fixef(lm0), tol=1e-4)
expect_equal(wmw$lnl, wm0$lnl, tol=1e-3)
expect_equal(coef(wmw), coef(wm0), tol=1e-4)
expect_equal(wmw$vars, wm0$vars, tol=1e-4)
# should be the same if not conditional as well
sleepstudy2$w1 <- sleepstudy2$w2
suppressWarnings(wmw <- mix(highR ~ Days + (1|Subject), data=sleepstudy2, weights=c("w1", "w2"),
family=binomial(link="probit"),
cWeights=FALSE))
expect_equal(wmw$lnl, wm0$lnl, tol=1e-3)
expect_equal(coef(wmw), coef(wm0), tol=1e-4)
expect_equal(wmw$vars, wm0$vars, tol=1e-4)
# sub 1 conditional weights work
sleepstudy2 <- sleepstudy
sleepstudy2$highR <- ifelse(sleepstudy2$Reaction > 340, 1, 0)
sleepstudy2$w1 <- 1
sleepstudy2$w2 <- 1
sleepstudy2$w3 <- 1
sleepstudy2$Subject <- as.character(sleepstudy2$Subject)
g1 <- subset(sleepstudy2, sleepstudy2$Subject %in% c("310", "309", "349", "335"))
# do not renaem subject
sleepstudyrep <- rbind(sleepstudy2, g1)
sleepstudyrep$w1 <- ifelse(sleepstudyrep$Subject %in% c("310", "309", "349", "335"), 0.5, 1)
# two level
suppressWarnings(wm0 <- mix(highR ~ Days + (1|Subject), data=sleepstudyrep, weights=c("w1", "w2"),
family=binomial(link="probit"),verbose=FALSE, run=TRUE, cWeights=TRUE))
suppressWarnings(wmw <- mix(highR ~ Days + (1|Subject), data=sleepstudy2, weights=c("w1", "w2"),
family=binomial(link="probit"),verbose=FALSE, run=TRUE))
lm0 <- glmer(highR ~ Days + (1|Subject), data=sleepstudy2, family=binomial(link="probit"),
nAGQ=13)
# check vars
lmevars1 <- data.frame(summary(lm0)$varcor)$sdcor
expect_equal(unname(wm0$vars),
expected = unname(lmevars1)^2,
tolerance = 1e-3)
expect_equal(as.numeric(wm0$lnl), as.numeric(logLik(lm0)), tol=1e-3)
expect_equal(coef(wm0), fixef(lm0), tol=1e-4)
expect_equal(wmw$lnl, wm0$lnl, tol=1e-3)
expect_equal(coef(wmw), coef(wm0), tol=1e-4)
expect_equal(wmw$vars, wm0$vars, tol=1e-4)
# should be the same if not conditional as well
suppressWarnings(wm0 <- mix(highR ~ Days + (1|Subject), data=sleepstudyrep, weights=c("w1", "w2"),
family=binomial(link="probit"),verbose=FALSE, run=TRUE))
expect_equal(wmw$lnl, wm0$lnl, tol=1e-3)
expect_equal(coef(wmw), coef(wm0), tol=1e-4)
expect_equal(wmw$vars, wm0$vars, tol=1e-4)
})
context("Weighted v unweighted replicated two level model, poisson")
test_that("Weighted v unweighted replicated two level model, poisson", {
skip_on_cran()
skip_if_not_installed("glmmTMB")
require(glmmTMB)
owls2 <- Owls
owls2$w1 <- 1
owls2$w2 <- 1
owls2$w3 <- 1
g1 <- subset(owls2, owls2$Nest %in% c("Jeuss", "Lully", "Seiry", "Oleyes"))
g1$Nest <- paste0("R", g1$Nest)
owlsrep <- rbind(owls2, g1)
owls2$w2 <- ifelse(owls2$Nest %in% c("Jeuss", "Lully", "Seiry", "Oleyes"), 2, 1)
# glmer seems to have a bug with nAGQ > 1 for poisson models
# (see https://stat.ethz.ch/pipermail/r-sig-mixed-models/2020q2/028606.html),
# so can only compare with nQuad = 1
suppressWarnings(wm0 <- mix(SiblingNegotiation ~ ArrivalTime + (1|Nest), data=owlsrep, weights=c("w1", "w2"),
family=poisson(), verbose=FALSE,nQuad=1))
suppressWarnings(wmw <- mix(SiblingNegotiation ~ ArrivalTime + (1|Nest), data=owls2, weights=c("w1", "w2"), nQuad = 1,
family=poisson(), verbose=FALSE, cWeights=TRUE))
lm0 <- glmer(SiblingNegotiation ~ ArrivalTime + (1|Nest), data=owlsrep,
family="poisson")
# check vars
lmevars1 <- data.frame(summary(lm0)$varcor)$sdcor
expect_equal(unname(wm0$vars),
expected = unname(lmevars1)^2,
tolerance = 1e-3)
expect_equal(as.numeric(wm0$lnl), as.numeric(logLik(lm0)), tol=1e-3)
expect_equal(coef(wm0), fixef(lm0), tol=1e-4)
expect_equal(wmw$lnl, wm0$lnl, tol=1e-3)
expect_equal(coef(wmw), coef(wm0), tol=1e-4)
expect_equal(wmw$vars, wm0$vars, tol=1e-4)
# should be the same if not conditional as well
owls2$w1 <- owls2$w2
suppressWarnings(wmw <- mix(SiblingNegotiation ~ ArrivalTime + (1|Nest), data=owls2, weights=c("w1", "w2"),
nQuad = 1,family=poisson(), verbose=FALSE, cWeights=FALSE))
expect_equal(wmw$lnl, wm0$lnl, tol=1e-3)
expect_equal(coef(wmw), coef(wm0), tol=1e-4)
expect_equal(wmw$vars, wm0$vars, tol=1e-4)
# sub 1 conditional weights work
owls2 <- Owls
owls2$w1 <- 1
owls2$w2 <- 1
owls2$w3 <- 1
owls2$Nest <- as.character(owls2$Nest)
g1 <- subset(owls2, owls2$Nest %in% c("Jeuss", "Lully", "Seiry", "Oleyes"))
# do not renaem subject
owlsrep <- rbind(owls2, g1)
owlsrep$w1 <- ifelse(owlsrep$Nest %in% c("Jeuss", "Lully", "Seiry", "Oleyes"), 0.5, 1)
# two level
suppressWarnings(wm0 <- mix(SiblingNegotiation ~ ArrivalTime + (1|Nest), data=owlsrep, weights=c("w1", "w2"),
family=poisson(), verbose=FALSE,cWeights=TRUE,nQuad=1))
suppressWarnings(wmw <- mix(SiblingNegotiation ~ ArrivalTime + (1|Nest), data=owls2, weights=c("w1", "w2"), nQuad = 1,
family=poisson(), verbose=FALSE))
lm0 <- glmer(SiblingNegotiation ~ ArrivalTime + (1|Nest), data=owls2,
family="poisson")
# check vars
lmevars1 <- data.frame(summary(lm0)$varcor)$sdcor
expect_equal(unname(wm0$vars),
expected = unname(lmevars1)^2,
tolerance = 1e-3)
expect_equal(as.numeric(wm0$lnl), as.numeric(logLik(lm0)), tol=1e-3)
expect_equal(coef(wm0), fixef(lm0), tol=1e-4)
expect_equal(wmw$lnl, wm0$lnl, tol=1e-3)
expect_equal(coef(wmw), coef(wm0), tol=1e-4)
expect_equal(wmw$vars, wm0$vars, tol=1e-4)
# should be the same if not conditional as well
suppressWarnings(wm0 <- mix(SiblingNegotiation ~ ArrivalTime + (1|Nest), data=owlsrep, weights=c("w1", "w2"),
family=poisson(), verbose=FALSE,nQuad=1))
expect_equal(wmw$lnl, wm0$lnl, tol=1e-3)
expect_equal(coef(wmw), coef(wm0), tol=1e-4)
expect_equal(wmw$vars, wm0$vars, tol=1e-4)
})
context("Quadrature works, 3 RE at level 2")
test_that("Quadrature works, 3 RE at level 2", {
skip_on_cran()
skip_if_not_installed("glmmTMB")
require(glmmTMB)
owls2 <- Owls
owls2$w1 <- 1
owls2$w2 <- 1
owls2$w3 <- 1
suppressWarnings(wm0 <- mix(SiblingNegotiation ~ ArrivalTime + (ArrivalTime + NegPerChick|Nest),
data=owls2, weights=c("w1", "w2"),
family=poisson(),verbose=FALSE,nQuad=7))
expect_is(wm0, "WeMixResults")
})
context("Three level model slash and colon")
test_that("Three level model slash and colon", {
sleepstudy2 <- sleepstudy
sleepstudy2$Group <- 1
sleepstudy2$Group <- ifelse(sleepstudy2$Subject %in% c("310", "309", "349", "335"), 2, sleepstudy2$Group)
sleepstudy2$Group <- ifelse(sleepstudy2$Subject %in% c("330", "352", "337", "369"), 3, sleepstudy2$Group)
sleepstudy2$Group <- ifelse(sleepstudy2$Subject %in% c("331", "332", "334", "372"), 4, sleepstudy2$Group)
sleepstudy2$Group <- factor(sleepstudy2$Group)
# make subj that restarts at 1 per group, so there are four Subjects with subj=1, one per group.
sleepstudy2 <- sleepstudy2[order(sleepstudy2$Group, sleepstudy2$Subject),]
sleepstudy2$subj <- factor(rep(c(1:6,1:4,1:4,1:4),each=10))
# table(sleepstudy2$subj, sleepstudy2$Group) # shows each group has a subject 1:4 and group 1 as a 5 and 6 too
ss2 <- sleepstudy2
ss2$w1 <- w1c <- w1 <- rep(1,180)
w2c <- w2 <- rep(1,18)
w3c <- w3 <- rep(1,4)
ss2$w2 <- rep(1, 180)
ss2$w3 <- rep(1, 180)
lmr <- lmer(Reaction ~ Days + (1|Group) + (1|Subject), data=sleepstudy2, REML=FALSE)
# next line is a bad specification: confounds group=1, subj=1 person with group=2, subj=2 person
# lmr <- lmer(Reaction ~ Days + (1|Group) + (1|subj), data=sleepstudy2, REML=FALSE)
wm0 <- mix(Reaction ~ Days + (1|Group:subj) + (1|Group) , data=ss2, weights=c("w1", "w2","w3"))
expect_equal(wm0$lnl, as.numeric(logLik(lmr)), tol=1e-6)
expect_equal(coef(wm0), fixef(lmr), tol=1e-6)
expect_equal(unname(wm0$vars[length(wm0$vars)]), unname(lmr@devcomp$cmp["sigmaML"]^2), tol=1e-4)
#group mean center Days to test group mean centering as well
# we know the mean is 4.5 for all groups because all individuals were tested on all days
sleepstudy2$gmc_days <- sleepstudy2$Days-4.5
lmr <- lmer(Reaction ~ gmc_days + (1|Group/subj), data=sleepstudy2, REML=FALSE)
# above is same as:
# lmr <- lmer(Reaction ~ Days + (1|Group) + (1|Group:subj), data=sleepstudy2, REML=FALSE)
wm0 <- mix(Reaction ~ Days + (1|Group/subj), data=ss2, weights=c("w1", "w2","w3"),center_group = list("Group/subj" = ~Days))
expect_equal(wm0$lnl, as.numeric(logLik(lmr)), tol=1e-6)
expect_equal(unname(coef(wm0)), unname(fixef(lmr)), tol=1e-6)
expect_equal(unname(wm0$vars[length(wm0$vars)]), unname(lmr@devcomp$cmp["sigmaML"]^2), tol=1e-4)
})
# check the format of summary output
context("Summary output format")
test_that("summary output format", {
skip_on_cran()
sleepstudy2 <- sleepstudy
sleepstudy2$Group <- 1
sleepstudy2$Group <- ifelse(sleepstudy2$Subject %in% c("310", "309", "349", "335"), 2, sleepstudy2$Group)
sleepstudy2$Group <- ifelse(sleepstudy2$Subject %in% c("330", "352", "337", "369"), 3, sleepstudy2$Group)
sleepstudy2$Group <- ifelse(sleepstudy2$Subject %in% c("331", "332", "334", "372"), 4, sleepstudy2$Group)
sleepstudy2$Group <- factor(sleepstudy2$Group)
# make subj that restarts at 1 per group, so there are four Subjects with subj=1, one per group.
sleepstudy2 <- sleepstudy2[order(sleepstudy2$Group, sleepstudy2$Subject),]
sleepstudy2$subj <- factor(rep(c(1:6,1:4,1:4,1:4),each=10))
sleepstudy2$carPr <- pnorm(sleepstudy2$Reaction -285,sd=50)
set.seed(2)
sleepstudy2$car <- runif(180) < sleepstudy2$carPr
# table(sleepstudy2$subj, sleepstudy2$Group) # shows each group has a subject 1:4 and group 1 as a 5 and 6 too
ss2 <- sleepstudy2
ss2$w1 <- w1c <- w1 <- rep(1,180)
w2c <- w2 <- rep(1,18)
w3c <- w3 <- rep(1,4)
ss2$w2 <- rep(1, 180)
ss2$w3 <- rep(1, 180)
co0 <- c("Call:",
"mix(formula = Reaction ~ Days + (Days + car | Subject), data = ss2, ",
" weights = c(\"w1\", \"w2\"))", "",
"Variance terms:",
" Level Group Name Variance Std. Error Std.Dev. Corr1 Corr2",
" 2 Subject (Intercept) 629.4 254.1 25.09 ",
" 2 Subject Days 36.7 12.2 6.06 0.17 ",
" 2 Subject carTRUE 648.5 309.3 25.46 -0.42 -0.39",
" 1 Residual 528.8 145.7 23.00 ",
"Groups:",
" Level Group n size mean wgt sum wgt",
" 2 Subject 18 1 18",
" 1 Obs 180 1 180",
"",
"Fixed Effects:",
" Estimate Std. Error t value",
"(Intercept) 252.72 6.53 38.70",
"Days 11.19 1.50 7.44",
"",
"lnl= -868.13 ",
"Intraclass Correlation= 0.713 ")
wm0 <- mix(Reaction ~ Days + (Days+car|Subject), data=ss2, weights=c("w1", "w2"))
withr::with_options(list(digits=2),
co <- capture.output(summary(wm0))
)
expect_equal(co, co0)
# this model is not accuractly estimated so is not a good test of output format.
# nevertheless, it uses data from this block, so appears here
expect_warning(wm0 <- mix(Reaction ~ Days + (car||Subject), data=ss2, weights=c("w1", "w2")), "standard error of variance terms")
varsref <- structure(list(Variance = c(698.125, 213.50 , 1618.92, 808.25 ),
`Std. Error` = c(385.8 , 232 , 850.5 , 162.0972),
Std.Dev. = c(26.4220, 14.6116, 40.2358, 28.42977)),
class = "data.frame",
row.names = c("Subject.(Intercept)", "Subject.carFALSE", "Subject.carTRUE", "Residual"))
coefref <- structure(c(236.842700778804, 9.59410909072588, 6.7061583949811,
1.56079000196884, 35.317194558968, 6.14695704010372), .Dim = 2:3,
.Dimnames = list(c("(Intercept)", "Days"),
c("Estimate", "Std. Error", "t value")))
sum0 <- summary(wm0)
expect_equal(sum0$coef, coefref, tolerance=1e-4)
expect_equal(sum0$vars, varsref, tolerance=4e-3)
})
context("Weighted three level model")
test_that("Weighted three level model", {
sleepstudy2 <- sleepstudy
sleepstudy2$Group <- 1
sleepstudy2$Group <- ifelse(sleepstudy2$Subject %in% c("310", "309", "349", "335"), 2, sleepstudy2$Group)
sleepstudy2$Group <- ifelse(sleepstudy2$Subject %in% c("330", "352", "337", "369"), 3, sleepstudy2$Group)
sleepstudy2$Group <- ifelse(sleepstudy2$Subject %in% c("331", "332", "334", "372"), 4, sleepstudy2$Group)
sleepstudy2$Group <- factor(sleepstudy2$Group)
#sleepstudy2$NewVar <- sleepstudy2$Reaction/100 + rnorm(180,10,5)
sleepstudy2$NewVar <- c(9.43457771946133, 12.7183652886321, 11.242490364576, 8.32049769489214, 21.4769136458094, 12.801691804759, 11.2944480319612, 23.2877967612343, 17.6554807720438, 15.6839576150421,
12.0415732799454, 12.6326929342708, 10.7497191412678, 8.67408880103278, 16.7872963953428, 14.6248001125594, 14.3553882674867, 11.8275049579491, 18.7525279773655, 5.57650843995549,
19.7623973348143, 19.1858339464439, 16.9845555399241, 5.79214586016569, 13.6760775160985, 8.20870365747105, 10.469622698121, 9.5294364080553, 0.474272552642729, 9.90998530696626,
7.85598926755497, 13.6908324775899, 8.363578606421, 12.4585008913942, 20.8371544197279, 6.95572306644593, 4.49104017816579, 13.529348727477, 8.31255622061886, 8.78452997885101,
17.216378291436, 13.6427335980228, 19.498592420836, 7.77906722601152, 15.1266474871889, 12.8030124499139, 6.62054433117513, 8.66974559938289, 11.9540398510289, 14.4271543736958,
10.0496588773239, 15.5177485124749, 9.12750122498998, 10.1801282848309, 14.2562515194416, 19.3249322019648, 10.7021090784807, 19.5336266451282, 7.58133671725254, 19.889529248361,
5.32630146259616, 13.6044477039579, 4.98570236262749, 23.8801945918399, 14.447634009947, 15.6043025989833, 22.411260279495, 16.66695588071, 13.2723030644539, 17.6752877679585,
13.4050771535649, 16.8476633007641, 14.4754165647234, 13.2597805356783, 12.2720489082021, 16.3515751846435, 18.2589186962967, 13.5608781488381, 22.128982372665, 14.0340540905538,
2.08454544532827, 10.8863647522786, 21.0103204678449, 18.4595691535641, 11.9914499297614, 9.96791953740012, 20.0728999616675, 10.7237225063556, 18.6117287544885, 15.6290844263721,
13.584186380461, 12.8188973592142, 8.73588740878364, 22.1049445716224, 8.82190889109436, 16.7501795316447, 11.7894955061621, 8.67751729426691, 8.28560204932531, 25.3725616103908,
14.0620375346182, 14.9882805901214, 17.181272042065, 19.6581190406994, 11.6757977643058, 16.638959090019, 9.68436466184703, 23.1282879805077, 20.7690736926973, 14.3401954001056,
16.2747349705209, 8.61220295948687, 19.869302293393, 5.78478905807136, 6.59988762289215, 20.489844013367, 27.239468499143, 20.5797203990527, 19.4544929129352, 12.5320006698243,
18.2664525619468, 16.6022500953422, 13.0216995414837, 12.4573456587016, 18.6922622366627, 7.37338081469625, 18.2016743123929, 12.9216600044986, 17.2543801572655, 20.1850819959469,
9.46342758905614, 16.1505649461445, 11.644647125078, 17.6992809329738, 13.5726089392145, 11.9006738563255, 9.11529530502925, 11.9473408244809, 13.7037965361481, 21.052192466607,
8.94718178506257, 15.165337112518, 8.23261412105177, 12.2188157585207, 14.1455895553682, 4.94479286361234, 16.9113105275149, 12.7165809453375, 5.58774136619933, 8.31670629685991,
6.67728085523989, 11.7741431392146, 7.91590113917014, 9.13428547474213, 18.5689544076891, 18.1995195536489, 13.3510327364336, 21.4047919488063, 17.7821335119331, 9.19515401988151,
10.7918240635605, 6.63952336871624, 15.3733302621749, 15.5262403960572, 13.4421915705244, 11.5980148647886, 17.3920393426159, 6.6883267354037, 10.6865661488425, 15.1104585129107,
9.99190979270527, 16.54240123493, 13.5687018760124, 13.5592897858269, 19.4019579232211, 15.0561754718535, 9.36036662501835, 4.27173388144395, 19.2040511963462, 6.36268714321505)
ss2 <- sleepstudy2
w1c <- w1 <- rep(1,180)
w2c <- w2 <- rep(1,18)
w3c <- w3 <- rep(1,4)
# unbalanced (non-identical within a Subject), level-1 (obs level) weights
w1c[1:4] <- w1[1:4] <- 2
uR <- sleepstudy2[1:4,]
sleepstudy2 <- rbind(sleepstudy2, uR)
# level-2 weights
w2c[1] <- w2[1] <- 2
w1[ss2$Subject=="308"] <- 2*w1[ss2$Subject=="308"]
sR <- subset(sleepstudy2, Subject == "308")
sR$Subject <- "S308"
sleepstudy2 <- rbind(sleepstudy2, sR)
# level-3 weights
w3c[1] <- w3[1] <- 2
w2[c(1,7,12,13,16,17)] <- 2*w2[c(1,7,12,13,16,17)]
w1[ss2$Group==1] <- 2*w1[ss2$Group==1]
gR <- subset(sleepstudy2, Group %in% 1)
gR$Subject <- paste0("G", gR$Subject)
gR$Group <- "11"
sleepstudy2 <- rbind(sleepstudy2, gR)
ss2$w1 <- w1
ss2$w2 <- rep(w2,each=10)
ss2$w3 <- ifelse(ss2$Subject %in% c("308", "333", "350", "351", "370", "371"),2,1)
# lmr for reference
lmr <- lmer(Reaction ~ Days + (1|Subject) + (1|Group), data=sleepstudy2, REML=FALSE)
#these are the conditional weights, used in the stata comparison
ss2$c1 <- w1c
ss2$c2 <- rep(w2c,each=10)
ss2$c3 <- ifelse(ss2$Subject %in% c("308", "333", "350", "351", "370", "371"),2,1)
#compare against mixed stata resutls
wm0 <- mix(Reaction ~ Days + (1|Subject) + (1|Group), data=ss2, weights=c("w1", "w2","w3"))
expect_equal(wm0$lnl, -1389.0983, tolerance=1e-3)
expect_equal(unname(wm0$coef), c(243.6831,12.67954),tolerance = 1e-3)
expect_equal(unname(wm0$vars), c(360.97 ,756.5857, 1153.704 ),tolerance = 1e-3)
expect_equal(unname( wm0$SE), c( 11.68881 , 2.468474 ),tolerance = 1e-3)
wm1 <- mix(Reaction ~ Days + (1 |Subject) + (1|Group)+ (0+Days|Group), data=ss2, weights=c("w1", "w2","w3"))
expect_equal(wm1$lnl, -1377.6876 , tolerance=1e-3)
expect_equal(unname(wm1$coef), c(247.6234 , 11.71211 ),tolerance = 1e-3)
expect_equal(unname(wm1$vars), c(397.6301 , 220.4182 , 16.86835 , 1022.88 ),tolerance = 1e-3)
expect_equal(unname( wm1$SE), c(9.97378 , 2.60348 ),tolerance = 1e-3)
})
context("Wald Test")
test_that("Test for Wald Tests using ", {
wm1<- mix(Reaction ~ Days + (0+Days|Subject), data=sleepstudyU, weights=c("weight1L1", "weight1L2"))
beta_1 <- waldTest(wm1, type="beta", coefs="Days", hypothesis=10.4073- 1.5458 *1.96)
#wald test and T test are the same for individual betas
expect_equal(beta_1$p,.05,.01)
beta_2<- waldTest(wm1, type="beta", coefs="(Intercept)", hypothesis=251.4051- 6.8246 *1.96)
#wald test and T test are the same for individual betas
expect_equal(beta_2$p,.05,.01)
# Ensure that tests for Lambda run
wm2<- mix(Reaction ~ Days + (1+Days|Subject), data=sleepstudyU, weights=c("weight1L1", "weight1L2"))
waldTest(wm2,type="Lambda")
waldTest(wm2,type="Lambda",coefs = c("Subject.Days.(Intercept)","Subject.Days" ))
})
context("Complex weighted three level model")
test_that("Complex weighted three level model", {
skip_on_cran()
sleepstudy2 <- sleepstudy
sleepstudy2$Group <- 1
sleepstudy2$Group <- ifelse(sleepstudy2$Subject %in% c("310", "309", "349", "335"), 2, sleepstudy2$Group)
sleepstudy2$Group <- ifelse(sleepstudy2$Subject %in% c("330", "352", "337", "369"), 3, sleepstudy2$Group)
sleepstudy2$Group <- ifelse(sleepstudy2$Subject %in% c("331", "332", "334", "372"), 4, sleepstudy2$Group)
sleepstudy2$Group <- factor(sleepstudy2$Group)
sleepstudy2$NewVar <- c(9.43457771946133, 12.7183652886321, 11.242490364576, 8.32049769489214, 21.4769136458094, 12.801691804759, 11.2944480319612, 23.2877967612343, 17.6554807720438, 15.6839576150421,
12.0415732799454, 12.6326929342708, 10.7497191412678, 8.67408880103278, 16.7872963953428, 14.6248001125594, 14.3553882674867, 11.8275049579491, 18.7525279773655, 5.57650843995549,
19.7623973348143, 19.1858339464439, 16.9845555399241, 5.79214586016569, 13.6760775160985, 8.20870365747105, 10.469622698121, 9.5294364080553, 0.474272552642729, 9.90998530696626,
7.85598926755497, 13.6908324775899, 8.363578606421, 12.4585008913942, 20.8371544197279, 6.95572306644593, 4.49104017816579, 13.529348727477, 8.31255622061886, 8.78452997885101,
17.216378291436, 13.6427335980228, 19.498592420836, 7.77906722601152, 15.1266474871889, 12.8030124499139, 6.62054433117513, 8.66974559938289, 11.9540398510289, 14.4271543736958,
10.0496588773239, 15.5177485124749, 9.12750122498998, 10.1801282848309, 14.2562515194416, 19.3249322019648, 10.7021090784807, 19.5336266451282, 7.58133671725254, 19.889529248361,
5.32630146259616, 13.6044477039579, 4.98570236262749, 23.8801945918399, 14.447634009947, 15.6043025989833, 22.411260279495, 16.66695588071, 13.2723030644539, 17.6752877679585,
13.4050771535649, 16.8476633007641, 14.4754165647234, 13.2597805356783, 12.2720489082021, 16.3515751846435, 18.2589186962967, 13.5608781488381, 22.128982372665, 14.0340540905538,
2.08454544532827, 10.8863647522786, 21.0103204678449, 18.4595691535641, 11.9914499297614, 9.96791953740012, 20.0728999616675, 10.7237225063556, 18.6117287544885, 15.6290844263721,
13.584186380461, 12.8188973592142, 8.73588740878364, 22.1049445716224, 8.82190889109436, 16.7501795316447, 11.7894955061621, 8.67751729426691, 8.28560204932531, 25.3725616103908,
14.0620375346182, 14.9882805901214, 17.181272042065, 19.6581190406994, 11.6757977643058, 16.638959090019, 9.68436466184703, 23.1282879805077, 20.7690736926973, 14.3401954001056,
16.2747349705209, 8.61220295948687, 19.869302293393, 5.78478905807136, 6.59988762289215, 20.489844013367, 27.239468499143, 20.5797203990527, 19.4544929129352, 12.5320006698243,
18.2664525619468, 16.6022500953422, 13.0216995414837, 12.4573456587016, 18.6922622366627, 7.37338081469625, 18.2016743123929, 12.9216600044986, 17.2543801572655, 20.1850819959469,
9.46342758905614, 16.1505649461445, 11.644647125078, 17.6992809329738, 13.5726089392145, 11.9006738563255, 9.11529530502925, 11.9473408244809, 13.7037965361481, 21.052192466607,
8.94718178506257, 15.165337112518, 8.23261412105177, 12.2188157585207, 14.1455895553682, 4.94479286361234, 16.9113105275149, 12.7165809453375, 5.58774136619933, 8.31670629685991,
6.67728085523989, 11.7741431392146, 7.91590113917014, 9.13428547474213, 18.5689544076891, 18.1995195536489, 13.3510327364336, 21.4047919488063, 17.7821335119331, 9.19515401988151,
10.7918240635605, 6.63952336871624, 15.3733302621749, 15.5262403960572, 13.4421915705244, 11.5980148647886, 17.3920393426159, 6.6883267354037, 10.6865661488425, 15.1104585129107,
9.99190979270527, 16.54240123493, 13.5687018760124, 13.5592897858269, 19.4019579232211, 15.0561754718535, 9.36036662501835, 4.27173388144395, 19.2040511963462, 6.36268714321505)
ss2 <- sleepstudy2
w1c <- w1 <- rep(1,180)
w2c <- w2 <- rep(1,18)
w3c <- w3 <- rep(1,4)
# unbalanced (non-identical within a Subject), level-1 (obs level) weights
w1c[1:4] <- w1[1:4] <- 2
uR <- sleepstudy2[1:4,]
sleepstudy2 <- rbind(sleepstudy2, uR)
# level-2 weights
w2c[1] <- w2[1] <- 2
w1[ss2$Subject=="308"] <- 2*w1[ss2$Subject=="308"]
sR <- subset(sleepstudy2, Subject == "308")
sR$Subject <- "S308"
sleepstudy2 <- rbind(sleepstudy2, sR)
# level-3 weights
w3c[1] <- w3[1] <- 2
w2[c(1,7,12,13,16,17)] <- 2*w2[c(1,7,12,13,16,17)]
w1[ss2$Group==1] <- 2*w1[ss2$Group==1]
gR <- subset(sleepstudy2, Group %in% 1)
gR$Subject <- paste0("G", gR$Subject)
gR$Group <- "11"
sleepstudy2 <- rbind(sleepstudy2, gR)
ss2$w1 <- w1
ss2$w2 <- rep(w2,each=10)
ss2$w3 <- ifelse(ss2$Subject %in% c("308", "333", "350", "351", "370", "371"),2,1)
ss2$n1 <- ss2$n2 <- ss2$n3 <- 1
#lme reference with duplicated data
lmr <- lmer(Reaction ~ Days + (1|Subject) + (1+Days|Group), data=sleepstudy2, REML=FALSE)
#Does WeMix match to duplicated lmr with one+two random effects
wmr <- mix(Reaction ~ Days + (1|Subject) + (1+Days|Group), data=ss2,weights=c("w1","w2","w3"))
expect_equal(wmr$lnl, as.numeric(logLik(lmr)), tol=1e-3)
expect_equal(coef(wmr), fixef(lmr), tol=1e-4)
lmevars1 <- data.frame(summary(lmr)$varcor)
vars <- lmevars1[is.na(lmevars1$var2),"vcov"]
expect_equal(unname(wmr$vars),
expected = unname(vars),
tolerance = 1e-3)
expect_equal(coef(wmr), fixef(lmr), tol=1e-4)
expect_equal(unname(wmr$vars[length(wmr$vars)]), unname(lmr@devcomp$cmp["sigmaML"]^2), tol=1e-4)
#lme reference with duplicated data
lmr2 <- lmer(Reaction ~ Days + (1 + NewVar |Subject) + (1+Days|Group),
data=sleepstudy2, control=lmerControl(optimizer="bobyqa"), REML=FALSE)
#wemix with weights
suppressWarnings(wmr2 <- mix(Reaction ~ Days + (1+ NewVar |Subject) + (1+Days|Group), data=ss2,weights=c("w1","w2","w3")))
#Does WeMix match to duplicated lmr with two correlated random effects at two levesl
expect_equal(wmr2$lnl, as.numeric(logLik(lmr2)), tol=1e-3)
expect_equal(coef(wmr2), fixef(lmr2), tol=1e-4)
lmevars2 <- data.frame(summary(lmr2)$varcor)
vars <- lmevars2[is.na(lmevars2$var2),"vcov"]
expect_equal(unname(wmr2$vars),
expected = unname(vars),
tolerance = 1e-3)
expect_equal(coef(wmr2), fixef(lmr2), tol=1e-4)
expect_equal(unname(wmr2$vars[length(wmr2$vars)]), unname(lmr2@devcomp$cmp["sigmaML"]^2), tol=1e-4)
})
context("Complex weighted three level model, binomial")
test_that("Complex weighted three level model, binomial", {
skip_on_cran()
sleepstudy2 <- sleepstudy
sleepstudy2$highR <- ifelse(sleepstudy2$Reaction > 340, 1, 0)
sleepstudy2$Group <- 1
sleepstudy2$Group <- ifelse(sleepstudy2$Subject %in% c("310", "309", "349", "335"), 2, sleepstudy2$Group)
sleepstudy2$Group <- ifelse(sleepstudy2$Subject %in% c("330", "352", "337", "369"), 3, sleepstudy2$Group)
sleepstudy2$Group <- ifelse(sleepstudy2$Subject %in% c("331", "332", "334", "372"), 4, sleepstudy2$Group)
sleepstudy2$Group <- factor(sleepstudy2$Group)
sleepstudy2$NewVar <- c(9.43457771946133, 12.7183652886321, 11.242490364576, 8.32049769489214, 21.4769136458094, 12.801691804759, 11.2944480319612, 23.2877967612343, 17.6554807720438, 15.6839576150421,
12.0415732799454, 12.6326929342708, 10.7497191412678, 8.67408880103278, 16.7872963953428, 14.6248001125594, 14.3553882674867, 11.8275049579491, 18.7525279773655, 5.57650843995549,
19.7623973348143, 19.1858339464439, 16.9845555399241, 5.79214586016569, 13.6760775160985, 8.20870365747105, 10.469622698121, 9.5294364080553, 0.474272552642729, 9.90998530696626,
7.85598926755497, 13.6908324775899, 8.363578606421, 12.4585008913942, 20.8371544197279, 6.95572306644593, 4.49104017816579, 13.529348727477, 8.31255622061886, 8.78452997885101,
17.216378291436, 13.6427335980228, 19.498592420836, 7.77906722601152, 15.1266474871889, 12.8030124499139, 6.62054433117513, 8.66974559938289, 11.9540398510289, 14.4271543736958,
10.0496588773239, 15.5177485124749, 9.12750122498998, 10.1801282848309, 14.2562515194416, 19.3249322019648, 10.7021090784807, 19.5336266451282, 7.58133671725254, 19.889529248361,
5.32630146259616, 13.6044477039579, 4.98570236262749, 23.8801945918399, 14.447634009947, 15.6043025989833, 22.411260279495, 16.66695588071, 13.2723030644539, 17.6752877679585,
13.4050771535649, 16.8476633007641, 14.4754165647234, 13.2597805356783, 12.2720489082021, 16.3515751846435, 18.2589186962967, 13.5608781488381, 22.128982372665, 14.0340540905538,
2.08454544532827, 10.8863647522786, 21.0103204678449, 18.4595691535641, 11.9914499297614, 9.96791953740012, 20.0728999616675, 10.7237225063556, 18.6117287544885, 15.6290844263721,
13.584186380461, 12.8188973592142, 8.73588740878364, 22.1049445716224, 8.82190889109436, 16.7501795316447, 11.7894955061621, 8.67751729426691, 8.28560204932531, 25.3725616103908,
14.0620375346182, 14.9882805901214, 17.181272042065, 19.6581190406994, 11.6757977643058, 16.638959090019, 9.68436466184703, 23.1282879805077, 20.7690736926973, 14.3401954001056,
16.2747349705209, 8.61220295948687, 19.869302293393, 5.78478905807136, 6.59988762289215, 20.489844013367, 27.239468499143, 20.5797203990527, 19.4544929129352, 12.5320006698243,
18.2664525619468, 16.6022500953422, 13.0216995414837, 12.4573456587016, 18.6922622366627, 7.37338081469625, 18.2016743123929, 12.9216600044986, 17.2543801572655, 20.1850819959469,
9.46342758905614, 16.1505649461445, 11.644647125078, 17.6992809329738, 13.5726089392145, 11.9006738563255, 9.11529530502925, 11.9473408244809, 13.7037965361481, 21.052192466607,
8.94718178506257, 15.165337112518, 8.23261412105177, 12.2188157585207, 14.1455895553682, 4.94479286361234, 16.9113105275149, 12.7165809453375, 5.58774136619933, 8.31670629685991,
6.67728085523989, 11.7741431392146, 7.91590113917014, 9.13428547474213, 18.5689544076891, 18.1995195536489, 13.3510327364336, 21.4047919488063, 17.7821335119331, 9.19515401988151,
10.7918240635605, 6.63952336871624, 15.3733302621749, 15.5262403960572, 13.4421915705244, 11.5980148647886, 17.3920393426159, 6.6883267354037, 10.6865661488425, 15.1104585129107,
9.99190979270527, 16.54240123493, 13.5687018760124, 13.5592897858269, 19.4019579232211, 15.0561754718535, 9.36036662501835, 4.27173388144395, 19.2040511963462, 6.36268714321505)
ss2 <- sleepstudy2
w1c <- w1 <- rep(1,180)
w2c <- w2 <- rep(1,18)
w3c <- w3 <- rep(1,4)
# unbalanced (non-identical within a Subject), level-1 (obs level) weights
w1c[1:4] <- w1[1:4] <- 2
uR <- sleepstudy2[1:4,]
sleepstudy2 <- rbind(sleepstudy2, uR)
# level-2 weights
w2c[1] <- w2[1] <- 2
w1[ss2$Subject=="308"] <- 2*w1[ss2$Subject=="308"]
sR <- subset(sleepstudy2, Subject == "308")
sR$Subject <- "S308"
sleepstudy2 <- rbind(sleepstudy2, sR)
# level-3 weights
w3c[1] <- w3[1] <- 2
w2[c(1,7,12,13,16,17)] <- 2*w2[c(1,7,12,13,16,17)]
w1[ss2$Group==1] <- 2*w1[ss2$Group==1]
gR <- subset(sleepstudy2, Group %in% 1)
gR$Subject <- paste0("G", gR$Subject)
gR$Group <- "11"
sleepstudy2 <- rbind(sleepstudy2, gR)
ss2$w1 <- w1
ss2$w2 <- rep(w2,each=10)
ss2$w3 <- ifelse(ss2$Subject %in% c("308", "333", "350", "351", "370", "371"),2,1)
ss2$n1 <- ss2$n2 <- ss2$n3 <- 1
#lme reference with duplicated data
glmr <- glmer(highR ~ Days + (1|Subject) + (1+Days|Group), data=sleepstudy2,
family=binomial())
#Does WeMix match to duplicated lmr with one+two random effects
suppressWarnings(wmr <- mix(highR ~ Days + (1|Subject) + (Days|Group), data=ss2,
weights=c("w1","w2","w3"),
family=binomial(),nQuad=1,verbose=FALSE))
expect_equal(as.numeric(wmr$lnl), as.numeric(logLik(glmr)), tol=1e-3)
expect_equal(coef(wmr), fixef(glmr), tol=1e-4)
lmevars1 <- data.frame(summary(glmr)$varcor)
vars <- lmevars1$vcov
expect_equal(unname(wmr$vars),
expected = unname(vars),
tolerance = 1e-3)
#lme reference with duplicated data
glmr2 <- glmer(highR ~ Days + (1 + NewVar|Subject) + (1+Days|Group),
data=sleepstudy2,
family=binomial(),control=glmerControl(optimizer="bobyqa"))
#wemix with weights
suppressWarnings(wmr2 <- mix(highR ~ Days + (1+ NewVar |Subject) + (1+Days|Group), data=ss2,
weights=c("w1","w2","w3"), family=binomial(),nQuad=1,
max_iteration = 20,verbose=FALSE))
#Does WeMix match to duplicated lmr with two correlated random effects at two levesl
expect_equal(as.numeric(wmr2$lnl), as.numeric(logLik(glmr2)), tol=1e-3)
expect_equal(coef(wmr2), fixef(glmr2), tol=1e-4)
lmevars2 <- data.frame(summary(glmr2)$varcor)
vars <- lmevars2$vcov
expect_equal(unname(wmr2$vars),
expected = unname(vars),
tolerance = 1e-3)
})
context("Model Matrix has a hard time with")
test_that("Model Matrix has a hard time with", {
skip_on_cran()
require(EdSurvey)
sdf <- readNAEP(system.file("extdata/data", "M36NT2PM.dat", package = "NAEPprimer"))
gg <- EdSurvey::getData(varnames=c("composite", "dsex", "b017451", "scrpsu", "origwt", "smsrswt"), data=sdf, returnJKreplicates=FALSE)
gg2 <- gg[gg$origwt > 0 & gg$smsrswt > 0,]
suppressMessages(m4 <- mix(mrpcm2 ~ dsex + b017451 + (1|scrpsu), data=gg2, weights=c("origwt", "smsrswt")))
expect_equal(m4$lnl, -81882.3634148408, tol=1e-5)
expect_equal(m4$coef, c(`(Intercept)` = 270.571817625202, dsexFemale = -2.14907551600309, `b017451Once every few weeks` = 3.84433128452533, `b017451About once a week` = 9.19954009166631, `b0174512 or 3 times a week` = 12.8701977366809, `b017451Every day` = 6.29635843233831))
expect_equal(m4$SE, c(`(Intercept)` = 1.14859092700564, dsexFemale = 0.636207386283148, `b017451Once every few weeks` = 1.05190060207538, `b017451About once a week` = 0.986062663493941, `b0174512 or 3 times a week` = 1.00366545196335, `b017451Every day` = 1.10858750669702))
m4varDF <- structure(list(grp = c("scrpsu", "Residual"),
var1 = c("(Intercept)", NA),
var2 = c(NA_character_, NA_character_),
vcov = c(300.695102536901, 969.501918646696),
ngrp = c(672, 16331),
level = c(2, 1),
SEvcov = c(25.0004604746421, 16.2937862202021),
fullGroup = c("scrpsu.(Intercept)", "Residual")),
row.names = c(NA, -2L),
class = "data.frame")
expect_equal(m4$varDF, m4varDF, tol=1e-5)
})
context("examples run")
test_that("examples run", {
skip_on_cran()
ss1 <- sleepstudy
# Create weights
ss1$W1 <- ifelse(ss1$Subject %in% c(308, 309, 310), 2, 1)
ss1$W2 <- 1
# Run random-intercept 2-level model
two_level <- mix(Reaction ~ Days + (1|Subject), data=ss1, weights=c("W1", "W2"))
expect_is(two_level, "WeMixResults")
#Run random-intercept 2-level model with group-mean centering
grp_centered <- mix(Reaction ~ Days + (1|Subject), data=ss1,
weights = c("W1", "W2"),
center_group = list("Subject" = ~Days))
expect_is(grp_centered, "WeMixResults")
#Run three level model with random slope and intercept.
#add group variables for 3 level model
ss1$Group <- 3
ss1$Group <- ifelse(as.numeric(ss1$Subject) %% 10 < 7, 2, ss1$Group)
ss1$Group <- ifelse(as.numeric(ss1$Subject) %% 10 < 4, 1, ss1$Group)
# level-3 weights
ss1$W3 <- ifelse(ss1$Group == 2, 2, 1)
three_level <- mix(Reaction ~ Days + (1|Subject) + (1+Days|Group), data=ss1,
weights=c("W1", "W2", "W3"))
expect_is(three_level, "WeMixResults")
})
Try the WeMix package in your browser
Any scripts or data that you put into this service are public.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.