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

## chunk 2
library(LMMstar)
library(ggplot2)

## * Univariate linear regression

## chunk 3
df1 <- data.frame(lifeExp = state.x77[,4],
                  illiteracy = state.x77[,3],
                  income = state.x77[,2]/1000,
                  murder = state.x77[,5],
                  edu = cut(state.x77[,6],c(0,50,60,100)))
head(df1,4)

## chunk 4
summarize(lifeExp + illiteracy + income + murder + edu ~ 1, data = df1,
          columns = c("observed","missing","mean","min","max","sd"))

## ** No interaction

## chunk 5
e.lm <- lmm(lifeExp ~ income + illiteracy + murder + edu, data = df1)
model.tables(e.lm)

## chunk 6
coef(e.lm) - coef(lm(lifeExp ~ income + illiteracy + murder + edu, data = df1))

## chunk 7
df1$pres <- residuals(e.lm, type = "partial", variable = c("(Intercept)","income"))
head(df1)

## chunk 9
c(69.05 - 0.17590 * 2.1 - (-0.27822) * 15.1,
  69.31 - 0.17590 * 1.5 - (-0.27822) * 11.3 - 0.77306)

## chunk 10
df1$pres2 <- residuals(e.lm, type = "partial", variable = c("(Intercept)","income"),
                       at = data.frame(illiteracy = 1.17, murder = 7.378, edu = "(50,60]"))
head(df1)

## chunk 11
c(69.05 - 0.17590 * (2.1-1.170) - (-0.27822) * (15.1-7.378) + 0.30141,
  69.31 - 0.17590 * (1.5-1.170) - (-0.27822) * (11.3-7.378) + 0.30141 - 0.77306)

## chunk 12
unique(df1$pres2 - df1$pres)

## chunk 13
gg.pres <- ggplot(df1) + geom_point(aes(x=income, y=pres))
gg.pres <- gg.pres + geom_abline(intercept = coef(e.lm)["(Intercept)"],
                                 slope = coef(e.lm)["income"])
gg.pres <- gg.pres + ggtitle("(B) partial residuals")
gg.pres

## chunk 14
plot(e.lm, type = "partial", variable = c("(Intercept)","income")) # C
plot(e.lm, type = "partial", variable = c("(Intercept)","income"),
     at = data.frame(illiteracy = 1.17, murder = 7.378, edu = "(50,60]")) # D

## chunk 15
gg.obs <- ggplot(df1) + geom_point(aes(x=income, y=lifeExp))
gg.obs <- gg.obs + ggtitle("(A) observed")
gg.obs

## chunk 17
ls.plot <- autoplot(e.lm, type = "partial", variable = c("(Intercept)","income"))
lapply(ls.plot, class)

## chunk 18
ls.plot$plot  + coord_cartesian(ylim=c(68,74))

## ** What about confidence intervals?

## chunk 19
plot(e.lm, type = "partial", variable = c("(Intercept)","income"), ci.alpha = 0.25) ## E

## chunk 20
pres.ci <- residuals(e.lm, type = "partial", variable = c("(Intercept)","income"),
                     keep.data = TRUE, fitted.ci = TRUE)
head(pres.ci)

## chunk 21
gg.pres + geom_ribbon(data = pres.ci, alpha = 0.25,
                      aes(ymin = fitted.lower, ymax = fitted.upper, x = income))

## chunk 23
plot(e.lm, type = "partial", variable = "income", ci.alpha = 0.25) ## F
plot(e.lm, type = "partial", variable = "murder", ci.alpha = 0.25) ## G

## ** Interaction with a categorical variable

## chunk 25
e.lmI <- lmm(lifeExp ~ income:edu + illiteracy + murder, data = df1)
model.tables(e.lmI)

## chunk 26
plot(e.lmI, type = "partial", variable = c("(Intercept)","income","edu")) ## H

## chunk 27
plot(e.lm, type = "partial", variable = c("(Intercept)","income","edu")) ## I

## chunk 29
residuals(e.lmI, type = "partial", variable = c("(Intercept)","income","edu"))[1:5]

## chunk 30
c(69.05 - 0.12870 * 2.1 - (-0.27940) * 15.1,
  69.31 - 0.12870 * 1.5 - (-0.27940) * 11.3)

## chunk 31
residuals(e.lmI, type = "partial", variable = c("(Intercept)","income"),
          at = data.frame(illiteracy = 1.17, murder = 7.378))[1:5]

## chunk 32
c(69.05 - 0.12870 * (2.1-1.170) - (-0.27940) * (15.1-7.378),
  69.31 - 0.12870 * (1.5-1.170) - (-0.27940) * (11.3-7.378))

## * Linear mixed model 

## chunk 33
data(armd.wide, package = "nlmeU")
armd.long <- reshape(armd.wide, direction ="long",
                     varying = paste0("visual",c(0,4,12,24,52)), times = c(0,4,12,24,52),
                     timevar = "week.num", v.names = "visual")
armd.long$week <- as.factor(armd.long$week.num)

## chunk 34
summarizeNA(armd.long)

## chunk 35
e.lmm <- lmm(visual ~ week*treat.f + lesion, data = armd.long, repetition =~week|subject)

## chunk 36
plot(e.lmm, facet = ~lesion, labeller = label_both, facet_ncol = 4)

## chunk 37
round(coef(e.lmm),2)

## chunk 38
plot(e.lmm, type = "partial", variable = c("(Intercept)","week","treat.f"),
     at = data.frame(lesion = 2))

## chunk 39
plot(e.lmm, type = "partial", variable = c("(Intercept)","week","treat.f"),
     facet =~week, facet_nrow = 1, time.var = "treat.f", color = FALSE,
     at = data.frame(lesion = 2))

## chunk 41
armd.long$pres <- residuals(e.lmm, type = "partial", 
                            variable = c("(Intercept)","week","treat.f"),
                            at = data.frame(lesion = 2))
head(armd.long)

## chunk 42
c(59 - (-3.19) * (3-2),
  65 - (-3.19) * (1-2))

## * R session

## chunk 43
sessionInfo()