%\VignetteIndexEntry{LMMstar: overview}
%\VignetteEngine{R.rsp::asis}
%\VignetteKeyword{PDF}
%\VignetteKeyword{vignette}
%\VignetteKeyword{package}

## chunk 2
library(LMMstar)

## chunk 3
utils::packageVersion("LMMstar")

## chunk 4
library(ggplot2) ## autoplot method
library(nlme) ## ranef method
library(lava) ## iid, information, manifest methods

## * Illustrative dataset

## chunk 5
data(gastricbypassL, package = "LMMstar")
head(gastricbypassL)

## chunk 6
gastricbypassL$group <- as.factor(as.numeric(gastricbypassL$id)%%2)

## chunk 7
gastricbypassL$baseline <- gastricbypassL$time<0

## chunk 8
data(gastricbypassW, package = "LMMstar")
head(gastricbypassW)

## chunk 9
gastricbypassW$group <- as.numeric(gastricbypassW$id)%%2

## chunk 10
gastricbypassL.NNA <- gastricbypassL[!is.na(gastricbypassL$glucagonAUC),]

## * Visualization & descriptive statistics
## ** Graphical display

## chunk 11
scatterplot(gastricbypassW, ## left panel
            columns = c("weight1","weight2","weight3","weight4")) 

## chunk 12
scatterplot(weight~time|id, data = gastricbypassL, ## right panel
            type.diag = "hist", group = "group")

## chunk 14
gg.spa <- ggplot(gastricbypassL, aes(x=time,y=weight,group=id,color=id))
gg.spa <- gg.spa + geom_point() + geom_line()
gg.spa

## ** Missing data patterns

## chunk 15
mp <- summarizeNA(gastricbypassL)
mp

## chunk 16
plot(mp)

## ** Summary statistics

## chunk 18
sss <- summarize(weight+glucagonAUC ~ time, data = gastricbypassL, na.rm = TRUE)
print(sss, digits = 3)

## chunk 19
sss2 <- summarize(weight ~ time|id, data = gastricbypassL, na.rm = TRUE)
print(sss2, digits = 3)

## chunk 20
plot(sss2, type = "mean") ## left panel
plot(sss2, type = "sd") ## middle panel
plot(sss2, type = "cor") ## right panel

## ** Correlation and partial correlations

## chunk 22
partialCor(weight + glucagonAUC ~ 1, by = "group", data = gastricbypassL)

## chunk 23
gastricbypassL.0 <- gastricbypassL[gastricbypassL$group==0,]
rho <- cor.test(gastricbypassL.0$weight, gastricbypassL.0$glucagonAUC)
c(rho$estimate, p.value = rho$p.value)

## chunk 24
partialCor(weight + glucagonAUC ~ 1, by = "group", effects = "Dunnett",
           data = gastricbypassL)

## chunk 25
partialCor(weight4 + glucagonAUC4 ~ weight1,
           data = gastricbypassW)

## chunk 26
partialCor(list(weight1 ~ glucagonAUC1, weight4 ~ glucagonAUC4),
           data = gastricbypassW)

## * Multiple Student's t-tests

## chunk 27
restt <- t.test(weight1 ~ group, data = gastricbypassW)
c(estimate = unname(diff(restt$estimate)), p.value = restt$p.value)

## chunk 28
## single step max-test adjustment (see help(confint.Wald_lmm) for details)
mt.test(weight1+weight2+weight3+weight4~group, data = gastricbypassW)

## chunk 29
## no adjustment
mt.test(weight1+weight2+weight3+weight4~group, data = gastricbypassW, method = "none")

## chunk 30
## bonferroni adjustment
mt.test(weight1+weight2+weight3+weight4~group, data = gastricbypassW, method = "bonferroni")

## * Linear mixed model (LMM)
## ** Classical covariance patterns

## chunk 31
eId.lmm <- lmm(glucagonAUC ~ visit*group, repetition = ~time|id, 
               structure = "ID", data = gastricbypassL)
eId.lmm
cat(" modeled residual variance-covariance: \n");sigma(eId.lmm)

## chunk 32
eInd.lmm <- lmm(glucagonAUC ~ visit*group, repetition = ~time|id, 
                structure = "IND", data = gastricbypassL)
eInd.lmm
cat(" modeled residual variance-covariance: \n");sigma(eInd.lmm)

## chunk 33
eCS.lmm <- lmm(glucagonAUC ~ visit*group, repetition = ~time|id,
               structure = "CS", data = gastricbypassL)
eCS.lmm
cat(" modeled residual variance-covariance: \n");sigma(eCS.lmm)

## chunk 34
eTOE.lmm <- lmm(glucagonAUC ~ visit*group, repetition = ~time|id,
                structure = "TOEPLITZ", data = gastricbypassL)
eTOE.lmm
cat(" modeled residual correlation: \n");cov2cor(sigma(eTOE.lmm))

## chunk 35
eUN.lmm <- lmm(glucagonAUC ~ visit*group, repetition = ~time|id,
               structure = "UN", data = gastricbypassL)
eUN.lmm
cat(" modeled residual variance-covariance: \n");sigma(eUN.lmm)

## chunk 36
eSCS.lmm <- lmm(glucagonAUC ~ visit*group, repetition = ~time|id,
                structure = CS(group~1), data = gastricbypassL)
eSCS.lmm

## chunk 37
eSUN.lmm <- lmm(glucagonAUC ~ visit*group, repetition = ~time|id,
                structure = UN(group~1), data = gastricbypassL)
eSUN.lmm

## chunk 38
sigma(eSCS.lmm)

## chunk 39
sigma(eSUN.lmm)

## chunk 40
eBCS.lmm <- lmm(glucagonAUC ~ visit*group, repetition = ~time|id,
                structure = CS(~baseline, type = "homogeneous"), data = gastricbypassL)
eBCS.lmm
cat(" modeled residual variance-covariance: \n");sigma(eBCS.lmm)

## chunk 42
eBUN.lmm <- lmm(glucagonAUC ~ visit*group, repetition = ~time|id,
                structure = CS(~baseline, type = "heterogeneous"), data = gastricbypassL)
eBUN.lmm
cat(" modeled residual variance-covariance: \n");sigma(eBUN.lmm)

## ** User-specific covariance patterns

## chunk 43
rho.2block <- function(p,n.time,X){
  rho <- matrix(1, nrow = n.time, ncol = n.time)
  rho[1,2] <- rho[2,1] <- rho[4,5] <- rho[5,4] <- p["rho1"]
  rho[1,3] <- rho[3,1] <- rho[4,6] <- rho[6,4] <- p["rho2"]
  rho[2,3] <- rho[3,2] <- rho[5,6] <- rho[6,5] <- p["rho3"]
  rho[4:6,1:3] <- rho[1:3,4:6] <- p["rho4"]
  return(rho)
}
Rho <- rho.2block(p = c(rho1=0.25,rho2=0.5,rho3=0.4,rho4=0.1),
                  n.time = 6)
Rho

## chunk 44
set.seed(11)
Y <- mvtnorm::rmvnorm(1000, mean = rep(0,6), sigma = Rho)
dfW <- cbind(id = 1:NROW(Y), as.data.frame(Y))
dfL <- reshape2::melt(dfW, id.vars = "id", variable.name = "time")
dfL[dfL$id %in% 1:2,]

## chunk 45
dfL[dfL$id==1,]

## chunk 46
dfL[dfL$id==2,]

## chunk 47
myStruct <- CUSTOM(~time,
                   FCT.sigma = function(p,n.time,X){rep(p,n.time)}, ## function f
                   init.sigma = c("sigma"=1),
                   FCT.rho = rho.2block, ## function g
                   init.rho = c("rho1"=0.25,"rho2"=0.25,"rho3"=0.25,"rho4"=0.25))

## chunk 48
e.lmmCUSTOM <- lmm(value~time, repetition=~time|id,
                   structure = myStruct, data=dfL,
                   df = FALSE) ## df = FALSE to save computation time
logLik(e.lmmCUSTOM)

## chunk 49
cov2cor(sigma(e.lmmCUSTOM))

## chunk 50
system.time(
  e.lmmDEFAULT.CS <- lmm(value~time, repetition = ~time|id,
                         structure = "CS", data = dfL,
                         df = FALSE)
)

## chunk 51
myCS <- CUSTOM(~1,
               FCT.sigma = function(p,n.time,X){rep(p,n.time)},
               init.sigma = c("sigma"=1), 
               FCT.rho = function(p,n.time,X){p+diag(1-p,n.time,n.time)},
               init.rho = c("rho"=0.5))

## chunk 52
system.time(
  e.lmmCUSTOM.CS <- lmm(value~time, repetition = ~time|id,
                        structure = myCS, data = dfL,
                        df = FALSE)
)

## chunk 53
logLik(e.lmmDEFAULT.CS)
logLik(e.lmmCUSTOM.CS)

## chunk 54
e.lmmDEFAULT.CS$opt$n.iter
e.lmmCUSTOM.CS$opt$n.iter

## chunk 55
myCS.wD <- CUSTOM(~1,
                  FCT.sigma = function(p,n.time,X){rep(p,n.time)},
                  dFCT.sigma = function(p,n.time,X){list(sigma = rep(1,n.time))},
                  d2FCT.sigma = function(p,n.time,X){list(sigma = rep(0,n.time))},
                  init.sigma = c("sigma"=1),
                  FCT.rho = function(p,n.time,X){p+diag(1-p,n.time,n.time)},
                  dFCT.rho = function(p,n.time,X){list(rho = 1-diag(1,n.time,n.time))},
                  d2FCT.rho = function(p,n.time,X){list(rho = matrix(0,n.time,n.time))},
                  init.rho = c("rho"=0.5))

## chunk 56
system.time(
  e.lmmCUSTOMwD.CS <- lmm(value~time,
                          repetition = ~time|id,
                          structure = myCS.wD, 
                          data = dfL, df = FALSE
                          )
)

## ** Estimation procedure

## chunk 57
eCS.lmm.bis <- update(eCS.lmm, control = list(trace = 2))

## chunk 58
init.all <- coef(eCS.lmm, effects = "all")
eCS.lmm.bis <- update(eCS.lmm, control = list(init = init.all, trace = 1))

## chunk 59
init.mean <- coef(eCS.lmm, effects = "mean")
eCS.lmm.bis <- update(eCS.lmm, control = list(init = init.mean, trace = 2))

## chunk 60
init.vcov <- sigma(eCS.lmm)
eCS.lmm.bis <- update(eCS.lmm, control = list(init = init.vcov, trace = 1))

## ** Model output

## chunk 61
summary(eUN.lmm)

## chunk 62
summary(eUN.lmm, hide.mean = TRUE)

## chunk 63
oo <- capture.output(summary(eUN.lmm, hide.fit = TRUE, hide.data = TRUE, hide.cor = TRUE, hide.var = TRUE, hide.sd = TRUE))
cat(sapply(oo[-(1:2)],paste0,"\n"))

## ** Extract estimated coefficients

## chunk 64
coef(eUN.lmm)

## chunk 65
coef(eUN.lmm, effects = "variance")

## chunk 66
coef(eUN.lmm, effects = "variance", transform.k = "sd")

## ** Extract estimated coefficient and associated uncertainty

## chunk 67
model.tables(eUN.lmm)

## chunk 68
model.tables(eUN.lmm, effect = c("variance","correlation"))

## chunk 69
model.tables(eUN.lmm, columns = c("estimate","p.value"))

## chunk 70
model.tables(eUN.lmm, columns = add("statistic"))

## ** Extract estimated residual variance-covariance structure

## chunk 71
Sigma <- sigma(eUN.lmm)
Sigma

## chunk 72
round(Sigma,3)

## chunk 73
cov2cor(Sigma)

## chunk 74
round(cov2cor(Sigma), 3)

## chunk 75
sigma(eUN.lmm, cluster = 5)

## chunk 76
newdata <- data.frame(id = "X", time = c("-13","-1","1","13"))
sigma(eUN.lmm, cluster = newdata)

## ** Marginal effects

## chunk 77
effects(eUN.lmm, variable = "group")

## chunk 78
effects(eUN.lmm, type = "change", variable = "group")

## chunk 79
effects(eUN.lmm, type = "auc", variable = "group")

## chunk 80
effects(eUN.lmm, type = "difference", variable = "group")

## chunk 81
effects(eUN.lmm, type = c("change","difference"), variable = "group")

## chunk 82
effects(eUN.lmm, type = c("auc","difference"), variable = "group")

## ** Random effects

## chunk 83
eRI.lmm <- lmm(glucagonAUC ~ visit*group + (1|id), data = gastricbypassL)
eRI.lmm

## chunk 84
eNRI.lmm <- lmm(glucagonAUC ~ visit*group + (1|id/baseline), data = gastricbypassL)
eNRI.lmm

## chunk 85
head(ranef(eRI.lmm, format = "wide"))

## chunk 86
head(ranef(eNRI.lmm, format = "wide"))

## chunk 88
ranef(eRI.lmm, effects = "variance",
      format = "wide")

## chunk 89
ranef(eRI.lmm, effects = "variance", format = "wide")
cat(" \n")
cat(" \n")

## chunk 90
ranef(eNRI.lmm, effects = "variance",
      format = "wide")

## chunk 91
head(ranef(eRI.lmm, se = TRUE))

## ** Sum of squares

## chunk 92
sigma2 <- coef(eCS.lmm, effect = "variance")^2
tau <- coef(eCS.lmm, effect = "correlation")*sigma2
delta <- unname(sigma2 - tau)

## chunk 93
df.res <- df.residual(eCS.lmm)
SSE <- df.res * delta
c(df.res = df.res, SSE = SSE)

## chunk 94
eBeta.lmm <- coef(eCS.lmm)
eVcov.lmm <- vcov(eCS.lmm, type.information = "expected")

## chunk 95
attr(model.matrix(eCS.lmm),"assign")

## chunk 96
SSRstar.time <- eBeta.lmm[2:4] %*% solve(eVcov.lmm[2:4,2:4]) %*% eBeta.lmm[2:4] 
SSRstar.group <- eBeta.lmm[5] %*% solve(eVcov.lmm[5,5]) %*% eBeta.lmm[5] 

## chunk 97
SSR.time <- as.double(SSRstar.time * delta)
SSR.group <- as.double(SSRstar.group * delta)
c(time = SSR.time, group = SSR.group)

## ** Proportion of explained variance and partial correlation

## chunk 99
c(SSR.time/ (SSR.time + SSE),
  SSR.group/ (SSR.group + SSE))

## chunk 100
eCS.R2 <- partialCor(eCS.lmm, R2 = TRUE)
summary(eCS.R2)

## chunk 101
aCS.aov <- anova(eCS.lmm)$multivariate
setNames(with(aCS.aov, statistic*df.num/(statistic*df.num+df.denom)), aCS.aov$test)

## ** Model diagnostic

## chunk 103
eUN.diagW <- residuals(eUN.lmm, type = "normalized", format = "wide")
colnames(eUN.diagW) <- gsub("normalized.","",colnames(eUN.diagW))
head(eUN.diagW)

## chunk 104
eUN.diagL <- residuals(eUN.lmm, type = "normalized", format = "long", keep.data = TRUE)
head(eUN.diagL)

## chunk 105
plot(eUN.lmm, type = "scatterplot")

## chunk 106
plot(eUN.lmm, type = "scatterplot2")

## chunk 107
plot(eUN.lmm, type = "correlation", type.residual = "response")
plot(eUN.lmm, type = "correlation", type.residual = "normalized")

## chunk 109
plot(eUN.lmm, type = "qqplot", engine = "qqtest",
     facet = ~time, labeller = "label_both", facet_nrow=1)
## Note: the qqtest package to be installed to use the argument engine.plot = "qqtest" 

## chunk 110
eUN.lmm_profile <- profile(eUN.lmm, effects = c("sigma","rho(-13,-1)"))
plot(eUN.lmm_profile)

## ** Visualize model fit

## chunk 111
## left panel
plot(eUN.lmm, type = "fit", color = "group", size.text = 20)

## chunk 112
## middle panel
plot(eUN.lmm, type = "fit", color = "group", ci.alpha = NA, size.text = 20)

## chunk 113
## right panel
plot(eUN.lmm, type = "fit", obs.alpha = 0.25, ci = FALSE, size.text = 20)

## chunk 114
## add baseline weight
gastricbypassLB <- merge(gastricbypassL, gastricbypassW[c("id","weight1")], by = "id")

eUN.lmmB <- lmm(glucagonAUC ~ weight1 + visit*group, repetition = ~time|id,
                structure = "UN", data = gastricbypassLB)

## chunk 115
## left panel
plot(eUN.lmmB, type = "fit", color = "group", ci = FALSE, size.text = 20)

## chunk 116
## middel panel
plot(eUN.lmmB, type = "fit", color = "group", ci = FALSE, size.text = 20,
     at = data.frame(weight1 = 150), obs.alpha = 0.2)

## chunk 117
## right panel
gg.traj <- autoplot(eUN.lmmB, type = "fit", color = "group", size.text = 20, facet =~id)
gg.traj$plot + theme(legend.position = "bottom")

## ** Partial residuals

## chunk 118
df.pres <- residuals(eUN.lmmB, type = "partial", variable = "weight1", keep.data = TRUE)
head(df.pres)

## chunk 119
## left panel
plot(eUN.lmmB, type = "partial", variable = "weight1")
## right panel
plot(eUN.lmmB, type = "partial", variable = c("(Intercept)","weight1"))

## ** Statistical inference (single model, linear)

## chunk 121
anova(eUN.lmm)

## chunk 122
anova(eUN.lmm, effects = c("visit3-visit2=0"))

## chunk 123
e.anova <- anova(eUN.lmm, effects = c("visit3-visit2=0","visit4-visit2=0"))
summary(e.anova)

## chunk 124
library(multcomp)
summary(anova(eUN.lmm, effects = mcp(visit = "Tukey")))

## chunk 125
try(
  anova(eUN.lmm,
        effects = c("log(k).-1=0","log(k).1=0","log(k).13=0"))
)

## chunk 126
name.coef <- rownames(confint(eUN.lmm, effects = "all"))
name.varcoef <- grep("^k",name.coef, value = TRUE)
C <- matrix(0, nrow = 3, ncol = length(name.coef), dimnames = list(name.varcoef, name.coef))
diag(C[name.varcoef,name.varcoef]) <- 1
C[,1:9]

## chunk 127
anova(eUN.lmm, effects = C)

## ** Statistical inference (multiple models, linear)

## chunk 128
Manova <- rbind(anova(eInd.lmm, effects = "visit3:group1 = 0", robust = FALSE),
                anova(eCS.lmm, effects = "visit3:group1 = 0", robust = FALSE),
                anova(eUN.lmm, effects = "visit3:group1 = 0", robust = FALSE),
                name = c("Ind","CS","UN"))
summary(Manova) 

## ** Statistical inference (single model, non-linear)

## chunk 129
e.ANCOVA <- lm(weight4 ~ weight1 + group, data = gastricbypassW)
summary(e.ANCOVA)$coef

## chunk 130
gastricbypassL14 <- gastricbypassL[gastricbypassL$visit %in% c(1,4),]
gastricbypassL14$visit <- droplevels(gastricbypassL14$visit)
e.lmmANCOVA <- lmm(weight ~ visit + visit:group, repetition = ~visit|id,
                   data = gastricbypassL14)

## chunk 131
lava::estimate(e.lmmANCOVA, f = function(p){
  c(Y1 = as.double(p["rho(1,4)"]*p["k.4"]),
    X1 = as.double(p["visit4:group1"]-p["rho(1,4)"]*p["k.4"]*p["visit1:group1"]))
})

## ** Baseline adjustment

## chunk 132
gastricbypassL$treat <- baselineAdjustment(gastricbypassL, variable = "group",
                                repetition = ~visit|id, constrain = c("1","2"),
                                new.level = "none")
table(treat = gastricbypassL$treat,
      visit = gastricbypassL$visit,
      group = gastricbypassL$group)

## chunk 133
table(treat = gastricbypassL$treat, visit = gastricbypassL$visit, group = gastricbypassL$group)[,,1,drop=FALSE]

## chunk 134
table(treat = gastricbypassL$treat, visit = gastricbypassL$visit, group = gastricbypassL$group)[,,2,drop=FALSE]

## chunk 135
gastricbypassL$treat2 <- baselineAdjustment(gastricbypassL, variable = "group",
                                 repetition = ~visit|id, constrain = c("1","2"))
table(treat = gastricbypassL$treat2,
      visit = gastricbypassL$visit,
      group = gastricbypassL$group)

## chunk 136
table(treat = gastricbypassL$treat2, visit = gastricbypassL$visit, group = gastricbypassL$group)[,,1,drop=FALSE]

## chunk 137
table(treat = gastricbypassL$treat2, visit = gastricbypassL$visit, group = gastricbypassL$group)[,,2,drop=FALSE]

## chunk 138
gastricbypassL$visitXtreat <- baselineAdjustment(gastricbypassL, variable = "group",
                                      repetition = ~visit|id, constrain = c("1","2"),
                                      collapse.time = ".")

table(treat = gastricbypassL$visitXtreat,
      visit = gastricbypassL$visit,
      group = gastricbypassL$group)

## chunk 139
table(treat = gastricbypassL$visitXtreat, visit = gastricbypassL$visit, group = gastricbypassL$group)[,,1,drop=FALSE]

## chunk 140
table(treat = gastricbypassL$visitXtreat, visit = gastricbypassL$visit, group = gastricbypassL$group)[,,2,drop=FALSE]

## chunk 141
eC.lmm <- lmm(glucagonAUC ~ visitXtreat, data = gastricbypassL,
              repetition = ~visit|id, structure = "UN")
coef(eC.lmm) ## change from baseline

## chunk 142
eC2.lmm <- lmm(glucagonAUC ~ 0 + visitXtreat, data = gastricbypassL,
              repetition = ~visit|id, structure = "UN")
coef(eC2.lmm) ## absolute value

## chunk 143
colnames(model.matrix(glucagonAUC ~ treat*visit, data = gastricbypassL))

## chunk 144
eC3.lmm <- lmm(glucagonAUC ~ treat2*visit, data = gastricbypassL,
               repetition = ~visit|id, structure = "UN")

## chunk 145
model.tables(eC3.lmm)

## chunk 146
plot(eC3.lmm, color = "group", ci = FALSE, size.text = 20, obs.alpha = 0.1)

## chunk 147
effects(eC3.lmm, variable = "treat2", type = "difference")

## ** Predictions

## chunk 148
news <- gastricbypassL[gastricbypassL$id==2,]
news$glucagon <- 0
predict(eUN.lmm, newdata = news, se = TRUE)

## chunk 149
X.12 <- model.matrix(formula(eUN.lmm), news)
X.12

## chunk 150
X.12 %*% coef(eUN.lmm)

## chunk 151
newd <- rbind(
  data.frame(id = 1, time = -13, visit = "1", group = 0, glucagonAUC = coef(eUN.lmm)["(Intercept)"]),
  data.frame(id = 1, time = 1, visit = "3", group = 0, glucagonAUC = NA),
  data.frame(id = 2, time = -13, visit = "1", group = 0, glucagonAUC = 50),
  data.frame(id = 2, time = 1, visit = "3", group = 0, glucagonAUC = NA)
)
predict(eUN.lmm, newdata = newd, type = "dynamic", keep.data = TRUE)

## chunk 152
mu1 <- coef(eUN.lmm)["(Intercept)"]
mu3 <- mu1 + coef(eUN.lmm)["visit3"]
Omega_11 <- sigma(eUN.lmm)[1,1]
Omega_31 <- sigma(eUN.lmm)[3,1]
as.double(mu3 + Omega_31 * (50 - mu1) / Omega_11)

## * Equivalence with other statistical methods
## ** Welch two sample t-test

## chunk 153
t.test(weight4 ~ group, data = gastricbypassW)

## chunk 154
e.ttest4 <- lmm(weight4 ~ group, structure = IND(~group), 
               data = gastricbypassW, trace = FALSE)
model.tables(e.ttest4)

## ** Paired t-test

## chunk 155
t.test(gastricbypassW$weight4, gastricbypassW$weight1, paired = TRUE)

## chunk 156
e.lmm2tt <- lmm(weight ~ visit, repetition = ~visit|id, structure = "UN",
                data = gastricbypassL)
model.tables(e.lmm2tt)["visit4",,drop=FALSE]

## ** Welch two sample t-test on the change

## chunk 157
gastricbypassW.0 <- gastricbypassW[gastricbypassW$group==0,]
gastricbypassW.1 <- gastricbypassW[gastricbypassW$group==1,]
t.test(gastricbypassW.0$weight4-gastricbypassW.0$weight1,
       gastricbypassW.1$weight4-gastricbypassW.1$weight1)

## chunk 158
e.lmm2tt2 <- lmm(weight ~ visit*group, repetition = ~visit|id, structure = UN(~group),
                 data = gastricbypassL)
model.tables(e.lmm2tt2)["visit4:group1",,drop=FALSE]

## ** Multiple Student's t-test

## chunk 159
e.ttest1 <- lmm(weight1 ~ group, structure = IND(~group), 
                data = gastricbypassW, trace = FALSE)
e.ttest2 <- lmm(weight2 ~ group, structure = IND(~group), 
                data = gastricbypassW, trace = FALSE)
e.ttest3 <- lmm(weight3 ~ group, structure = IND(~group), 
                data = gastricbypassW, trace = FALSE)

## chunk 160
e.mttest <- rbind(anova(e.ttest1, effects = "group=0"),
                  anova(e.ttest2, effects = "group=0"),
                  anova(e.ttest3, effects = "group=0"),
                  anova(e.ttest4, effects = "group=0"))
model.tables(e.mttest, method = "bonferroni")

## chunk 161
e.mttest2 <- mlmm(weight ~ group, structure = IND(~group),
                  data = gastricbypassL, trace = FALSE,
                  effects = "group1=0", by = "time", repetition = ~time|id)
model.tables(e.mttest2, method = "single-step2")

## chunk 162
mt.test(weight1+weight2+weight3+weight4~group, data = gastricbypassW)

## ** Linear regression on the change 

## chunk 163
gastricbypassW$changeG41 <- gastricbypassW$glucagonAUC4-gastricbypassW$glucagonAUC1
e.change41 <- lm(changeG41 ~ weight1, data = gastricbypassW)
summary(e.change41)$coef

## chunk 164
gastricbypassL41 <- gastricbypassL[gastricbypassL$visit %in% c(1,4),]
gastricbypassL41$visit <- droplevels(gastricbypassL41$visit)
gastricbypassL41$weight1 <- gastricbypassW$weight1[gastricbypassL41$id]

e.lmm41 <- lmm(glucagonAUC ~ visit + visit*weight1,
               repetition =~ visit|id, structure = "UN",
               data = gastricbypassL41)
model.tables(e.lmm41)

## chunk 165
index.missing41 <- which(is.na(gastricbypassW$changeG41))
index.missing41

## ** Correlation between changes 

## chunk 166
gastricbypassW$changeG41 <- gastricbypassW$glucagonAUC4-gastricbypassW$glucagonAUC1
gastricbypassW$changeW41 <- gastricbypassW$weight4-gastricbypassW$weight1

## chunk 167
cor.test(gastricbypassW$changeW41, gastricbypassW$changeG41)

## chunk 168
e2.change41 <- lm(changeG41 ~ changeW41, data = gastricbypassW)
summary(e2.change41)$coef

## chunk 169
keep.col <- c("id","weight1","weight4","glucagonAUC1","glucagonAUC4")
gastricbypassL4 <- reshape(gastricbypassW[,keep.col], direction = "long",
                           idvar = "id", varying = 2:5, timevar = "type", v.names = "value")
gastricbypassL4$type <- factor(gastricbypassL4$type, labels = keep.col[-1])
gastricbypassL4 <- gastricbypassL4[order(gastricbypassL4$id),]
head(gastricbypassL4)

## chunk 170
e.lmm4 <- lmm(value ~ type,
              repetition = ~type|id, structure = "UN",
              data = gastricbypassL4)

## chunk 171
sigma.lmm4 <- sigma(e.lmm4)
sigma.lmm4

## chunk 172
Mcon <- cbind(c(-1,1,0,0),c(0,0,-1,1))
sigmeChange.lmm4 <- t(Mcon) %*% sigma.lmm4 %*% Mcon
dimnames(sigmeChange.lmm4) <- list(c("d.weight","d.glucagonAUC"),
                                   c("d.weight","d.glucagonAUC"))
sigmeChange.lmm4

## chunk 173
cov2cor(sigmeChange.lmm4)[1,2]
sigmeChange.lmm4[1,2]/sigmeChange.lmm4[1,1]

## chunk 174
estimate(e.lmm4, function(p){
  Sigma.change <- t(Mcon) %*% sigma(e.lmm4, p = p) %*% Mcon
  c(cor = cov2cor(Sigma.change)[1,2],
    beta = Sigma.change[1,2]/Sigma.change[1,1])
})

## * Missing values and imputation

## chunk 175
sss <- summarize(glucagonAUC ~ time, data = gastricbypassL, na.rm = TRUE)
cbind(sss[,1:4], pc = paste0(100 * sss$missing / (sss$missing + sss$observed), "%"))

## chunk 176
summarizeNA(data = gastricbypassL, repetition = ~ time|id)

## chunk 177
## long format
gastricbypassL32 <- gastricbypassL[gastricbypassL$visit %in% c(3,2),]
gastricbypassL32$visit <- droplevels(gastricbypassL32$visit)
gastricbypassL32$weight1 <- gastricbypassW$weight1[gastricbypassL32$id]
## wide format
gastricbypassW$changeG32 <- gastricbypassW$glucagonAUC3-gastricbypassW$glucagonAUC2

## ** Full information approach

## chunk 178
e.lmm32 <- lmm(glucagonAUC ~ visit + visit*weight1,
               repetition =~ visit|id, structure = "UN",
               data = gastricbypassL32)
model.tables(e.lmm32)

## chunk 179
e.change32 <- lm(changeG32 ~ weight1, data = gastricbypassW)
summary(e.change32)$coef

## chunk 180
coef(lm(changeG32 ~ weight1, data = gastricbypassW[-c(5,15),]))

## chunk 181
gastricbypassWA <- fitted(e.lmm32, type = "outcome", format = "wide")
gastricbypassWA$change32 <- gastricbypassWA$glucagonAUC_3 - gastricbypassWA$glucagonAUC_2
gastricbypassWA$weight1 <- gastricbypassW$weight1[match(gastricbypassW$id,gastricbypassWA$id)]
coef(lm(change32 ~ weight1, data = gastricbypassWA))

## ** Complete case approach

## chunk 182
e.lmmCC <- lmmCC(e.change32, repetition = changeG32 ~ glucagonAUC3-glucagonAUC2|id)
model.tables(e.lmmCC)

## chunk 183
summary(e.change32)$coef

## chunk 184
gastricbypassW$changeW32 <- gastricbypassW$weight3 - gastricbypassW$weight2

e2g.change32 <- lm(changeG32 ~ changeW32 + group, data = gastricbypassW)
summary(e2g.change32)$coef

## chunk 185
e2.lmmCC <-  lmmCC(e2g.change32, repetition = list(changeG32 ~ glucagonAUC3-glucagonAUC2|id,
                                                   changeW32 ~ weight3-weight2|id))
model.tables(e2.lmmCC)

## ** Imputation 

## chunk 186
eUN.lmmNA <- lmm(glucagonAUC ~ time, repetition = ~time|id, data = gastricbypassL)
nobs(eUN.lmmNA)

## chunk 187
eData <- fitted(eUN.lmmNA, type = "outcome", keep.data = TRUE)
eData <- eData[order(eData$id,eData$time),]
eData[eData$id %in% eData[eData$impute,"id"],c("id","visit","time","glucagonAUC","impute")]

## chunk 188
ggplot(eData, aes(x=time,y=glucagonAUC, group=id)) + geom_line() + geom_point(aes(color=impute))

## chunk 190
index.na <- which(is.na(gastricbypassL$glucagonAUC))
set.seed(1)
fitted(eUN.lmmNA, type = "impute", se = c(TRUE,TRUE))[index.na]
set.seed(2)
fitted(eUN.lmmNA, type = "impute", se = c(TRUE,TRUE))[index.na]
set.seed(3)
fitted(eUN.lmmNA, type = "impute", se = c(TRUE,TRUE))[index.na]

## ** Multiple imputation

## chunk 191
data(gastricbypassW, package = "LMMstar")
colSums(is.na(gastricbypassW))

## chunk 192
library(mice)
set.seed(10)
gastricbypassW.mice <- mice(gastricbypassW, m = 5, printFlag = FALSE)
gastricbypassW.NNA <- complete(gastricbypassW.mice, action = "long")
table(gastricbypassW.NNA$.imp)

## chunk 193
e.mlmm <- mlmm(glucagonAUC3~glucagonAUC2+weight2, data=gastricbypassW.NNA,
               by = ".imp", effects = "weight2=0", trace = FALSE)
model.tables(e.mlmm)

## chunk 194
model.tables(e.mlmm, method = "pool.rubin")

## chunk 195
e.mice <- with(data=gastricbypassW.mice,exp=lm(glucagonAUC3~glucagonAUC2+weight2))
summary(pool(e.mice))

## chunk 196
plot(e.mlmm, method = c("pool.rubin","none"))

## * Data generation

## chunk 198
set.seed(10) ## ensure reproductibility
n.obs <- 100
n.times <- 4
mu <- rep(0,4)
gamma <- matrix(0, nrow = n.times, ncol = 10) ## add interaction
gamma[,6] <- c(0,1,1.5,1.5)
dW <- sampleRem(n.obs, n.times = n.times, mu = mu, gamma = gamma, format = "wide")
head(round(dW,3))

## chunk 199
set.seed(10) ## ensure reproductibility
dL <- sampleRem(n.obs, n.times = n.times, mu = mu, gamma = gamma, format = "long")
head(dL)

## * Modifying default options

## chunk 200
LMMstar.options("type.information")

## chunk 201
LMMstar.options(type.information = "expected")

## chunk 202
LMMstar.options(reinitialise = TRUE)

## * R session

## chunk 203
sessionInfo()

## * References
## * Likelihood in a linear mixed model
## ** Log-likelihood
## ** Score
## ** Hessian
## ** Degrees of freedom
## * Likelihood ratio test with the REML criterion

## chunk 204
## data(dfL, package = "LMMstar")
dfTest <- gastricbypassL[!is.na(gastricbypassL$glucagonAUC),]
dfTest$gluc <- dfTest$glucagonAUC
dfTest$gluc2 <- dfTest$glucagonAUC*2

## chunk 205
eML.UN <- lmm(weight ~ time+gluc, data = dfTest, repetition = ~time|id, method = "ML")
eML.UN2 <- lmm(weight ~ time+gluc, data = dfTest, repetition = ~time|id, method = "ML")
c(logLik(eML.UN), logLik(eML.UN2), logLik(eML.UN) - logLik(eML.UN2))

## chunk 206
eREML.UN <- lmm(weight ~ time + gluc, data = dfTest, repetition = ~time|id, method = "REML")
eREML.UN2 <- lmm(weight ~ time + gluc2, data = dfTest, repetition = ~time|id, method = "REML")
c(logLik(eREML.UN), logLik(eREML.UN2), logLik(eREML.UN) - logLik(eREML.UN2), log(2))

## chunk 207
set.seed(5) 
dfTest$ff <- rbinom(NROW(dfTest), size = 1, prob = 0.5)
logLik(lmm(weight ~ time+gluc, data = dfTest, repetition = ~time|id, method = "REML"))
logLik(lmm(weight ~ time+gluc*ff, data = dfTest, repetition = ~time|id, method = "REML"))

## chunk 208
logLik(lmm(weight ~ time + gluc, data = dfTest, repetition = ~time|id, method = "ML"))
logLik(lmm(weight ~ time + gluc*ff, data = dfTest, repetition = ~time|id, method = "ML"))

## * Sum of squares in a linear mixed model

## chunk 209
df.aov <- gastricbypassL[!is.na(gastricbypassL$glucagon),]

## chunk 210
e.lm <- lm(weight ~ visit + glucagonAUC, data = df.aov)
car::Anova(e.lm, type = "II")

## chunk 212
e.lmm <- lmm(weight ~ visit + glucagonAUC, data = df.aov)

## chunk 213
SSEstar <- crossprod(residuals(e.lmm, type = "normalized"))
c(SSEstar = SSEstar, SSE = SSEstar * sigma(e.lmm))

## chunk 214
df.residual(e.lmm)

## chunk 215
eBeta.lmm <- coef(e.lmm)
eVcov.lmm <- vcov(e.lmm, type.information = "expected")

SSRstar.glucagon <- eBeta.lmm[5] %*% solve(eVcov.lmm[5,5]) %*% eBeta.lmm[5] 
SSRstar.time <- eBeta.lmm[2:4] %*% solve(eVcov.lmm[2:4,2:4]) %*% eBeta.lmm[2:4] 
c(SSR.glucagon = SSRstar.glucagon * sigma(e.lmm),
  SSR.time = SSRstar.time * sigma(e.lmm),
  F.glucagon = SSRstar.glucagon,
  F.time = SSRstar.time/3)

## chunk 216
R2.glucagon <- SSRstar.glucagon/(SSRstar.glucagon+SSEstar)
R2.glucagon

## chunk 217
sign(coef(e.lmm)["glucagonAUC"])*sqrt(R2.glucagon)

## chunk 218
summary(partialCor(e.lmm, R2 = TRUE))

## * Equivalence with other R packages
## ** nlme package

## chunk 219
library(nlme)

## chunk 220
eRI.lmm <- lmm(weight ~ visit*group, structure = "RE",
               data = gastricbypassL, repetition = ~visit|id)
eCS.gls <- gls(weight ~ visit*group, correlation = corCompSymm(form=~visit|id),
               data = gastricbypassL, na.action = na.omit)
eCS.lme <- lme(weight ~ visit*group, random = ~1|id,
               data = gastricbypassL, na.action = na.omit)
logLik(eRI.lmm)
logLik(eCS.lme)
logLik(eCS.gls)

## chunk 221
range(ranef(eRI.lmm)-ranef(eCS.lme))

## chunk 222
eUN.gls <- gls(glucagonAUC ~ visit*group,
               correlation = corSymm(form=~as.numeric(visit)|id),
               weights = varIdent(form=~1|visit),
               data = gastricbypassL, na.action = na.omit)
logLik(eUN.gls)
logLik(eUN.lmm)

## ** lme4 package

## chunk 223
library(lme4)
library(lmerTest)

## chunk 224
eRI.lmer <- lmer(weight ~ visit*group + (1|id),
                 data = gastricbypassL)
logLik(eRI.lmer)
logLik(eRI.lmm)

## chunk 225
range(ranef(eRI.lmm)-ranef(eRI.lmer)$id)

## chunk 226
eNRI.lmm <- lmm(weight ~ visit*group, structure = RE(~(1|id/baseline)),
               data = gastricbypassL, repetition = ~visit|id)
eNRI.lmer <- lmer(weight ~ visit*group + (1|id/baseline),
                  data = gastricbypassL)
logLik(eNRI.lmer)
logLik(eNRI.lmm)

## chunk 227
eRanefNRI.lmm <- ranef(eNRI.lmm, format = "wide")
eRanefNRI.lmer <- ranef(eNRI.lmer)
## id
range(eRanefNRI.lmm$estimate-eRanefNRI.lmer$id)
## baseline
range(c(eRanefNRI.lmm$estimate.FALSE,eRanefNRI.lmm$estimate.TRUE)-ranef(eNRI.lmer)$`baseline:id`)

## chunk 228
eUN.lmer <- lmer(glucagonAUC ~ visit*group + (0 + visit|id),
                 data = gastricbypassL,
                 control = lmerControl(check.nobs.vs.nRE = "ignore"))
logLik(eUN.lmer)
logLik(eUN.lmm)

## chunk 229
anova(eUN.lmm)

## chunk 230
## only the last line is comparable
anova(eUN.lmer)

## chunk 231
data("Penicillin")
eCRI.lmer <- lmer(diameter ~ 1 + (1|plate) + (1|sample), Penicillin)
logLik(eCRI.lmer)

## chunk 232
Penicillin$index <- paste(Penicillin$sample,Penicillin$plate,sep=".")
Penicillin$id <- 1

eCRI.lmm <- lmm(diameter ~ 1 + (1|plate) + (1|sample), data = Penicillin)
logLik(eCRI.lmm)

## chunk 233
range(ranef(eCRI.lmm)$estimate-rbind(ranef(eCRI.lmer)$plate,ranef(eCRI.lmer)$sample))

## ** mmrm package

## chunk 234
library(mmrm)
e.mmrm <- mmrm(
  formula = FEV1 ~ RACE + SEX + ARMCD * AVISIT + us(AVISIT | USUBJID),
  data = fev_data
)

## chunk 235
e.lmm <- lmm(
  formula = FEV1 ~ RACE + SEX + ARMCD * AVISIT,
  repetition = ~ AVISIT | USUBJID, structure = "UN",
  data = fev_data, type.information = "expected"
)

## chunk 236
logLik(e.mmrm) - logLik(e.lmm)
range(coef(e.mmrm) - coef(e.lmm))
range(vcov(e.mmrm) - vcov(e.lmm))

## ** emmeans package

## chunk 240
gastricbypassLB$group2 <- gastricbypassLB$weight1>150

## chunk 241
eCS.lmm_2 <- lmm(glucagonAUC ~ visit*group2, repetition =~visit|id, structure = "CS", data = gastricbypassLB)
logLik(eCS.lmm_2)

## chunk 242
eRI.lmer_2 <- lmer(glucagonAUC ~ visit*group2 + (1|id), data = gastricbypassLB)
logLik(eRI.lmer_2)

## chunk 243
effects(eCS.lmm_2, variable = NULL)

## chunk 244
library(emmeans)
emmeans(eRI.lmer_2, specs=~visit)

## chunk 245
table(gastricbypassLB$group2)/NROW(gastricbypassLB)

## chunk 246
eCS.elmm_2 <- model.tables(effects(eCS.lmm_2, variable = "group2"))
eCS.elmm_2

## chunk 247
0.5*eCS.elmm_2[eCS.elmm_2$group2==FALSE,"estimate"]+0.5*eCS.elmm_2[eCS.elmm_2$group2==TRUE,"estimate"]

## chunk 248
0.8*eCS.elmm_2[eCS.elmm_2$group2==FALSE,"estimate"]+0.2*eCS.elmm_2[eCS.elmm_2$group2==TRUE,"estimate"]

## chunk 249
mu.group1 <-  as.double(coef(e.group)["(Intercept)"])
mu.group2 <-  as.double(coef(e.group)["(Intercept)"] + coef(e.group)["group2TRUE"])
p.group1 <- 14/20          ; p.group2 <- 6/20
c(emmeans = (mu.group1+mu.group2)/2, predict = mu.group1 * p.group1 + mu.group2 * p.group2)

## ** effectsize package (\(R^2\) or \(\eta^2\))

## chunk 250
library(effectsize)
eta_squared(eCS.lmer)
cat("\n")

## chunk 251
eCS.Wald <- anova(eCS.lmm)$multivariate
eCS.Wald$df.num*eCS.Wald$statistic/(eCS.Wald$df.num*eCS.Wald$statistic+eCS.Wald$df.denom)

## chunk 252
eUN.Wald <- anova(eUN.lmm)$multivariate
eUN.Wald$df.num*eUN.Wald$statistic/(eUN.Wald$df.num*eUN.Wald$statistic+eUN.Wald$df.denom)

## chunk 253
eta_squared(eUN.lmer)
cat("\n")

## ** MuMIn package (\(R^2\))

## chunk 254
library(MuMIn)
r.squaredGLMM(eCS.lmer)
cat("\n")

## chunk 255
sigmaW <- sigma(eCS.lmm)[1,1]-sigma(eCS.lmm)[1,2]

## chunk 256
sigmaB <- sigma(eCS.lmm)[1,2]

## chunk 257
sigma2_XB <- var(fitted(eCS.lmm))

## chunk 258
c(R2m = sigma2_XB/(sigmaW + sigmaB + sigma2_XB),
  R2c = (sigma2_XB + sigmaB)/(sigmaW + sigmaB + sigma2_XB))

## ** stats package (partial residuals)

## chunk 264
gastricbypassW$group <- as.factor(as.numeric(gastricbypassW$id)%%2)
eIID.lm <- lm(weight4 ~ group + weight1, data = gastricbypassW)
pRes.lm <- residuals(eIID.lm, type = "partial")
head(pRes.lm)

## chunk 265
eIID.lmm <- lmm(weight4 ~ group + weight1, data = gastricbypassW)
(residuals(eIID.lmm, type = "partial", variable = "group") - pRes.lm[,"group"])
(residuals(eIID.lmm, type = "partial", variable = "weight1") - pRes.lm[,"weight1"])

## chunk 266
coef(eIID.lm)["group1"] * mean(gastricbypassW$group=="1")
coef(eIID.lm)["weight1"] * mean(gastricbypassW$weight1)

## chunk 267
(residuals(eIID.lmm, type = "partial-center", variable = "weight1") - pRes.lm[,"weight1"])