R/remixTest.R
In ClaireKelley/WeMix: Weighted Mixed-Effects Models using Multilevel Pseudo Maximum Likelihood Estimation
if(FALSE) {
require("HSAUR2")
g2 <- glm(outcome ~ treatment * visit, data=toenail, family=binomial(link="logit"))
y <- 0 + (toenail$outcome=="moderate or severe")
b <- coef(g2)
ll <- binomial(link="logit")
pred0 <- predict(g2, type="response")
precW <- pred0 * (1-pred0)
W <- diag(precW)
X <- model.matrix(outcome ~ treatment *visit, data=toenail)
etaHat <- ll$linkfun(pred0)
z <- etaHat + (y - pred0) * ll$mu.eta(etaHat)
# as in gAnalyticSolve
W12 <- diag(sqrt(precW))
qq <- qr(W12 %*% X)
cbind(solve(t(X) %*% W %*% X) %*% t(X) %*% W %*% z, qr.solve(qq, W12 %*% z), b)
# manually build RE
ll <- binomial(link="logit")
gm2 <- glmer(outcome~treatment*visit+(1|patientID),
data=toenail,
family=binomial,nAGQ=20)
# fixed
Lambda <- getME(gm2, "Lambda")
u <- u0 <- rep(0, 294)
b <- b0 <- rep(0,4)
Z <- model.matrix(~0+patientID, data=toenail)
X <- model.matrix(outcome ~ treatment *visit, data=toenail)
disc <- function(y, mu, W, u) {
d <- y - mu
W <- diag(sqrt(mu*(1-mu)))
t(d) %*% W %*% d + sum(u^2)
}
etaHat <- as.vector(X%*%b + Z%*%Lambda %*%u )
pred0 <- ll$linkinv(etaHat)
precW <- 1/(pred0 * (1-pred0))
W <- diag(as.vector(precW))
drb <- disc(y=y, mu=ll$linkinv(etaHat), W=W,u=u)
# start update
bprev <- b0 + 3
while(sum(abs(bprev - b)) > 0.01) {
etaHat <- as.vector(X%*%b + Z%*%Lambda %*%u )
pred0 <- ll$linkinv(etaHat)
zaug <- c((y - pred0), -1*u)
precW <- 1/(pred0 * (1-pred0))
W12dmde <- diag(as.vector(sqrt(precW)*ll$mu.eta(etaHat)), nrow=length(precW))
U <- W12dmde %*% Z %*% Lambda
V <- W12dmde %*% X
A <- rbind(cbind(U, V ),
cbind(diag(1,ncol(Z)), matrix(0,ncol=ncol(X),nrow=ncol(Z))))
qrA <- qr(A)
dub <- qr.solve(qrA, zaug)
W <- diag(as.vector(precW))
ss <- 1
dnew <- disc(y=y, mu=ll$linkinv(as.vector(X%*%(b+ss*dub[-(1:294)])+ Z %*% Lambda %*%(u+ss*dub[1:294]))), W=W,u=u+ss*dub[1:294])
d0 <- disc(y=y, mu=ll$linkinv(as.vector(X%*%b+ Z %*% Lambda %*%u)), W=W, u=u)
cat("dold=",d0, "\n")
while(dnew > d0) {
cat("ss=",ss, "dnew=",dnew,"\n")
ss <- ss * 0.5
if(abs(ss) < 1e-6) {
ss <- -1
}
dnew <- disc(y=y, mu=ll$linkinv(as.vector(X%*%(b+ss*dub[-(1:294)])+ Z %*% Lambda %*% (u+ss*dub[1:294]))), W=W,u=u+ss*dub[1:294])
}
cat("ss=",ss, "dnew=",dnew,"\n")
bprev <- b
u <- u+ss*dub[1:294]
b <- b+ss*dub[-(1:294)]
print(rbind(cbind(b, rb), c(dnew, drb)))
}
#A <- rbind(cbind( W12 %*% Z %*% Lambda, W12 %*% X), cbind(diag(1, nrow=294), matrix(0, nrow=294, ncol=4)))
rb <- fixef(gm2)
rq <- qr(A)
ufit <- getME(gm2, "u")
uW12 <- ll$linkinv(ufit * (1-ufit))
W12aug <- diag(c(sqrt(precW), rep(1,294)))
Waug <- diag(c(precW, rep(1,294)))
# unit weight
cbind((solve(t(A) %*% Waug %*% A) %*% t(A) %*% Waug %*% zaug)[-(1:294)],
qr.solve(rq, W12aug %*% zaug)[-(1:294)],
qr.solve(rq, zaug)[-(1:294)],
rb)
# plot REs
plot((solve(t(A) %*% Waug %*% A) %*% t(A) %*% Waug %*% zaug)[1:294], getME(gm2, "u"))
plot(qr.solve(rq, W12aug %*% zaug)[1:294], getME(gm2, "u"))
# ufit weights
W12aug <- diag(c(sqrt(precW), ufit))
cbind(qr.solve(rq, W12aug %*% zaug)[-(1:294)], rb)
# plot REs
plot(qr.solve(rq, W12aug %*% zaug)[1:294], getME(gm2, "u"))
# try with no W12agu
cbind(qr.solve(rq, zaug)[-(1:294)], rb)
# plot REs
plot(qr.solve(rq, zaug)[1:294], getME(gm2, "u"))
gm2 <- glmer(outcome~treatment*visit+(1|patientID),
data=toenail,
family=binomial,nAGQ=20)
pred <- predict(gm2, type="response")
precW <- (pred * (1-pred))
X <- model.matrix(outcome ~ treatment *visit, data=toenail)
Zlist <- list(model.matrix(~0+patientID, data=toenail))
weights <- list(rep(1,nrow(X)), rep(1, 294))
groupID <- matrix(as.numeric(toenail$patientID))
lmeVarDF <- data.frame(summary(gm2)$varcor)
lmeVarDF$sdcor <- NULL
lmeVarDF$level <- 2
lmeVarDF$ngrp <- 294
lmeVarDF <- rbind(lmeVarDF, c("Residual", NA, NA, vcov=1, sdcor=1, level=1, ngrp=nrow(X)))
Zlevels <- 2
g0 <- gAnalyticSolve(y, X, Zlist, Zlevels, weights=weights, pWeights=precW, groupID=groupID, lmeVarDF=lmeVarDF, pred0=pred, family=ll)
gg <- g0(v=getME(gm2, "theta"))
cbind(gg,fixef(gm2))
}
if(FALSE){
temp_var_form <- function(lmesummary,data){
ngrp <- lmesummary$ngrps
lmeVarDF <- data.frame(lmesummary$varcor)
# prepare variance data frame
# rename groups in var-cov matrix to be all distinct
ind <- 1
while(sum(grepl(paste0("\\.", ind, "$"), lmeVarDF$grp)) > 0) {
lmeVarDF$grp <- sub(paste0("\\.", ind, "$"), "", lmeVarDF$grp)
ind <- ind + 1
}
lmeVarDF$sdcor <- NULL # remove extraneous column
lmeVarDF$ngrp <- NA
lmeVarDF$grp <- gsub(".", ":", lmeVarDF$grp, fixed=TRUE)
# add column showing how many elements are in each group
for(vari in 1:nrow(lmeVarDF)) {
if(lmeVarDF$grp[vari] == "Residual") {
lmeVarDF$ngrp[vari] <- nrow(data)
} else {
lmeVarDF$ngrp[vari] <- ngrp[names(ngrp) == lmeVarDF$grp[vari]]
}
}
# add column indicating which model level each group belongs to
ngrp2 <- rev(sort(unique(lmeVarDF$ngrp)))
for(grpi in 1:length(ngrp2)) {
lmeVarDF$level[ngrp2[grpi] == lmeVarDF$ngrp] <- grpi + ifelse("Residual" %in% lmeVarDF$grp, 0, 1) # add 1 if no residual
}
return(lmeVarDF)
}
require(lme4)
require(testthat)
###################################################
# TEST 1, theta > 1, unweighted #
#################################
w <- rep(1,180)
w2 <- rep(1,18)
# lmer
lmr <- lmer(Reaction ~ Days + (1|Subject), data=sleepstudy, REML=FALSE)
# make matrixes needed to fit with bw
y <- sleepstudy$Reaction
X <- model.matrix(Reaction ~ Days, data=sleepstudy)
Z1 <- sparse.model.matrix(Reaction ~ 0+Subject, data=sleepstudy)
#this duplicates code from adaptive quad that prepares variance covariance data frame (lines 314-352 and can be removed when integrated)
lmr_sum <- summary(lmr)
lmeVarDF <- temp_var_form(lmr_sum,sleepstudy)
bsq <- bw(y, X, Zlist=list(Z1), Zlevels=c(2), weights=list(w,w2),
groupID=matrix(as.numeric(sleepstudy$Subject),ncol=1),
lmeVarDF = lmeVarDF,
method="qr")
#parsing theta from lmr with names
V0 <- getME(lmr,"theta")
# level-1 replicate
bhatq <- bsq(V0)
# lnl agrees
expect_equal(logLik(lmr)[[1]], bhatq$lnl)
expect_equal(lmr@beta, unname(bhatq$b))
# optimize
bsr <- lmer(Reaction ~ Days + (1|Subject)+(0+Days|Subject), data=sleepstudy, REML=FALSE, devFunOnly=TRUE)
bsqG <- devG(y=y, X=X, Zlist=list(Z1), Zlevels=c(2), weights=list(w,w2), groupID=matrix(as.numeric(sleepstudy$Subject),ncol=1),lmeVarDF = lmeVarDF, method="qr")
#names of theta are important
base_params <- c(.5)
names(base_params) <- names(V0)
system.time(optq <- optim(fn=bsqG, par=base_params, method="Nelder-Mead"))
optq$par
bsq <- bsqG(c(0,0), getBS=TRUE) # recover the bw() object
expect_equal(unname(optq$par), lmr@theta, tolerance =1e-3)
###################################################
# TEST 2: three levels, unweighted #
####################################
w <- rep(1,180)
w2 <- rep(1,18)
# lmer
lmr <- lmer(Reaction ~ Days + (1|Subject)+(0+Days|Subject), data=sleepstudy, REML=FALSE)
# make matrixes needed to fit with bw
y <- sleepstudy$Reaction
X <- model.matrix(Reaction ~ Days, data=sleepstudy)
Z1 <- sparse.model.matrix(Reaction ~ 0+Subject, data=sleepstudy)
Z2 <- sparse.model.matrix(Reaction ~ 0+Days:Subject, data=sleepstudy)
#this duplicates code from adaptive quad that prepares variance covariance data frame (lines 314-352 and can be removed when integrated)
lmr_sum <- summary(lmr)
lmeVarDF <- temp_var_form(lmr_sum,sleepstudy)
bsq <- bw(y, X, Zlist=list(Z1, Z2), Zlevels=c(2,2), weights=list(w,w2),
groupID=matrix(as.numeric(sleepstudy$Subject),ncol=1),
lmeVarDF = lmeVarDF,
method="qr")
#parsing theta from lmr with names
V0 <- getME(lmr,"theta")
# level-1 replicate
bhatq <- bsq(V0)
# lnl agrees
expect_equal(logLik(lmr)[[1]], bhatq$lnl)
expect_equal(lmr@beta, unname(bhatq$b))
# optimize
bsr <- lmer(Reaction ~ Days + (1|Subject)+(0+Days|Subject), data=sleepstudy, REML=FALSE, devFunOnly=TRUE)
bsqG <- devG(y=y, X=X, Zlist=list(Z1, Z2), Zlevels=c(2,2), weights=list(w,w2), groupID=matrix(as.numeric(sleepstudy$Subject),ncol=1),lmeVarDF = lmeVarDF, method="qr")
#names of theta are important
base_params <- c(.5,.5)
names(base_params) <- names(V0)
system.time(optq <- optim(fn=bsqG, par=base_params, method="Nelder-Mead"))
optq$par
bsq <- bsqG(c(0,0), getBS=TRUE) # recover the bw() object
expect_equal(unname(optq$par), lmr@theta, tolerance =1e-3)
###################################################
# TEST 2: two levels, correlations #
####################################
w <- rep(1,180)
w2 <- rep(1,18)
# lmer
lmr <- lmer(Reaction ~ Days + (1+Days|Subject), data=sleepstudy, REML=FALSE)
# make matrixes needed to fit with bw
y <- sleepstudy$Reaction
X <- model.matrix(Reaction ~ Days, data=sleepstudy)
Z1 <- sparse.model.matrix(Reaction ~ 0+Subject, data=sleepstudy)
Z2 <- sparse.model.matrix(Reaction ~ 0+Days:Subject, data=sleepstudy)
Z <- cbind(Z1,Z2)
#this duplicates code from adaptive quad that prepares variance covariance data frame (lines 314-352 and can be removed when integrated)
lmr_sum <- summary(lmr)
lmeVarDF <- temp_var_form(lmr_sum,sleepstudy)
bsq <- bw(y, X, Zlist=list(Z), Zlevels=c(2), weights=list(w,w2),
groupID=matrix(as.numeric(sleepstudy$Subject),ncol=1),
lmeVarDF = lmeVarDF,
method="qr")
#parsing theta from lmr with names
V0 <- getME(lmr,"theta")
# level-1 replicate
bhatq <- bsq(V0)
# lnl agrees
expect_equal(logLik(lmr)[[1]], bhatq$lnl)
expect_equal(lmr@beta, unname(bhatq$b))
# optimize
bsr <- lmer(Reaction ~ Days + (1+Days|Subject), data=sleepstudy, REML=FALSE, devFunOnly=TRUE)
bsqG <- devG(y=y, X=X, Zlist=list(Z1, Z2), Zlevels=c(2,2), weights=list(w,w2), groupID=matrix(as.numeric(sleepstudy$Subject),ncol=1),lmeVarDF = lmeVarDF, method="qr")
#names of theta are important
base_params <- c(.5,.5,.2)
names(base_params) <- names(V0)
system.time(optq <- optim(fn=bsqG, par=base_params, method="Nelder-Mead"))
optq$par
bsq <- bsqG(c(0,0), getBS=TRUE) # recover the bw() object
expect_equal(abs(unname(optq$par)), abs(lmr@theta), tolerance =1e-3)
##################################
# level-1 weights vs replication #
##################################
w <- rep(1,180)
w[1:10] <- 2
sleepstudy2 <- rbind(sleepstudy[1:10,], sleepstudy)
# lmer
lmr <- lmer(Reaction ~ Days + (1|Subject)+(0+Days|Subject), data=sleepstudy2, REML=FALSE)
lmr_sum <- summary(lmr)
lmeVarDF <- temp_var_form(lmr_sum,sleepstudy2)
# bw
bsq <- bw(y, X, Zlist=list(Z1, Z2), Zlevels=c(2,2), weights=list(w,rep(1,18)),
groupID=matrix(as.numeric(sleepstudy$Subject),ncol=1),
lmeVarDF= lmeVarDF,
method="qr")
V0 <- getME(lmr,"theta")
# level-1 replicate
yD <- sleepstudy2$Reaction
XD <- model.matrix(Reaction ~ Days, data=sleepstudy2)
Z1D <- sparse.model.matrix(Reaction ~ 0+Subject, data=sleepstudy2)
Z2D <- sparse.model.matrix(Reaction ~ 0+Days:Subject, data=sleepstudy2)
bsc <- bw(yD, XD, Zlist=list(Z1D, Z2D), Zlevels=c(2,2), weights=list(rep(1,190), rep(1,18)),
groupID=matrix(as.numeric(sleepstudy2$Subject),ncol=1),
lmeVarDF= lmeVarDF,
method="qr")
bhatc <- bsc(V0)
bhatq <- bsq(V0)
expect_equal(bhatc$lnl, bhatq$lnl)
expect_equal(logLik(lmr)[[1]], bhatq$lnl)
expect_equal(logLik(lmr)[[1]], bhatc$lnl) # redundant but helps when only one fails
expect_equal(bhatc$b, bhatq$b)
expect_equal(lmr@beta, unname(bhatq$b))
expect_equal(lmr@beta, unname(bhatc$b)) # redundant but helps when only one fails
##################################
# level-2 weights vs replication #
##################################
w <- rep(1,180)
w[1:10] <- 2
w2 <- rep(1,18)
w2[1] <- 2
sleepstudy2 <- rbind(sleepstudy[1:10,], sleepstudy)
sleepstudy2$Subject <- as.character(sleepstudy2$Subject)
sleepstudy2$Subject[11:20] <- "999"
sleepstudy2$Subject <- factor(sleepstudy2$Subject)
lmr <- lmer(Reaction ~ Days + (1|Subject)+(0+Days|Subject), data=sleepstudy2, REML=FALSE)
lmr_sum <- summary(lmr)
lmeVarDF <- temp_var_form(lmr_sum,sleepstudy2)
#need to get the variance from the non double case becasue it has the right group size
lmr_u <- lmer(Reaction ~ Days + (1|Subject)+(0+Days|Subject), data=sleepstudy, REML=FALSE)
lmr_sum_u <- summary(lmr_u)
lmeVarDF_u <- temp_var_form(lmr_sum_u,sleepstudy2)
bsq <- bw(y, X, Zlist=list(Z1, Z2), Zlevels=c(2,2), weights=list(w,w2),
groupID=matrix(as.numeric(sleepstudy$Subject),ncol=1), method="qr",
lmeVarDF = lmeVarDF_u, REML=FALSE)
bsqG <- devG(y, X, Zlist=list(Z1, Z2), Zlevels=c(2,2), weights=list(w,w2),
groupID=matrix(as.numeric(sleepstudy$Subject),ncol=1),method="qr", lmeVarDF = lmeVarDF_u ,REML=FALSE)
start <- c(0.5,0.5)
names(start) <- names(V0)
system.time(optq <- optim(fn=bsqG,par=start , method="BFGS"))
expect_equal(unname(optq$par), lmr@theta, tolerance=1e-3)
bsq <- bsqG(c(0,0),getBS=TRUE)
bhatq <- bsq(optq$par)
expect_equal(bhatq$lnl, logLik(lmr)[[1]])
V0 <- getME(lmr,"theta")
bhatq <- bsq(V0)
yD <- sleepstudy2$Reaction
XD <- model.matrix(Reaction ~ Days, data=sleepstudy2)
Z1D <- sparse.model.matrix(Reaction ~ 0+Subject, data=sleepstudy2)
Z2D <- sparse.model.matrix(Reaction ~ 0+Days:Subject, data=sleepstudy2)
bscG <- devG(y=yD, X=XD, Zlist=list(Z1D, Z2D), Zlevels=c(2,2), weights=list(rep(1,190), rep(1,19)),
groupID=matrix(as.numeric(sleepstudy2$Subject),ncol=1),
lmeVarDF=lmeVarDF,
method="qr")
start <- c(0.5,0.5)
names(start) <- names(V0)
system.time(optc <- optim(fn=bscG, par=start, method="BFGS"))
bsc <- bscG(c(0,0), getBS=TRUE)
bhatc <- bsc(optc$par)
expect_equal(unname(optc$par), lmr@theta, tolerance=1e-3)
bhatq <- bsq(optq$par)
expect_equal(bhatc$lnl, logLik(lmr)[[1]])
expect_equal(bhatc$lnl, bhatq$lnl)
expect_equal(logLik(lmr)[[1]], bhatq$lnl)
expect_equal(logLik(lmr)[[1]], bhatc$lnl) # redundant but helps when only one fails
expect_equal(bhatc$b, bhatq$b)
expect_equal(lmr@beta, unname(bhatq$b))
expect_equal(lmr@beta, unname(bhatc$b)) # redundant but helps when only one fails
##############################################################################
# level-1 (non-homogenious within group) and level-2 weights vs replication #
#############################################################################
w <- rep(1,180)
w[1:10] <- c(rep(4,4), rep(2,6))
w2 <- rep(1,18)
w2[1] <- 2
sleepstudy2 <- rbind(sleepstudy[1:10,], sleepstudy)
sleepstudy2$Subject <- as.character(sleepstudy2$Subject)
sleepstudy2$Subject[11:20] <- "999"
sleepstudy2 <- sleepstudy2[c(1:4,1:10,11:14,11:20,21:nrow(sleepstudy2)),]
#
sleepstudy2$Subject <- factor(sleepstudy2$Subject)
lmr <- lmer(Reaction ~ Days + (1|Subject)+(0+Days|Subject), data=sleepstudy2, REML=FALSE)
lmr_sum <- summary(lmr)
lmeVarDF <- temp_var_form(lmr_sum,sleepstudy2)
#
bsq <- bw(y, X, Zlist=list(Z1, Z2), Zlevels=c(2,2), weights=list(w, w2),
groupID=matrix(as.numeric(sleepstudy$Subject),ncol=1),
lmeVarDF=lmeVarDF_u,
method="qr", REML=FALSE)
V0 <- getME(lmr,"theta")
bhatq <- bsq(V0)
yD <- sleepstudy2$Reaction
XD <- model.matrix(Reaction ~ Days, data=sleepstudy2)
Z1D <- sparse.model.matrix(Reaction ~ 0+Subject, data=sleepstudy2)
Z2D <- sparse.model.matrix(Reaction ~ 0+Days:Subject, data=sleepstudy2)
bsc <- bw(yD, XD, Zlist=list(Z1D, Z2D), Zlevels=c(2,2), weights=list(rep(1,198), rep(1,19)),
groupID=matrix(as.numeric(sleepstudy2$Subject),ncol=1),
lmeVarDF= lmeVarDF,
method="qr")
bhatc <- bsc(V0)
# tests
expect_equal(bhatc$lnl, bhatq$lnl)
expect_equal(logLik(lmr)[[1]], bhatq$lnl)
expect_equal(logLik(lmr)[[1]], bhatc$lnl) # redundant but helps when only one fails
expect_equal(bhatc$b, bhatq$b)
expect_equal(lmr@beta, unname(bhatq$b))
expect_equal(lmr@beta, unname(bhatc$b)) # redundant but helps when only one fails
#####################################
# un-weighted 3-level model w/colon #
#####################################
sleepstudy2 <- sleepstudy
sleepstudy2$Group <- 1
# group 1 c("331", "333", "350", "351", "370", "371", )
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)
# just subject level variance
lmr <- lmer(Reaction ~ Days + (1|Subject:Group), data=sleepstudy2, REML=FALSE)
#############
# : and / #
#############
sleepstudy2 <- sleepstudy
sleepstudy2$Group <- 1
# group 1 c("331", "333", "350", "351", "370", "371", )
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$Reaction <- sleepstudy2$Reaction + as.numeric(sleepstudy2$Group) * 20
sleepstudy2 <- sleepstudy2[order(sleepstudy2$Group, sleepstudy2$Subject),]
sleepstudy2$subj <- 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
lmr <- lmer(Reaction ~ Days + (1|Group) + (1|Group:subj), 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)
lmr <- lmer(Reaction ~ Days + (1|Group/Subject), data=sleepstudy2, REML=FALSE)
# above is same as:
# lmr <- lmer(Reaction ~ Days + (1|Group) + (1|Group:subj), data=sleepstudy2, REML=FALSE)
#################################
# un-weighted 3-level model #
#################################
sleepstudy2 <- sleepstudy
sleepstudy2$Group <- 1
# group 1 c("331", "333", "350", "351", "370", "371", )
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)
lmr <- lmer(Reaction ~ Days + (1|Subject) + (1|Group), data=sleepstudy2, REML=FALSE)
lmr@theta
lmr_sum <- summary(lmr)
lmeVarDF <- temp_var_form(lmr_sum,sleepstudy2)
yD <- sleepstudy2$Reaction
XD <- model.matrix(Reaction ~ Days, data=sleepstudy2)
Z1D <- sparse.model.matrix(Reaction ~ 0+Subject, data=sleepstudy2)
Z2D <- sparse.model.matrix(Reaction ~ 0+Group, data=sleepstudy2)
grp <- matrix(c(as.numeric(sleepstudy2$Subject), as.numeric(sleepstudy2$Group)),ncol=2)
bsc <- bw(yD, XD, Zlist=list(Z1D, Z2D), Zlevels=c(2,3), weights=list(rep(1,180), rep(1,18), rep(1,4)),
lmeVarDF = lmeVarDF,
groupID=grp, method="qr")
bhatc <- bsc(getME(lmr,"theta"), verbose=FALSE)
expect_equal(bhatc$lnl, as.numeric(logLik(lmr)))
expect_equal(bhatc$b, fixef(lmr))
expect_equal(unname(bhatc$ranef[[2]]), unname(as.matrix(ranef(lmr)[[1]])))
expect_equal(bhatc$sigma, getME(lmr, name="sigma"))
Dc <- devG(yD, XD, Zlist=list(Z1D, Z2D), Zlevels=c(2,3), weights=list(rep(1,180), rep(1,18), rep(1,4)),
lmeVarDF = lmeVarDF,
groupID=grp, method="qr")
start <- c(0.5,0.5)
names(start) <- names(getME(lmr, "theta"))
system.time(opt <- optim(par=start,fn=Dc, method="BFGS"))
expect_equal(opt$par, getME(lmr, "theta"), tolerance=1E-3)
##############################################################################
# weighted 3-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)
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)
# table to show replication suceeded
with(sleepstudy2, table(Subject, Group))
# lmr for reference
lmr <- lmer(Reaction ~ Days + (1|Subject) + (1|Group), data=sleepstudy2, REML=FALSE)
lmr@theta
lmr_sum <- summary(lmr)
lmeVarDF <- temp_var_form(lmr_sum,sleepstudy2)
yD <- sleepstudy2$Reaction
XD <- model.matrix(Reaction ~ Days, data=sleepstudy2)
Z1D <- sparse.model.matrix(Reaction ~ 0+Subject, data=sleepstudy2)
Z2D <- sparse.model.matrix(Reaction ~ 0+Group, data=sleepstudy2)
grp <- matrix(c(as.numeric(sleepstudy2$Subject), as.numeric(sleepstudy2$Group)),ncol=2)
bsc <- bw(yD, XD, Zlist=list(Z1D, Z2D), Zlevels=c(2,3),
weights=list(rep(1,276), rep(1,26), rep(1,5)),
lmeVarDF = lmeVarDF,
groupID=grp, method="qr")
bhatc <- bsc(getME(lmr,"theta"), verbose=FALSE)
expect_equal(bhatc$lnl, as.numeric(logLik(lmr)))
expect_equal(bhatc$b, fixef(lmr))
expect_equal(unname(bhatc$ranef[[2]]), unname(as.matrix(ranef(lmr)[[1]])))
expect_equal(bhatc$sigma, getME(lmr, name="sigma"))
# ss2 is after assignment of group, before other manipulations
lmr0 <- lmer(Reaction ~ Days + (1|Subject) + (1|Group), data=ss2, REML=FALSE)
lmeVarDF0 <- temp_var_form(summary(lmr0),ss2)
y <- ss2$Reaction
X <- model.matrix(Reaction ~ Days, data=ss2)
Z1 <- sparse.model.matrix(Reaction ~ 0+Subject, data=ss2)
Z2 <- sparse.model.matrix(Reaction ~ 0+Group, data=ss2)
grp <- matrix(c(as.numeric(ss2$Subject), as.numeric(ss2$Group)),ncol=2)
bsq <- bw(y, X, Zlist=list(Z1, Z2), Zlevels=c(2,3),
weights=list(w1, w2, w3),
weightsC=list(w1c, w2c, w3c),
lmeVarDF = lmeVarDF0,
groupID=grp, method="qr")
bhatc <- bsc(getME(lmr, "theta"), verbose=FALSE)
bhatq <- bsq(getME(lmr, "theta"), verbose=FALSE)
expect_equal(bhatc$sigma, bhatq$sigma)
expect_equal(bhatq$lnl, as.numeric(logLik(lmr)))
Dc <- devG(y, X, Zlist=list(Z1, Z2), Zlevels=c(2,3),
weights=list(w1, w2, w3),
weightsC=list(w1c, w2c, w3c),
lmeVarDF = lmeVarDF0,
groupID=grp, method="qr")
# test that the optimum is correct
start <- c(0.5,0.5)
names(start) <- names(getME(lmr, "theta"))
system.time(opt <- optim(par=start,fn=Dc, method="BFGS"))
expect_equal(opt$par, getME(lmr, "theta"), tolerance=2E-3)
##add test for more than one random effect in 3 level model
# lmr for reference
lmr <- lmer(Reaction ~ Days + (1|Subject)+ (1 | Group) + (0+Days|Group), data=sleepstudy2, REML=FALSE)
lmr@theta
lmr_sum <- summary(lmr)
lmeVarDF <- temp_var_form(lmr_sum,sleepstudy2)
yD <- sleepstudy2$Reaction
XD <- model.matrix(Reaction ~ Days, data=sleepstudy2)
Z1D <- sparse.model.matrix(Reaction ~ 0+Subject, data=sleepstudy2)
Z22D <- sparse.model.matrix(Reaction ~ 0+Group, data=sleepstudy2)
Z2D <- sparse.model.matrix(Reaction ~ 0+Days:Group, data=sleepstudy2)
grp <- matrix(c(as.numeric(sleepstudy2$Subject), as.numeric(sleepstudy2$Group)),ncol=2)
bsc <- bw(yD, XD, Zlist=list(Z1D, Z22D, Z2D), Zlevels=c(2,3,3),
weights=list(rep(1,276), rep(1,26), rep(1,5)),
lmeVarDF = lmeVarDF,
groupID=grp, method="qr")
bhatc <- bsc(getME(lmr,"theta"), verbose=FALSE)
expect_equal(bhatc$lnl, as.numeric(logLik(lmr)))
expect_equal(bhatc$b, fixef(lmr))
expect_equal(unname(bhatc$ranef[[2]]), unname(as.matrix(ranef(lmr)[[1]])))
expect_equal(bhatc$sigma, getME(lmr, name="sigma"))
Dc <- devG(yD, XD, Zlist=list(Z1D, Z22D, Z2D), Zlevels=c(2,3,3),
weights=list(rep(1,276), rep(1,26), rep(1,5)),
lmeVarDF = lmeVarDF,
groupID=grp, method="qr")
# test that the optimum is correct
start <- c(0.5,0.5,.5)
names(start) <- names(getME(lmr, "theta"))
system.time(opt <- optim(par=start,fn=Dc, method="BFGS"))
expect_equal(opt$par, getME(lmr, "theta"), tolerance=2E-3)
################################################
# weighted 3-level model with 3 refs, covars #
################################################
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)
# table to show replication suceeded
with(sleepstudy2, table(Subject, Group))
# lmr for reference
lmr <- lmer(Reaction ~ Days + (1|Subject) + (1+Days|Group), data=sleepstudy2, REML=FALSE)
lmr@theta
lmr_sum <- summary(lmr)
lmeVarDF <- temp_var_form(lmr_sum,sleepstudy2)
yD <- sleepstudy2$Reaction
XD <- model.matrix(Reaction ~ Days, data=sleepstudy2)
Z1D <- sparse.model.matrix(Reaction ~ 0+Subject, data=sleepstudy2)
Z2D <- sparse.model.matrix(Reaction ~ 0+Group, data=sleepstudy2)
Z22D <- sparse.model.matrix(Reaction ~ 0+Days:Group, data=sleepstudy2)
grp <- matrix(c(as.numeric(sleepstudy2$Subject), as.numeric(sleepstudy2$Group)),ncol=2)
bsc <- bw(yD, XD, Zlist=list(Z1D, cbind(Z2D,Z22D)), Zlevels=c(2,3),
weights=list(rep(1,276), rep(1,26), rep(1,5)),
lmeVarDF = lmeVarDF,
groupID=grp, method="qr")
bhatc <- bsc(getME(lmr,"theta"), verbose=FALSE)
expect_equal(bhatc$lnl, as.numeric(logLik(lmr)))
expect_equal(bhatc$b, fixef(lmr))
expect_equal(unname(bhatc$ranef[[2]]), unname(as.matrix(ranef(lmr)[[1]])))
expect_equal(bhatc$sigma, getME(lmr, name="sigma"))
}
ClaireKelley/WeMix documentation built on May 21, 2019, 6:46 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.
if(FALSE) {
require("HSAUR2")
g2 <- glm(outcome ~ treatment * visit, data=toenail, family=binomial(link="logit"))
y <- 0 + (toenail$outcome=="moderate or severe")
b <- coef(g2)
ll <- binomial(link="logit")
pred0 <- predict(g2, type="response")
precW <- pred0 * (1-pred0)
W <- diag(precW)
X <- model.matrix(outcome ~ treatment *visit, data=toenail)
etaHat <- ll$linkfun(pred0)
z <- etaHat + (y - pred0) * ll$mu.eta(etaHat)
# as in gAnalyticSolve
W12 <- diag(sqrt(precW))
qq <- qr(W12 %*% X)
cbind(solve(t(X) %*% W %*% X) %*% t(X) %*% W %*% z, qr.solve(qq, W12 %*% z), b)
# manually build RE
ll <- binomial(link="logit")
gm2 <- glmer(outcome~treatment*visit+(1|patientID),
data=toenail,
family=binomial,nAGQ=20)
# fixed
Lambda <- getME(gm2, "Lambda")
u <- u0 <- rep(0, 294)
b <- b0 <- rep(0,4)
Z <- model.matrix(~0+patientID, data=toenail)
X <- model.matrix(outcome ~ treatment *visit, data=toenail)
disc <- function(y, mu, W, u) {
d <- y - mu
W <- diag(sqrt(mu*(1-mu)))
t(d) %*% W %*% d + sum(u^2)
}
etaHat <- as.vector(X%*%b + Z%*%Lambda %*%u )
pred0 <- ll$linkinv(etaHat)
precW <- 1/(pred0 * (1-pred0))
W <- diag(as.vector(precW))
drb <- disc(y=y, mu=ll$linkinv(etaHat), W=W,u=u)
# start update
bprev <- b0 + 3
while(sum(abs(bprev - b)) > 0.01) {
etaHat <- as.vector(X%*%b + Z%*%Lambda %*%u )
pred0 <- ll$linkinv(etaHat)
zaug <- c((y - pred0), -1*u)
precW <- 1/(pred0 * (1-pred0))
W12dmde <- diag(as.vector(sqrt(precW)*ll$mu.eta(etaHat)), nrow=length(precW))
U <- W12dmde %*% Z %*% Lambda
V <- W12dmde %*% X
A <- rbind(cbind(U, V ),
cbind(diag(1,ncol(Z)), matrix(0,ncol=ncol(X),nrow=ncol(Z))))
qrA <- qr(A)
dub <- qr.solve(qrA, zaug)
W <- diag(as.vector(precW))
ss <- 1
dnew <- disc(y=y, mu=ll$linkinv(as.vector(X%*%(b+ss*dub[-(1:294)])+ Z %*% Lambda %*%(u+ss*dub[1:294]))), W=W,u=u+ss*dub[1:294])
d0 <- disc(y=y, mu=ll$linkinv(as.vector(X%*%b+ Z %*% Lambda %*%u)), W=W, u=u)
cat("dold=",d0, "\n")
while(dnew > d0) {
cat("ss=",ss, "dnew=",dnew,"\n")
ss <- ss * 0.5
if(abs(ss) < 1e-6) {
ss <- -1
}
dnew <- disc(y=y, mu=ll$linkinv(as.vector(X%*%(b+ss*dub[-(1:294)])+ Z %*% Lambda %*% (u+ss*dub[1:294]))), W=W,u=u+ss*dub[1:294])
}
cat("ss=",ss, "dnew=",dnew,"\n")
bprev <- b
u <- u+ss*dub[1:294]
b <- b+ss*dub[-(1:294)]
print(rbind(cbind(b, rb), c(dnew, drb)))
}
#A <- rbind(cbind( W12 %*% Z %*% Lambda, W12 %*% X), cbind(diag(1, nrow=294), matrix(0, nrow=294, ncol=4)))
rb <- fixef(gm2)
rq <- qr(A)
ufit <- getME(gm2, "u")
uW12 <- ll$linkinv(ufit * (1-ufit))
W12aug <- diag(c(sqrt(precW), rep(1,294)))
Waug <- diag(c(precW, rep(1,294)))
# unit weight
cbind((solve(t(A) %*% Waug %*% A) %*% t(A) %*% Waug %*% zaug)[-(1:294)],
qr.solve(rq, W12aug %*% zaug)[-(1:294)],
qr.solve(rq, zaug)[-(1:294)],
rb)
# plot REs
plot((solve(t(A) %*% Waug %*% A) %*% t(A) %*% Waug %*% zaug)[1:294], getME(gm2, "u"))
plot(qr.solve(rq, W12aug %*% zaug)[1:294], getME(gm2, "u"))
# ufit weights
W12aug <- diag(c(sqrt(precW), ufit))
cbind(qr.solve(rq, W12aug %*% zaug)[-(1:294)], rb)
# plot REs
plot(qr.solve(rq, W12aug %*% zaug)[1:294], getME(gm2, "u"))
# try with no W12agu
cbind(qr.solve(rq, zaug)[-(1:294)], rb)
# plot REs
plot(qr.solve(rq, zaug)[1:294], getME(gm2, "u"))
gm2 <- glmer(outcome~treatment*visit+(1|patientID),
data=toenail,
family=binomial,nAGQ=20)
pred <- predict(gm2, type="response")
precW <- (pred * (1-pred))
X <- model.matrix(outcome ~ treatment *visit, data=toenail)
Zlist <- list(model.matrix(~0+patientID, data=toenail))
weights <- list(rep(1,nrow(X)), rep(1, 294))
groupID <- matrix(as.numeric(toenail$patientID))
lmeVarDF <- data.frame(summary(gm2)$varcor)
lmeVarDF$sdcor <- NULL
lmeVarDF$level <- 2
lmeVarDF$ngrp <- 294
lmeVarDF <- rbind(lmeVarDF, c("Residual", NA, NA, vcov=1, sdcor=1, level=1, ngrp=nrow(X)))
Zlevels <- 2
g0 <- gAnalyticSolve(y, X, Zlist, Zlevels, weights=weights, pWeights=precW, groupID=groupID, lmeVarDF=lmeVarDF, pred0=pred, family=ll)
gg <- g0(v=getME(gm2, "theta"))
cbind(gg,fixef(gm2))
}
if(FALSE){
temp_var_form <- function(lmesummary,data){
ngrp <- lmesummary$ngrps
lmeVarDF <- data.frame(lmesummary$varcor)
# prepare variance data frame
# rename groups in var-cov matrix to be all distinct
ind <- 1
while(sum(grepl(paste0("\\.", ind, "$"), lmeVarDF$grp)) > 0) {
lmeVarDF$grp <- sub(paste0("\\.", ind, "$"), "", lmeVarDF$grp)
ind <- ind + 1
}
lmeVarDF$sdcor <- NULL # remove extraneous column
lmeVarDF$ngrp <- NA
lmeVarDF$grp <- gsub(".", ":", lmeVarDF$grp, fixed=TRUE)
# add column showing how many elements are in each group
for(vari in 1:nrow(lmeVarDF)) {
if(lmeVarDF$grp[vari] == "Residual") {
lmeVarDF$ngrp[vari] <- nrow(data)
} else {
lmeVarDF$ngrp[vari] <- ngrp[names(ngrp) == lmeVarDF$grp[vari]]
}
}
# add column indicating which model level each group belongs to
ngrp2 <- rev(sort(unique(lmeVarDF$ngrp)))
for(grpi in 1:length(ngrp2)) {
lmeVarDF$level[ngrp2[grpi] == lmeVarDF$ngrp] <- grpi + ifelse("Residual" %in% lmeVarDF$grp, 0, 1) # add 1 if no residual
}
return(lmeVarDF)
}
require(lme4)
require(testthat)
###################################################
# TEST 1, theta > 1, unweighted #
#################################
w <- rep(1,180)
w2 <- rep(1,18)
# lmer
lmr <- lmer(Reaction ~ Days + (1|Subject), data=sleepstudy, REML=FALSE)
# make matrixes needed to fit with bw
y <- sleepstudy$Reaction
X <- model.matrix(Reaction ~ Days, data=sleepstudy)
Z1 <- sparse.model.matrix(Reaction ~ 0+Subject, data=sleepstudy)
#this duplicates code from adaptive quad that prepares variance covariance data frame (lines 314-352 and can be removed when integrated)
lmr_sum <- summary(lmr)
lmeVarDF <- temp_var_form(lmr_sum,sleepstudy)
bsq <- bw(y, X, Zlist=list(Z1), Zlevels=c(2), weights=list(w,w2),
groupID=matrix(as.numeric(sleepstudy$Subject),ncol=1),
lmeVarDF = lmeVarDF,
method="qr")
#parsing theta from lmr with names
V0 <- getME(lmr,"theta")
# level-1 replicate
bhatq <- bsq(V0)
# lnl agrees
expect_equal(logLik(lmr)[[1]], bhatq$lnl)
expect_equal(lmr@beta, unname(bhatq$b))
# optimize
bsr <- lmer(Reaction ~ Days + (1|Subject)+(0+Days|Subject), data=sleepstudy, REML=FALSE, devFunOnly=TRUE)
bsqG <- devG(y=y, X=X, Zlist=list(Z1), Zlevels=c(2), weights=list(w,w2), groupID=matrix(as.numeric(sleepstudy$Subject),ncol=1),lmeVarDF = lmeVarDF, method="qr")
#names of theta are important
base_params <- c(.5)
names(base_params) <- names(V0)
system.time(optq <- optim(fn=bsqG, par=base_params, method="Nelder-Mead"))
optq$par
bsq <- bsqG(c(0,0), getBS=TRUE) # recover the bw() object
expect_equal(unname(optq$par), lmr@theta, tolerance =1e-3)
###################################################
# TEST 2: three levels, unweighted #
####################################
w <- rep(1,180)
w2 <- rep(1,18)
# lmer
lmr <- lmer(Reaction ~ Days + (1|Subject)+(0+Days|Subject), data=sleepstudy, REML=FALSE)
# make matrixes needed to fit with bw
y <- sleepstudy$Reaction
X <- model.matrix(Reaction ~ Days, data=sleepstudy)
Z1 <- sparse.model.matrix(Reaction ~ 0+Subject, data=sleepstudy)
Z2 <- sparse.model.matrix(Reaction ~ 0+Days:Subject, data=sleepstudy)
#this duplicates code from adaptive quad that prepares variance covariance data frame (lines 314-352 and can be removed when integrated)
lmr_sum <- summary(lmr)
lmeVarDF <- temp_var_form(lmr_sum,sleepstudy)
bsq <- bw(y, X, Zlist=list(Z1, Z2), Zlevels=c(2,2), weights=list(w,w2),
groupID=matrix(as.numeric(sleepstudy$Subject),ncol=1),
lmeVarDF = lmeVarDF,
method="qr")
#parsing theta from lmr with names
V0 <- getME(lmr,"theta")
# level-1 replicate
bhatq <- bsq(V0)
# lnl agrees
expect_equal(logLik(lmr)[[1]], bhatq$lnl)
expect_equal(lmr@beta, unname(bhatq$b))
# optimize
bsr <- lmer(Reaction ~ Days + (1|Subject)+(0+Days|Subject), data=sleepstudy, REML=FALSE, devFunOnly=TRUE)
bsqG <- devG(y=y, X=X, Zlist=list(Z1, Z2), Zlevels=c(2,2), weights=list(w,w2), groupID=matrix(as.numeric(sleepstudy$Subject),ncol=1),lmeVarDF = lmeVarDF, method="qr")
#names of theta are important
base_params <- c(.5,.5)
names(base_params) <- names(V0)
system.time(optq <- optim(fn=bsqG, par=base_params, method="Nelder-Mead"))
optq$par
bsq <- bsqG(c(0,0), getBS=TRUE) # recover the bw() object
expect_equal(unname(optq$par), lmr@theta, tolerance =1e-3)
###################################################
# TEST 2: two levels, correlations #
####################################
w <- rep(1,180)
w2 <- rep(1,18)
# lmer
lmr <- lmer(Reaction ~ Days + (1+Days|Subject), data=sleepstudy, REML=FALSE)
# make matrixes needed to fit with bw
y <- sleepstudy$Reaction
X <- model.matrix(Reaction ~ Days, data=sleepstudy)
Z1 <- sparse.model.matrix(Reaction ~ 0+Subject, data=sleepstudy)
Z2 <- sparse.model.matrix(Reaction ~ 0+Days:Subject, data=sleepstudy)
Z <- cbind(Z1,Z2)
#this duplicates code from adaptive quad that prepares variance covariance data frame (lines 314-352 and can be removed when integrated)
lmr_sum <- summary(lmr)
lmeVarDF <- temp_var_form(lmr_sum,sleepstudy)
bsq <- bw(y, X, Zlist=list(Z), Zlevels=c(2), weights=list(w,w2),
groupID=matrix(as.numeric(sleepstudy$Subject),ncol=1),
lmeVarDF = lmeVarDF,
method="qr")
#parsing theta from lmr with names
V0 <- getME(lmr,"theta")
# level-1 replicate
bhatq <- bsq(V0)
# lnl agrees
expect_equal(logLik(lmr)[[1]], bhatq$lnl)
expect_equal(lmr@beta, unname(bhatq$b))
# optimize
bsr <- lmer(Reaction ~ Days + (1+Days|Subject), data=sleepstudy, REML=FALSE, devFunOnly=TRUE)
bsqG <- devG(y=y, X=X, Zlist=list(Z1, Z2), Zlevels=c(2,2), weights=list(w,w2), groupID=matrix(as.numeric(sleepstudy$Subject),ncol=1),lmeVarDF = lmeVarDF, method="qr")
#names of theta are important
base_params <- c(.5,.5,.2)
names(base_params) <- names(V0)
system.time(optq <- optim(fn=bsqG, par=base_params, method="Nelder-Mead"))
optq$par
bsq <- bsqG(c(0,0), getBS=TRUE) # recover the bw() object
expect_equal(abs(unname(optq$par)), abs(lmr@theta), tolerance =1e-3)
##################################
# level-1 weights vs replication #
##################################
w <- rep(1,180)
w[1:10] <- 2
sleepstudy2 <- rbind(sleepstudy[1:10,], sleepstudy)
# lmer
lmr <- lmer(Reaction ~ Days + (1|Subject)+(0+Days|Subject), data=sleepstudy2, REML=FALSE)
lmr_sum <- summary(lmr)
lmeVarDF <- temp_var_form(lmr_sum,sleepstudy2)
# bw
bsq <- bw(y, X, Zlist=list(Z1, Z2), Zlevels=c(2,2), weights=list(w,rep(1,18)),
groupID=matrix(as.numeric(sleepstudy$Subject),ncol=1),
lmeVarDF= lmeVarDF,
method="qr")
V0 <- getME(lmr,"theta")
# level-1 replicate
yD <- sleepstudy2$Reaction
XD <- model.matrix(Reaction ~ Days, data=sleepstudy2)
Z1D <- sparse.model.matrix(Reaction ~ 0+Subject, data=sleepstudy2)
Z2D <- sparse.model.matrix(Reaction ~ 0+Days:Subject, data=sleepstudy2)
bsc <- bw(yD, XD, Zlist=list(Z1D, Z2D), Zlevels=c(2,2), weights=list(rep(1,190), rep(1,18)),
groupID=matrix(as.numeric(sleepstudy2$Subject),ncol=1),
lmeVarDF= lmeVarDF,
method="qr")
bhatc <- bsc(V0)
bhatq <- bsq(V0)
expect_equal(bhatc$lnl, bhatq$lnl)
expect_equal(logLik(lmr)[[1]], bhatq$lnl)
expect_equal(logLik(lmr)[[1]], bhatc$lnl) # redundant but helps when only one fails
expect_equal(bhatc$b, bhatq$b)
expect_equal(lmr@beta, unname(bhatq$b))
expect_equal(lmr@beta, unname(bhatc$b)) # redundant but helps when only one fails
##################################
# level-2 weights vs replication #
##################################
w <- rep(1,180)
w[1:10] <- 2
w2 <- rep(1,18)
w2[1] <- 2
sleepstudy2 <- rbind(sleepstudy[1:10,], sleepstudy)
sleepstudy2$Subject <- as.character(sleepstudy2$Subject)
sleepstudy2$Subject[11:20] <- "999"
sleepstudy2$Subject <- factor(sleepstudy2$Subject)
lmr <- lmer(Reaction ~ Days + (1|Subject)+(0+Days|Subject), data=sleepstudy2, REML=FALSE)
lmr_sum <- summary(lmr)
lmeVarDF <- temp_var_form(lmr_sum,sleepstudy2)
#need to get the variance from the non double case becasue it has the right group size
lmr_u <- lmer(Reaction ~ Days + (1|Subject)+(0+Days|Subject), data=sleepstudy, REML=FALSE)
lmr_sum_u <- summary(lmr_u)
lmeVarDF_u <- temp_var_form(lmr_sum_u,sleepstudy2)
bsq <- bw(y, X, Zlist=list(Z1, Z2), Zlevels=c(2,2), weights=list(w,w2),
groupID=matrix(as.numeric(sleepstudy$Subject),ncol=1), method="qr",
lmeVarDF = lmeVarDF_u, REML=FALSE)
bsqG <- devG(y, X, Zlist=list(Z1, Z2), Zlevels=c(2,2), weights=list(w,w2),
groupID=matrix(as.numeric(sleepstudy$Subject),ncol=1),method="qr", lmeVarDF = lmeVarDF_u ,REML=FALSE)
start <- c(0.5,0.5)
names(start) <- names(V0)
system.time(optq <- optim(fn=bsqG,par=start , method="BFGS"))
expect_equal(unname(optq$par), lmr@theta, tolerance=1e-3)
bsq <- bsqG(c(0,0),getBS=TRUE)
bhatq <- bsq(optq$par)
expect_equal(bhatq$lnl, logLik(lmr)[[1]])
V0 <- getME(lmr,"theta")
bhatq <- bsq(V0)
yD <- sleepstudy2$Reaction
XD <- model.matrix(Reaction ~ Days, data=sleepstudy2)
Z1D <- sparse.model.matrix(Reaction ~ 0+Subject, data=sleepstudy2)
Z2D <- sparse.model.matrix(Reaction ~ 0+Days:Subject, data=sleepstudy2)
bscG <- devG(y=yD, X=XD, Zlist=list(Z1D, Z2D), Zlevels=c(2,2), weights=list(rep(1,190), rep(1,19)),
groupID=matrix(as.numeric(sleepstudy2$Subject),ncol=1),
lmeVarDF=lmeVarDF,
method="qr")
start <- c(0.5,0.5)
names(start) <- names(V0)
system.time(optc <- optim(fn=bscG, par=start, method="BFGS"))
bsc <- bscG(c(0,0), getBS=TRUE)
bhatc <- bsc(optc$par)
expect_equal(unname(optc$par), lmr@theta, tolerance=1e-3)
bhatq <- bsq(optq$par)
expect_equal(bhatc$lnl, logLik(lmr)[[1]])
expect_equal(bhatc$lnl, bhatq$lnl)
expect_equal(logLik(lmr)[[1]], bhatq$lnl)
expect_equal(logLik(lmr)[[1]], bhatc$lnl) # redundant but helps when only one fails
expect_equal(bhatc$b, bhatq$b)
expect_equal(lmr@beta, unname(bhatq$b))
expect_equal(lmr@beta, unname(bhatc$b)) # redundant but helps when only one fails
##############################################################################
# level-1 (non-homogenious within group) and level-2 weights vs replication #
#############################################################################
w <- rep(1,180)
w[1:10] <- c(rep(4,4), rep(2,6))
w2 <- rep(1,18)
w2[1] <- 2
sleepstudy2 <- rbind(sleepstudy[1:10,], sleepstudy)
sleepstudy2$Subject <- as.character(sleepstudy2$Subject)
sleepstudy2$Subject[11:20] <- "999"
sleepstudy2 <- sleepstudy2[c(1:4,1:10,11:14,11:20,21:nrow(sleepstudy2)),]
#
sleepstudy2$Subject <- factor(sleepstudy2$Subject)
lmr <- lmer(Reaction ~ Days + (1|Subject)+(0+Days|Subject), data=sleepstudy2, REML=FALSE)
lmr_sum <- summary(lmr)
lmeVarDF <- temp_var_form(lmr_sum,sleepstudy2)
#
bsq <- bw(y, X, Zlist=list(Z1, Z2), Zlevels=c(2,2), weights=list(w, w2),
groupID=matrix(as.numeric(sleepstudy$Subject),ncol=1),
lmeVarDF=lmeVarDF_u,
method="qr", REML=FALSE)
V0 <- getME(lmr,"theta")
bhatq <- bsq(V0)
yD <- sleepstudy2$Reaction
XD <- model.matrix(Reaction ~ Days, data=sleepstudy2)
Z1D <- sparse.model.matrix(Reaction ~ 0+Subject, data=sleepstudy2)
Z2D <- sparse.model.matrix(Reaction ~ 0+Days:Subject, data=sleepstudy2)
bsc <- bw(yD, XD, Zlist=list(Z1D, Z2D), Zlevels=c(2,2), weights=list(rep(1,198), rep(1,19)),
groupID=matrix(as.numeric(sleepstudy2$Subject),ncol=1),
lmeVarDF= lmeVarDF,
method="qr")
bhatc <- bsc(V0)
# tests
expect_equal(bhatc$lnl, bhatq$lnl)
expect_equal(logLik(lmr)[[1]], bhatq$lnl)
expect_equal(logLik(lmr)[[1]], bhatc$lnl) # redundant but helps when only one fails
expect_equal(bhatc$b, bhatq$b)
expect_equal(lmr@beta, unname(bhatq$b))
expect_equal(lmr@beta, unname(bhatc$b)) # redundant but helps when only one fails
#####################################
# un-weighted 3-level model w/colon #
#####################################
sleepstudy2 <- sleepstudy
sleepstudy2$Group <- 1
# group 1 c("331", "333", "350", "351", "370", "371", )
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)
# just subject level variance
lmr <- lmer(Reaction ~ Days + (1|Subject:Group), data=sleepstudy2, REML=FALSE)
#############
# : and / #
#############
sleepstudy2 <- sleepstudy
sleepstudy2$Group <- 1
# group 1 c("331", "333", "350", "351", "370", "371", )
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$Reaction <- sleepstudy2$Reaction + as.numeric(sleepstudy2$Group) * 20
sleepstudy2 <- sleepstudy2[order(sleepstudy2$Group, sleepstudy2$Subject),]
sleepstudy2$subj <- 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
lmr <- lmer(Reaction ~ Days + (1|Group) + (1|Group:subj), 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)
lmr <- lmer(Reaction ~ Days + (1|Group/Subject), data=sleepstudy2, REML=FALSE)
# above is same as:
# lmr <- lmer(Reaction ~ Days + (1|Group) + (1|Group:subj), data=sleepstudy2, REML=FALSE)
#################################
# un-weighted 3-level model #
#################################
sleepstudy2 <- sleepstudy
sleepstudy2$Group <- 1
# group 1 c("331", "333", "350", "351", "370", "371", )
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)
lmr <- lmer(Reaction ~ Days + (1|Subject) + (1|Group), data=sleepstudy2, REML=FALSE)
lmr@theta
lmr_sum <- summary(lmr)
lmeVarDF <- temp_var_form(lmr_sum,sleepstudy2)
yD <- sleepstudy2$Reaction
XD <- model.matrix(Reaction ~ Days, data=sleepstudy2)
Z1D <- sparse.model.matrix(Reaction ~ 0+Subject, data=sleepstudy2)
Z2D <- sparse.model.matrix(Reaction ~ 0+Group, data=sleepstudy2)
grp <- matrix(c(as.numeric(sleepstudy2$Subject), as.numeric(sleepstudy2$Group)),ncol=2)
bsc <- bw(yD, XD, Zlist=list(Z1D, Z2D), Zlevels=c(2,3), weights=list(rep(1,180), rep(1,18), rep(1,4)),
lmeVarDF = lmeVarDF,
groupID=grp, method="qr")
bhatc <- bsc(getME(lmr,"theta"), verbose=FALSE)
expect_equal(bhatc$lnl, as.numeric(logLik(lmr)))
expect_equal(bhatc$b, fixef(lmr))
expect_equal(unname(bhatc$ranef[[2]]), unname(as.matrix(ranef(lmr)[[1]])))
expect_equal(bhatc$sigma, getME(lmr, name="sigma"))
Dc <- devG(yD, XD, Zlist=list(Z1D, Z2D), Zlevels=c(2,3), weights=list(rep(1,180), rep(1,18), rep(1,4)),
lmeVarDF = lmeVarDF,
groupID=grp, method="qr")
start <- c(0.5,0.5)
names(start) <- names(getME(lmr, "theta"))
system.time(opt <- optim(par=start,fn=Dc, method="BFGS"))
expect_equal(opt$par, getME(lmr, "theta"), tolerance=1E-3)
##############################################################################
# weighted 3-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)
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)
# table to show replication suceeded
with(sleepstudy2, table(Subject, Group))
# lmr for reference
lmr <- lmer(Reaction ~ Days + (1|Subject) + (1|Group), data=sleepstudy2, REML=FALSE)
lmr@theta
lmr_sum <- summary(lmr)
lmeVarDF <- temp_var_form(lmr_sum,sleepstudy2)
yD <- sleepstudy2$Reaction
XD <- model.matrix(Reaction ~ Days, data=sleepstudy2)
Z1D <- sparse.model.matrix(Reaction ~ 0+Subject, data=sleepstudy2)
Z2D <- sparse.model.matrix(Reaction ~ 0+Group, data=sleepstudy2)
grp <- matrix(c(as.numeric(sleepstudy2$Subject), as.numeric(sleepstudy2$Group)),ncol=2)
bsc <- bw(yD, XD, Zlist=list(Z1D, Z2D), Zlevels=c(2,3),
weights=list(rep(1,276), rep(1,26), rep(1,5)),
lmeVarDF = lmeVarDF,
groupID=grp, method="qr")
bhatc <- bsc(getME(lmr,"theta"), verbose=FALSE)
expect_equal(bhatc$lnl, as.numeric(logLik(lmr)))
expect_equal(bhatc$b, fixef(lmr))
expect_equal(unname(bhatc$ranef[[2]]), unname(as.matrix(ranef(lmr)[[1]])))
expect_equal(bhatc$sigma, getME(lmr, name="sigma"))
# ss2 is after assignment of group, before other manipulations
lmr0 <- lmer(Reaction ~ Days + (1|Subject) + (1|Group), data=ss2, REML=FALSE)
lmeVarDF0 <- temp_var_form(summary(lmr0),ss2)
y <- ss2$Reaction
X <- model.matrix(Reaction ~ Days, data=ss2)
Z1 <- sparse.model.matrix(Reaction ~ 0+Subject, data=ss2)
Z2 <- sparse.model.matrix(Reaction ~ 0+Group, data=ss2)
grp <- matrix(c(as.numeric(ss2$Subject), as.numeric(ss2$Group)),ncol=2)
bsq <- bw(y, X, Zlist=list(Z1, Z2), Zlevels=c(2,3),
weights=list(w1, w2, w3),
weightsC=list(w1c, w2c, w3c),
lmeVarDF = lmeVarDF0,
groupID=grp, method="qr")
bhatc <- bsc(getME(lmr, "theta"), verbose=FALSE)
bhatq <- bsq(getME(lmr, "theta"), verbose=FALSE)
expect_equal(bhatc$sigma, bhatq$sigma)
expect_equal(bhatq$lnl, as.numeric(logLik(lmr)))
Dc <- devG(y, X, Zlist=list(Z1, Z2), Zlevels=c(2,3),
weights=list(w1, w2, w3),
weightsC=list(w1c, w2c, w3c),
lmeVarDF = lmeVarDF0,
groupID=grp, method="qr")
# test that the optimum is correct
start <- c(0.5,0.5)
names(start) <- names(getME(lmr, "theta"))
system.time(opt <- optim(par=start,fn=Dc, method="BFGS"))
expect_equal(opt$par, getME(lmr, "theta"), tolerance=2E-3)
##add test for more than one random effect in 3 level model
# lmr for reference
lmr <- lmer(Reaction ~ Days + (1|Subject)+ (1 | Group) + (0+Days|Group), data=sleepstudy2, REML=FALSE)
lmr@theta
lmr_sum <- summary(lmr)
lmeVarDF <- temp_var_form(lmr_sum,sleepstudy2)
yD <- sleepstudy2$Reaction
XD <- model.matrix(Reaction ~ Days, data=sleepstudy2)
Z1D <- sparse.model.matrix(Reaction ~ 0+Subject, data=sleepstudy2)
Z22D <- sparse.model.matrix(Reaction ~ 0+Group, data=sleepstudy2)
Z2D <- sparse.model.matrix(Reaction ~ 0+Days:Group, data=sleepstudy2)
grp <- matrix(c(as.numeric(sleepstudy2$Subject), as.numeric(sleepstudy2$Group)),ncol=2)
bsc <- bw(yD, XD, Zlist=list(Z1D, Z22D, Z2D), Zlevels=c(2,3,3),
weights=list(rep(1,276), rep(1,26), rep(1,5)),
lmeVarDF = lmeVarDF,
groupID=grp, method="qr")
bhatc <- bsc(getME(lmr,"theta"), verbose=FALSE)
expect_equal(bhatc$lnl, as.numeric(logLik(lmr)))
expect_equal(bhatc$b, fixef(lmr))
expect_equal(unname(bhatc$ranef[[2]]), unname(as.matrix(ranef(lmr)[[1]])))
expect_equal(bhatc$sigma, getME(lmr, name="sigma"))
Dc <- devG(yD, XD, Zlist=list(Z1D, Z22D, Z2D), Zlevels=c(2,3,3),
weights=list(rep(1,276), rep(1,26), rep(1,5)),
lmeVarDF = lmeVarDF,
groupID=grp, method="qr")
# test that the optimum is correct
start <- c(0.5,0.5,.5)
names(start) <- names(getME(lmr, "theta"))
system.time(opt <- optim(par=start,fn=Dc, method="BFGS"))
expect_equal(opt$par, getME(lmr, "theta"), tolerance=2E-3)
################################################
# weighted 3-level model with 3 refs, covars #
################################################
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)
# table to show replication suceeded
with(sleepstudy2, table(Subject, Group))
# lmr for reference
lmr <- lmer(Reaction ~ Days + (1|Subject) + (1+Days|Group), data=sleepstudy2, REML=FALSE)
lmr@theta
lmr_sum <- summary(lmr)
lmeVarDF <- temp_var_form(lmr_sum,sleepstudy2)
yD <- sleepstudy2$Reaction
XD <- model.matrix(Reaction ~ Days, data=sleepstudy2)
Z1D <- sparse.model.matrix(Reaction ~ 0+Subject, data=sleepstudy2)
Z2D <- sparse.model.matrix(Reaction ~ 0+Group, data=sleepstudy2)
Z22D <- sparse.model.matrix(Reaction ~ 0+Days:Group, data=sleepstudy2)
grp <- matrix(c(as.numeric(sleepstudy2$Subject), as.numeric(sleepstudy2$Group)),ncol=2)
bsc <- bw(yD, XD, Zlist=list(Z1D, cbind(Z2D,Z22D)), Zlevels=c(2,3),
weights=list(rep(1,276), rep(1,26), rep(1,5)),
lmeVarDF = lmeVarDF,
groupID=grp, method="qr")
bhatc <- bsc(getME(lmr,"theta"), verbose=FALSE)
expect_equal(bhatc$lnl, as.numeric(logLik(lmr)))
expect_equal(bhatc$b, fixef(lmr))
expect_equal(unname(bhatc$ranef[[2]]), unname(as.matrix(ranef(lmr)[[1]])))
expect_equal(bhatc$sigma, getME(lmr, name="sigma"))
}
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.