diff --git a/code/depression.R b/code/depression.R new file mode 100644 index 0000000..df163f6 --- /dev/null +++ b/code/depression.R @@ -0,0 +1,157 @@ +# depression.R +# +# content: (1) Read data +# (2) Visualize data +# (3) Fit growth curve model +# +# input: data/reisby.txt +# output: -- +# +# last mod: 2026-04-20, NW + +library("lme4") +library("lattice") + +#----- (1) Read data ---------------------------------------------------------- +dat <- read.table("data/reisby.txt", header = TRUE) +dat$id <- factor(dat$id) +dat$diag <- factor(dat$diag, levels = c("nonen", "endog")) +dat <- na.omit(dat) # drop missing values + +boxplot(hamd ~ week, dat, + col = "#3CB4DC", + ylab = "HDRS score", + xlab = "Time (week)") + +#----- (2) Fitting random slope model ----------------------------------------- + +# random slope model +lme2 <- lmer(hamd ~ week + (week | id), dat, REML = FALSE) + +xyplot(hamd + predict(lme2) ~ week | id, data = dat, + type = c("p", "l", "g"), + ratio = "xy", + distribute.type = TRUE, + layout = c(11, 6), + ylab = "HDRS score", + xlab = "Time (week)") + +#--------------- (3) Fitting quadratic model ---------------------------------- + +# model with quadratic time trend +lme3 <- lmer(hamd ~ week + I(week^2) + (week + I(week^2) | id), dat, + REML = FALSE) + +xyplot(hamd + predict(lme3) ~ week | id, data = dat, + type = c("p", "l", "g"), + ratio = "xy", + distribute.type = TRUE, + layout = c(11, 6), + ylab = "HDRS score", + xlab = "Time (week)") + +#--------------- (4) Check random effects structure --------------- + +# Catterpillar plots + +dotplot(ranef(lme3), col = "#3CB4DC", + scales = list(x = list(relation = "free")))[[1]] + +# Shrinkage plots + +df <- coef(lmList(hamd ~ week_c + I(week_c^2) | id, dat)) +cc1 <- as.data.frame(coef(lme3reml)$id) +names(cc1) <- c("A", "B", "C") + +df <- cbind(df, cc1) +ff <- fixef(lme3reml) + + +## shrinkage intercept and week +with(df, + xyplot(`(Intercept)` ~ week_c, aspect = 1, + x1 = B, y1 = A, + panel = function(x, y, x1, y1, subscripts, ...) { + panel.grid(h = -1, v = -1) + x1 <- x1[subscripts] + y1 <- y1[subscripts] + larrows(x, y, x1, y1, type = "closed", length = 0.1, fill = "black", + angle = 15, ...) + lpoints(x, y, + pch = 16, + col = trellis.par.get("superpose.symbol")$col[2]) + lpoints(x1, y1, + pch = 16, + col = trellis.par.get("superpose.symbol")$col[1]) + lpoints(ff[2], ff[1], + pch = 16, + col = trellis.par.get("superpose.symbol")$col[3]) + }, + xlab = "week_c", + ylab = "(Intercept)", + key = list(space = "top", columns = 3, + text = list(c("Mixed model", "Within-subject", "Population")), + points = list(col = trellis.par.get("superpose.symbol")$col[1:3], + pch = 16)) + ) +) + +## shrinkage intercept and week^2 +with(df, + xyplot(`(Intercept)` ~ `I(week_c^2)`, aspect = 1, + x1 = C, y1 = A, + panel = function(x, y, x1, y1, subscripts, ...) { + panel.grid(h = -1, v = -1) + x1 <- x1[subscripts] + y1 <- y1[subscripts] + larrows(x, y, x1, y1, type = "closed", length = 0.1, fill = "black", + angle = 15, ...) + lpoints(x, y, + pch = 16, + col = trellis.par.get("superpose.symbol")$col[2]) + lpoints(x1, y1, + pch = 16, + col = trellis.par.get("superpose.symbol")$col[1]) + lpoints(ff[3], ff[1], + pch = 16, + col = trellis.par.get("superpose.symbol")$col[3]) + }, + xlab = expression(week_c^2), + ylab = "(Intercept)", + key = list(space = "top", columns = 3, + text = list(c("Mixed model", "Within-subject", "Population")), + points = list(col = trellis.par.get("superpose.symbol")$col[1:3], + pch = 16)) + ) +) + +## shrinkage week and week^2 +with(df, + xyplot(week_c ~ `I(week_c^2)`, aspect = 1, + x1 = C, y1 = B, + panel = function(x, y, x1, y1, subscripts, ...) { + panel.grid(h = -1, v = -1) + x1 <- x1[subscripts] + y1 <- y1[subscripts] + larrows(x, y, x1, y1, type = "closed", length = 0.1, fill = "black", + angle = 15, ...) + lpoints(x, y, + pch = 16, + col = trellis.par.get("superpose.symbol")$col[2]) + lpoints(x1, y1, + pch = 16, + col = trellis.par.get("superpose.symbol")$col[1]) + lpoints(ff[3], ff[2], + pch = 16, + col = trellis.par.get("superpose.symbol")$col[3]) + }, + xlab = expression(week_c^2), + ylab = "week_c", + key = list(space = "top", columns = 3, + text = list(c("Mixed model", "Within-subject", "Population")), + points = list(col = trellis.par.get("superpose.symbol")$col[1:3], + pch = 16)) + ) +) + + diff --git a/data/moeller.csv b/data/moeller.csv new file mode 100644 index 0000000..99a5695 --- /dev/null +++ b/data/moeller.csv @@ -0,0 +1,2010 @@ +# Data from Moeller et al., International Journal of +# Neuropsychopharmacology (2008), 11, 985-997 +id,treat,week,pans +1,risp,0,85 +1,risp,1,91 +1,risp,2,86 +1,risp,3,77 +1,risp,4,69 +1,risp,5,59 +1,risp,6,69 +1,risp,7,60 +1,risp,8,64 +2,risp,0,77 +2,risp,1,69 +2,risp,2,59 +2,risp,3,52 +2,risp,4,67 +2,risp,5,59 +2,risp,6,56 +2,risp,7,56 +2,risp,8,61 +3,halo,0,66 +3,halo,1,59 +3,halo,2,66 +3,halo,3,55 +3,halo,4,66 +3,halo,5,61 +3,halo,6,48 +3,halo,7,57 +3,halo,8,57 +4,risp,0,99 +4,risp,1,86 +4,risp,2,66 +4,risp,3,58 +4,risp,4,57 +4,risp,5,57 +4,risp,6,57 +4,risp,7,57 +4,risp,8,57 +5,halo,0,64 +5,halo,1,57 +5,halo,2,46 +5,halo,3,50 +5,halo,4,52 +5,halo,5,52 +5,halo,6,52 +5,halo,7,52 +5,halo,8,52 +6,risp,0,88 +6,risp,1,81 +6,risp,2,83 +6,risp,3,71 +6,risp,4,49 +6,risp,5,72 +6,risp,6,57 +6,risp,7,74 +6,risp,8,74 +7,halo,0,112 +7,halo,1,79 +7,halo,2,93 +7,halo,3,93 +7,halo,4,93 +7,halo,5,93 +7,halo,6,93 +7,halo,7,93 +7,halo,8,93 +8,halo,0,68 +8,halo,1,30 +8,halo,2,82 +8,halo,3,74 +8,halo,4,103 +8,halo,5,103 +8,halo,6,103 +8,halo,7,103 +8,halo,8,103 +9,risp,0,105 +9,risp,1,126 +9,risp,2,126 +9,risp,3,126 +9,risp,4,126 +9,risp,5,126 +9,risp,6,126 +9,risp,7,126 +9,risp,8,126 +10,halo,0,84 +10,halo,1,56 +10,halo,2,64 +10,halo,3,55 +10,halo,4,53 +10,halo,5,44 +10,halo,6,44 +10,halo,7,44 +10,halo,8,44 +11,risp,0,88 +11,risp,1,88 +11,risp,2,93 +11,risp,3,94 +11,risp,4,97 +11,risp,5,99 +11,risp,6,99 +11,risp,7,99 +11,risp,8,99 +12,halo,0,121 +12,halo,1,109 +12,halo,2,105 +12,halo,3,82 +12,halo,4,75 +12,halo,5,75 +12,halo,6,75 +12,halo,7,73 +12,halo,8,73 +13,risp,0,111 +13,risp,2,52 +13,risp,3,45 +13,risp,4,45 +13,risp,5,45 +13,risp,6,45 +13,risp,7,45 +13,risp,8,57 +14,risp,0,118 +14,risp,1,106 +14,risp,2,106 +14,risp,3,106 +14,risp,4,106 +14,risp,5,106 +14,risp,6,106 +14,risp,7,106 +14,risp,8,106 +15,halo,0,102 +15,halo,1,94 +15,halo,2,90 +15,halo,3,89 +15,halo,4,96 +15,halo,5,96 +15,halo,6,96 +15,halo,7,96 +15,halo,8,96 +16,halo,0,96 +16,halo,1,63 +16,halo,2,50 +16,halo,3,73 +16,halo,4,81 +16,halo,5,92 +16,halo,6,92 +16,halo,7,92 +16,halo,8,92 +17,halo,0,73 +17,halo,1,41 +17,halo,2,41 +17,halo,3,41 +17,halo,4,41 +17,halo,5,41 +17,halo,6,41 +17,halo,7,41 +17,halo,8,41 +18,halo,0,84 +18,halo,1,40 +18,halo,2,45 +18,halo,3,35 +18,halo,4,34 +18,halo,5,37 +18,halo,6,32 +18,halo,7,32 +18,halo,8,32 +19,halo,0,74 +19,halo,1,65 +19,halo,2,38 +19,halo,3,31 +19,halo,4,33 +19,halo,5,33 +19,halo,6,33 +19,halo,7,33 +19,halo,8,33 +20,risp,0,64 +20,risp,1,30 +20,risp,2,36 +20,risp,3,34 +20,risp,4,39 +20,risp,5,45 +20,risp,6,41 +20,risp,7,45 +20,risp,8,41 +21,halo,0,52 +21,halo,1,38 +21,halo,2,42 +21,halo,3,41 +21,halo,4,36 +21,halo,5,36 +21,halo,6,43 +21,halo,7,41 +21,halo,8,41 +22,risp,0,94 +22,risp,1,72 +22,risp,2,73 +22,risp,3,72 +22,risp,4,69 +22,risp,5,74 +22,risp,6,81 +22,risp,7,89 +22,risp,8,87 +23,risp,0,89 +23,risp,1,35 +23,risp,2,34 +23,risp,3,52 +23,risp,4,58 +23,risp,5,54 +23,risp,6,54 +23,risp,7,54 +23,risp,8,54 +24,halo,0,72 +24,halo,1,64 +24,halo,2,63 +24,halo,3,65 +24,halo,4,71 +24,halo,5,73 +24,halo,6,62 +24,halo,7,62 +24,halo,8,62 +25,halo,0,87 +25,halo,1,95 +25,halo,2,77 +25,halo,3,88 +25,halo,4,70 +25,halo,5,57 +25,halo,6,70 +25,halo,7,71 +25,halo,8,72 +26,risp,0,119 +26,risp,1,132 +26,risp,2,114 +26,risp,3,88 +26,risp,4,83 +26,risp,5,72 +26,risp,6,65 +26,risp,7,65 +26,risp,8,69 +27,halo,0,109 +27,halo,1,93 +27,halo,2,84 +27,halo,3,90 +27,halo,4,82 +27,halo,5,91 +27,halo,6,92 +27,halo,7,92 +27,halo,8,92 +28,halo,0,106 +28,halo,1,107 +28,halo,2,94 +28,halo,3,70 +28,halo,4,70 +28,halo,5,70 +28,halo,6,70 +28,halo,7,70 +28,halo,8,70 +29,halo,0,92 +29,halo,1,90 +29,halo,2,96 +29,halo,3,96 +29,halo,4,96 +29,halo,5,96 +29,halo,6,96 +29,halo,7,96 +29,halo,8,96 +30,halo,0,77 +30,halo,1,88 +30,halo,2,63 +30,halo,3,67 +30,halo,4,67 +30,halo,5,67 +30,halo,6,67 +30,halo,7,67 +30,halo,8,67 +31,risp,0,62 +31,risp,1,47 +31,risp,2,47 +31,risp,3,47 +31,risp,4,39 +31,risp,5,41 +31,risp,6,41 +31,risp,7,41 +31,risp,8,41 +32,risp,0,59 +32,risp,1,52 +32,risp,2,51 +32,risp,3,49 +32,risp,4,43 +32,risp,5,43 +32,risp,6,39 +32,risp,7,36 +32,risp,8,36 +33,halo,0,100 +33,halo,1,99 +33,halo,2,90 +33,halo,3,68 +33,halo,4,74 +33,halo,5,64 +33,halo,6,55 +33,halo,7,45 +33,halo,8,49 +34,risp,0,99 +34,risp,2,115 +34,risp,3,72 +34,risp,4,68 +34,risp,5,68 +34,risp,6,68 +34,risp,7,56 +34,risp,8,55 +35,risp,0,97 +35,risp,1,66 +35,risp,2,64 +35,risp,3,68 +35,risp,4,53 +35,risp,5,56 +35,risp,6,56 +35,risp,7,53 +35,risp,8,54 +36,halo,0,93 +36,halo,1,71 +36,halo,2,71 +36,halo,3,69 +36,halo,4,64 +36,halo,5,65 +36,halo,6,72 +36,halo,7,72 +36,halo,8,72 +37,risp,0,98 +37,risp,1,91 +37,risp,2,97 +37,risp,3,97 +37,risp,4,97 +37,risp,5,97 +37,risp,6,97 +37,risp,7,97 +37,risp,8,97 +38,halo,0,113 +38,halo,1,105 +38,halo,2,55 +38,halo,3,42 +38,halo,4,37 +38,halo,5,38 +38,halo,6,35 +38,halo,7,35 +38,halo,8,35 +39,halo,0,74 +39,halo,1,58 +39,halo,2,65 +39,halo,3,54 +39,halo,4,49 +39,halo,5,49 +39,halo,6,49 +39,halo,7,49 +39,halo,8,49 +40,risp,0,117 +40,risp,1,113 +40,risp,2,109 +40,risp,3,96 +40,risp,4,84 +40,risp,5,83 +40,risp,6,77 +40,risp,7,70 +40,risp,8,69 +41,halo,0,90 +41,halo,1,51 +41,halo,2,49 +41,halo,3,49 +41,halo,4,74 +41,halo,5,74 +41,halo,6,74 +41,halo,7,74 +41,halo,8,74 +42,halo,0,79 +42,halo,1,89 +42,halo,2,69 +42,halo,3,66 +42,halo,4,57 +42,halo,5,57 +42,halo,6,55 +42,halo,7,59 +42,halo,8,40 +43,risp,0,94 +43,risp,1,76 +43,risp,2,72 +43,risp,3,62 +43,risp,4,50 +43,risp,5,39 +43,risp,6,37 +43,risp,7,32 +43,risp,8,32 +44,risp,0,88 +44,risp,1,85 +44,risp,2,92 +44,risp,3,78 +44,risp,4,75 +44,risp,5,64 +44,risp,6,68 +44,risp,7,70 +44,risp,8,68 +45,risp,0,74 +45,risp,1,49 +45,risp,2,48 +45,risp,3,48 +45,risp,4,40 +45,risp,5,43 +45,risp,6,49 +45,risp,7,53 +45,risp,8,47 +46,halo,0,80 +46,halo,1,73 +46,halo,2,74 +46,halo,3,60 +46,halo,4,65 +46,halo,5,49 +46,halo,6,46 +46,halo,7,49 +46,halo,8,46 +47,halo,0,82 +47,halo,1,53 +47,halo,2,48 +47,halo,3,45 +47,halo,4,56 +47,halo,5,43 +47,halo,6,39 +47,halo,7,39 +47,halo,8,39 +48,risp,0,98 +48,risp,1,95 +48,risp,2,48 +48,risp,3,40 +48,risp,4,53 +48,risp,5,50 +48,risp,6,50 +48,risp,7,56 +48,risp,8,59 +49,halo,0,78 +49,halo,1,87 +49,halo,2,77 +49,halo,3,58 +49,halo,4,44 +49,halo,5,48 +49,halo,6,45 +49,halo,7,39 +49,halo,8,41 +50,halo,0,54 +50,halo,1,44 +50,halo,2,43 +50,halo,3,54 +50,halo,4,47 +50,halo,5,45 +50,halo,6,40 +50,halo,7,50 +50,halo,8,43 +51,risp,0,93 +51,risp,1,83 +51,risp,2,79 +51,risp,3,77 +51,risp,4,55 +51,risp,5,61 +51,risp,6,67 +51,risp,7,56 +51,risp,8,61 +52,risp,0,60 +52,risp,1,63 +52,risp,2,51 +52,risp,3,48 +52,risp,4,41 +52,risp,5,38 +52,risp,6,36 +52,risp,7,37 +52,risp,8,40 +53,halo,0,89 +53,halo,1,89 +53,halo,2,62 +53,halo,3,62 +53,halo,4,62 +53,halo,5,62 +53,halo,6,75 +53,halo,7,75 +53,halo,8,75 +54,risp,0,68 +54,risp,1,53 +54,risp,2,52 +54,risp,3,38 +54,risp,4,35 +54,risp,5,37 +54,risp,6,35 +54,risp,7,37 +54,risp,8,34 +55,halo,0,88 +55,halo,1,110 +55,halo,2,119 +55,halo,3,110 +55,halo,4,115 +55,halo,5,115 +55,halo,6,100 +55,halo,7,109 +55,halo,8,109 +56,risp,0,42 +56,risp,1,45 +56,risp,2,43 +56,risp,3,36 +56,risp,4,34 +56,risp,5,35 +56,risp,6,37 +56,risp,7,39 +56,risp,8,39 +57,halo,0,58 +57,halo,1,39 +57,halo,2,30 +57,halo,3,30 +57,halo,4,30 +57,halo,5,30 +57,halo,6,30 +57,halo,7,30 +57,halo,8,30 +58,risp,0,32 +58,risp,1,30 +58,risp,2,35 +58,risp,3,35 +58,risp,4,35 +58,risp,5,35 +58,risp,6,35 +58,risp,7,35 +58,risp,8,35 +59,halo,0,41 +59,halo,1,38 +59,halo,2,32 +59,halo,3,34 +59,halo,4,32 +59,halo,5,30 +59,halo,6,30 +59,halo,7,30 +59,halo,8,30 +60,risp,0,40 +60,risp,1,69 +60,risp,2,32 +60,risp,3,30 +60,risp,4,30 +60,risp,5,31 +60,risp,6,31 +60,risp,7,31 +60,risp,8,31 +61,halo,0,73 +61,halo,1,64 +61,halo,2,62 +61,halo,3,63 +61,halo,4,61 +61,halo,5,63 +61,halo,6,62 +61,halo,7,62 +61,halo,8,62 +62,halo,0,85 +62,halo,1,72 +62,halo,2,61 +62,halo,3,40 +62,halo,4,91 +62,halo,5,99 +62,halo,6,99 +62,halo,7,101 +62,halo,8,101 +63,risp,0,55 +63,risp,1,42 +63,risp,2,44 +63,risp,3,53 +63,risp,4,53 +63,risp,5,58 +63,risp,6,51 +63,risp,7,55 +63,risp,8,73 +64,risp,0,59 +64,risp,1,56 +64,risp,2,78 +64,risp,3,78 +64,risp,4,78 +64,risp,5,78 +64,risp,6,78 +64,risp,7,78 +64,risp,8,78 +65,risp,0,59 +65,risp,1,56 +65,risp,2,58 +65,risp,3,52 +65,risp,4,66 +65,risp,5,64 +65,risp,6,66 +65,risp,7,59 +65,risp,8,57 +66,halo,0,63 +66,halo,1,62 +66,halo,2,69 +66,halo,3,60 +66,halo,4,60 +66,halo,5,71 +66,halo,6,71 +66,halo,7,71 +66,halo,8,71 +67,risp,0,92 +67,risp,1,76 +67,risp,2,86 +67,risp,3,86 +67,risp,4,80 +67,risp,5,68 +67,risp,6,75 +67,risp,7,73 +67,risp,8,73 +68,halo,0,98 +68,halo,1,77 +68,halo,2,85 +68,halo,3,78 +68,halo,4,79 +68,halo,5,79 +68,halo,6,83 +68,halo,7,84 +68,halo,8,84 +69,risp,0,58 +69,risp,1,32 +69,risp,2,49 +69,risp,3,49 +69,risp,4,49 +69,risp,5,49 +69,risp,6,49 +69,risp,7,49 +69,risp,8,49 +70,risp,0,65 +70,risp,1,54 +70,risp,2,54 +70,risp,3,54 +70,risp,4,54 +70,risp,5,54 +70,risp,6,54 +70,risp,7,54 +70,risp,8,54 +71,halo,0,64 +71,halo,1,51 +71,halo,2,48 +71,halo,3,38 +71,halo,4,38 +71,halo,5,38 +71,halo,6,38 +71,halo,7,33 +71,halo,8,35 +72,halo,0,51 +72,halo,1,48 +72,halo,2,48 +72,halo,3,46 +72,halo,4,40 +72,halo,5,40 +72,halo,6,73 +72,halo,7,73 +72,halo,8,73 +73,risp,0,41 +73,risp,1,36 +73,risp,2,37 +73,risp,3,49 +73,risp,4,47 +73,risp,5,47 +73,risp,6,42 +73,risp,7,42 +73,risp,8,42 +74,risp,0,38 +74,risp,1,32 +74,risp,2,32 +74,risp,3,33 +74,risp,4,33 +74,risp,5,33 +74,risp,6,33 +74,risp,7,33 +74,risp,8,33 +75,risp,0,61 +75,risp,3,49 +75,risp,4,42 +75,risp,5,42 +75,risp,6,42 +75,risp,7,57 +75,risp,8,39 +76,risp,0,87 +76,risp,1,81 +76,risp,2,78 +76,risp,3,71 +76,risp,4,60 +76,risp,5,57 +76,risp,6,56 +76,risp,7,56 +76,risp,8,46 +77,risp,0,82 +77,risp,1,52 +77,risp,2,50 +77,risp,3,40 +77,risp,4,45 +77,risp,5,40 +77,risp,6,34 +77,risp,7,34 +77,risp,8,41 +78,halo,0,101 +78,halo,1,98 +78,halo,2,93 +78,halo,3,78 +78,halo,4,78 +78,halo,5,75 +78,halo,6,45 +78,halo,7,36 +78,halo,8,47 +79,risp,0,74 +79,risp,1,61 +79,risp,2,42 +79,risp,3,43 +79,risp,4,43 +79,risp,5,38 +79,risp,6,38 +79,risp,7,38 +79,risp,8,38 +80,halo,0,82 +80,halo,1,58 +80,halo,2,46 +80,halo,3,40 +80,halo,4,43 +80,halo,5,42 +80,halo,6,42 +80,halo,7,39 +80,halo,8,34 +81,halo,0,47 +81,halo,1,45 +81,halo,2,46 +81,halo,3,50 +81,halo,4,43 +81,halo,5,58 +81,halo,6,47 +81,halo,7,45 +81,halo,8,39 +82,risp,0,72 +82,risp,1,43 +82,risp,2,42 +82,risp,3,48 +82,risp,4,48 +82,risp,5,41 +82,risp,6,45 +82,risp,7,44 +82,risp,8,44 +83,halo,0,67 +83,halo,1,52 +83,halo,2,51 +83,halo,3,42 +83,halo,4,37 +83,halo,5,55 +83,halo,6,56 +83,halo,7,56 +83,halo,8,56 +84,risp,0,82 +84,risp,1,58 +84,risp,2,53 +84,risp,3,42 +84,risp,4,36 +84,risp,5,39 +84,risp,6,39 +84,risp,7,39 +84,risp,8,39 +85,halo,0,60 +85,halo,1,43 +85,halo,2,48 +85,halo,3,46 +85,halo,4,47 +85,halo,5,42 +85,halo,6,42 +85,halo,7,42 +85,halo,8,42 +86,risp,0,54 +86,risp,1,47 +86,risp,2,38 +86,risp,3,40 +86,risp,4,36 +86,risp,5,34 +86,risp,6,34 +86,risp,7,34 +86,risp,8,34 +87,halo,0,57 +87,halo,1,52 +87,halo,2,42 +87,halo,3,58 +87,halo,4,58 +87,halo,5,58 +87,halo,6,58 +87,halo,7,58 +87,halo,8,58 +88,risp,0,69 +88,risp,1,51 +88,risp,2,42 +88,risp,3,36 +88,risp,4,38 +88,risp,5,40 +88,risp,6,40 +88,risp,7,40 +88,risp,8,40 +89,risp,0,64 +89,risp,1,61 +89,risp,2,67 +89,risp,3,74 +89,risp,4,70 +89,risp,5,77 +89,risp,6,70 +89,risp,7,78 +89,risp,8,71 +90,risp,0,92 +90,risp,2,82 +90,risp,3,82 +90,risp,4,92 +90,risp,5,83 +90,risp,6,83 +90,risp,7,83 +90,risp,8,83 +91,halo,0,84 +91,halo,1,66 +91,halo,2,83 +91,halo,3,69 +91,halo,4,69 +91,halo,5,70 +91,halo,6,51 +91,halo,7,54 +91,halo,8,55 +92,risp,0,63 +92,risp,1,40 +92,risp,2,36 +92,risp,3,36 +92,risp,4,37 +92,risp,5,37 +92,risp,6,37 +92,risp,7,37 +92,risp,8,37 +93,risp,0,48 +93,risp,1,46 +93,risp,2,50 +93,risp,3,40 +93,risp,4,40 +93,risp,5,40 +93,risp,6,39 +93,risp,7,38 +93,risp,8,42 +94,halo,0,67 +94,halo,1,65 +94,halo,2,76 +94,halo,3,63 +94,halo,4,53 +94,halo,5,51 +94,halo,6,51 +94,halo,7,51 +94,halo,8,51 +95,halo,0,66 +95,halo,1,44 +95,halo,2,46 +95,halo,3,50 +95,halo,4,58 +95,halo,5,60 +95,halo,6,60 +95,halo,7,60 +95,halo,8,60 +96,risp,0,45 +96,risp,1,47 +96,risp,2,43 +96,risp,3,41 +96,risp,4,52 +96,risp,5,42 +96,risp,6,38 +96,risp,7,38 +96,risp,8,43 +97,risp,0,68 +97,risp,1,63 +97,risp,2,54 +97,risp,3,54 +97,risp,4,65 +97,risp,5,68 +97,risp,6,49 +97,risp,7,48 +97,risp,8,48 +98,halo,0,92 +98,halo,1,86 +98,halo,2,71 +98,halo,3,53 +98,halo,4,60 +98,halo,5,66 +98,halo,6,66 +98,halo,7,66 +98,halo,8,66 +99,halo,0,52 +99,halo,1,49 +99,halo,2,42 +99,halo,3,52 +99,halo,4,48 +99,halo,5,57 +99,halo,6,58 +99,halo,7,58 +99,halo,8,58 +100,risp,0,91 +100,risp,1,74 +100,risp,2,60 +100,risp,3,50 +100,risp,4,55 +100,risp,5,49 +100,risp,6,52 +100,risp,7,53 +100,risp,8,50 +101,halo,0,50 +101,halo,1,44 +101,halo,2,44 +101,halo,3,41 +101,halo,4,36 +101,halo,5,32 +101,halo,6,32 +101,halo,7,42 +101,halo,8,34 +102,halo,0,70 +102,halo,1,44 +102,halo,2,57 +102,halo,3,57 +102,halo,4,57 +102,halo,5,62 +102,halo,6,60 +102,halo,7,59 +102,halo,8,54 +103,halo,0,56 +103,halo,1,56 +103,halo,2,54 +103,halo,3,57 +103,halo,4,53 +103,halo,5,47 +103,halo,6,50 +103,halo,7,50 +103,halo,8,50 +104,halo,0,78 +104,halo,1,40 +104,halo,2,41 +104,halo,3,37 +104,halo,4,33 +104,halo,5,39 +104,halo,6,44 +104,halo,7,42 +104,halo,8,45 +105,risp,0,66 +105,risp,1,44 +105,risp,2,43 +105,risp,3,33 +105,risp,4,30 +105,risp,5,30 +105,risp,6,30 +105,risp,7,30 +105,risp,8,30 +106,risp,0,70 +106,risp,1,53 +106,risp,2,56 +106,risp,3,50 +106,risp,4,37 +106,risp,5,37 +106,risp,6,40 +106,risp,7,38 +106,risp,8,43 +107,halo,0,67 +107,halo,1,44 +107,halo,2,59 +107,halo,3,58 +107,halo,4,43 +107,halo,5,44 +107,halo,6,37 +107,halo,7,34 +107,halo,8,36 +108,risp,0,71 +108,risp,1,80 +108,risp,2,66 +108,risp,3,52 +108,risp,4,54 +108,risp,5,54 +108,risp,6,54 +108,risp,7,54 +108,risp,8,54 +109,halo,0,79 +109,halo,1,40 +109,halo,2,34 +109,halo,3,36 +109,halo,4,32 +109,halo,5,32 +109,halo,6,32 +109,halo,7,32 +109,halo,8,30 +110,halo,0,54 +110,halo,1,45 +110,halo,2,42 +110,halo,3,43 +110,halo,4,44 +110,halo,5,42 +110,halo,6,43 +110,halo,7,41 +110,halo,8,41 +111,risp,0,103 +111,risp,1,68 +111,risp,2,67 +111,risp,3,61 +111,risp,4,49 +111,risp,5,56 +111,risp,6,51 +111,risp,7,51 +111,risp,8,49 +112,halo,0,84 +112,halo,1,71 +112,halo,2,59 +112,halo,3,58 +112,halo,4,58 +112,halo,5,55 +112,halo,6,46 +112,halo,7,42 +112,halo,8,46 +113,risp,0,68 +113,risp,1,56 +113,risp,2,40 +113,risp,3,45 +113,risp,4,44 +113,risp,5,44 +113,risp,6,42 +113,risp,7,42 +113,risp,8,44 +114,halo,0,47 +114,halo,1,58 +114,halo,2,47 +114,halo,3,43 +114,halo,4,43 +114,halo,5,46 +114,halo,6,42 +114,halo,7,37 +114,halo,8,40 +115,risp,0,67 +115,risp,1,68 +115,risp,2,60 +115,risp,3,59 +115,risp,4,47 +115,risp,5,50 +115,risp,6,49 +115,risp,7,57 +115,risp,8,53 +116,risp,0,68 +116,risp,1,52 +116,risp,2,49 +116,risp,3,46 +116,risp,4,40 +116,risp,5,40 +116,risp,6,40 +116,risp,7,40 +116,risp,8,32 +117,halo,0,44 +117,halo,1,42 +117,halo,2,43 +117,halo,3,54 +117,halo,4,48 +117,halo,5,42 +117,halo,6,42 +117,halo,7,40 +117,halo,8,38 +118,risp,0,76 +118,risp,1,37 +118,risp,2,36 +118,risp,3,34 +118,risp,4,34 +118,risp,5,32 +118,risp,6,32 +118,risp,7,32 +118,risp,8,30 +119,risp,0,81 +119,risp,1,74 +119,risp,2,69 +119,risp,3,69 +119,risp,4,69 +119,risp,5,69 +119,risp,6,69 +119,risp,7,69 +119,risp,8,69 +120,halo,0,48 +120,halo,1,34 +120,halo,2,30 +120,halo,3,30 +120,halo,4,30 +120,halo,5,30 +120,halo,6,30 +120,halo,7,30 +120,halo,8,30 +121,halo,0,56 +121,halo,1,47 +121,halo,2,36 +121,halo,3,32 +121,halo,4,38 +121,halo,5,31 +121,halo,6,31 +121,halo,7,33 +121,halo,8,32 +122,risp,0,76 +122,risp,1,54 +122,risp,2,71 +122,risp,3,51 +122,risp,4,37 +122,risp,5,36 +122,risp,6,41 +122,risp,7,40 +122,risp,8,40 +123,halo,0,67 +123,halo,1,72 +123,halo,2,46 +123,halo,3,39 +123,halo,4,36 +123,halo,5,32 +123,halo,6,34 +123,halo,7,35 +123,halo,8,37 +124,halo,0,59 +124,halo,1,42 +124,halo,2,38 +124,halo,3,38 +124,halo,4,34 +124,halo,5,32 +124,halo,6,32 +124,halo,7,30 +124,halo,8,30 +125,risp,0,62 +125,risp,1,56 +125,risp,2,53 +125,risp,3,42 +125,risp,4,42 +125,risp,5,41 +125,risp,6,50 +125,risp,7,44 +125,risp,8,44 +126,halo,0,70 +126,halo,1,43 +126,halo,2,54 +126,halo,3,36 +126,halo,4,40 +126,halo,5,42 +126,halo,6,41 +126,halo,7,41 +126,halo,8,56 +127,risp,0,60 +127,risp,1,54 +127,risp,2,51 +127,risp,3,59 +127,risp,4,59 +127,risp,5,59 +127,risp,6,59 +127,risp,7,59 +127,risp,8,59 +128,halo,0,47 +128,halo,1,40 +128,halo,2,36 +128,halo,3,37 +128,halo,4,34 +128,halo,5,34 +128,halo,6,34 +128,halo,7,34 +128,halo,8,34 +129,risp,0,63 +129,risp,1,56 +129,risp,2,43 +129,risp,3,44 +129,risp,4,44 +129,risp,5,44 +129,risp,6,44 +129,risp,7,44 +129,risp,8,44 +130,risp,0,56 +130,risp,1,47 +130,risp,2,40 +130,risp,3,48 +130,risp,4,48 +130,risp,5,44 +130,risp,6,52 +130,risp,7,56 +130,risp,8,43 +131,risp,0,67 +131,risp,1,66 +131,risp,2,58 +131,risp,3,45 +131,risp,4,52 +131,risp,5,47 +131,risp,6,46 +131,risp,7,38 +131,risp,8,41 +132,risp,0,50 +132,risp,1,50 +132,risp,2,46 +132,risp,3,65 +132,risp,4,61 +132,risp,5,55 +132,risp,6,62 +132,risp,7,64 +132,risp,8,63 +133,halo,0,76 +133,halo,1,59 +133,halo,2,54 +133,halo,3,44 +133,halo,4,44 +133,halo,5,44 +133,halo,6,44 +133,halo,7,44 +133,halo,8,44 +134,risp,0,57 +134,risp,1,46 +134,risp,2,49 +134,risp,3,56 +134,risp,4,68 +134,risp,5,61 +134,risp,6,53 +134,risp,7,48 +134,risp,8,50 +135,risp,0,65 +135,risp,1,67 +135,risp,2,44 +135,risp,3,55 +135,risp,4,44 +135,risp,5,44 +135,risp,6,30 +135,risp,7,50 +135,risp,8,43 +136,halo,0,66 +136,halo,1,54 +136,halo,2,44 +136,halo,3,42 +136,halo,4,42 +136,halo,5,42 +136,halo,6,60 +136,halo,7,54 +136,halo,8,50 +137,halo,0,45 +137,halo,1,38 +137,halo,2,35 +137,halo,3,31 +137,halo,4,31 +137,halo,5,38 +137,halo,6,38 +137,halo,7,38 +137,halo,8,38 +138,risp,0,59 +138,risp,1,50 +138,risp,2,44 +138,risp,3,50 +138,risp,4,66 +138,risp,5,66 +138,risp,6,66 +138,risp,7,66 +138,risp,8,66 +139,halo,0,82 +139,halo,1,82 +139,halo,2,80 +139,halo,3,68 +139,halo,4,64 +139,halo,5,52 +139,halo,6,46 +139,halo,7,46 +139,halo,8,46 +140,risp,0,59 +140,risp,1,44 +140,risp,2,35 +140,risp,3,41 +140,risp,4,53 +140,risp,5,63 +140,risp,6,67 +140,risp,7,66 +140,risp,8,66 +141,halo,0,75 +141,halo,1,59 +141,halo,2,59 +141,halo,3,61 +141,halo,4,54 +141,halo,5,40 +141,halo,6,40 +141,halo,7,35 +141,halo,8,43 +142,risp,0,44 +142,risp,1,45 +142,risp,2,45 +142,risp,3,45 +142,risp,4,45 +142,risp,5,45 +142,risp,6,45 +142,risp,7,45 +142,risp,8,45 +143,halo,0,62 +143,halo,1,56 +143,halo,2,60 +143,halo,3,57 +143,halo,4,56 +143,halo,5,58 +143,halo,6,56 +143,halo,7,51 +143,halo,8,58 +144,risp,0,63 +144,risp,1,56 +144,risp,2,53 +144,risp,3,58 +144,risp,4,52 +144,risp,5,51 +144,risp,6,69 +144,risp,7,49 +144,risp,8,55 +145,halo,0,48 +145,halo,1,46 +145,halo,2,42 +145,halo,3,46 +145,halo,4,44 +145,halo,5,44 +145,halo,6,44 +145,halo,7,44 +145,halo,8,44 +146,risp,0,45 +146,risp,1,45 +146,risp,2,31 +146,risp,3,30 +146,risp,4,30 +146,risp,5,30 +146,risp,6,30 +146,risp,7,30 +146,risp,8,30 +147,halo,0,72 +147,halo,1,33 +147,halo,2,33 +147,halo,3,33 +147,halo,4,33 +147,halo,5,33 +147,halo,6,33 +147,halo,7,33 +147,halo,8,33 +148,risp,0,40 +148,risp,1,43 +148,risp,2,42 +148,risp,3,40 +148,risp,4,39 +148,risp,5,39 +148,risp,6,39 +148,risp,7,39 +148,risp,8,39 +149,halo,0,62 +149,halo,1,54 +149,halo,2,53 +149,halo,3,51 +149,halo,4,46 +149,halo,5,46 +149,halo,6,46 +149,halo,7,46 +149,halo,8,46 +150,halo,0,58 +150,halo,1,30 +150,halo,2,40 +150,halo,3,40 +150,halo,4,40 +150,halo,5,40 +150,halo,6,40 +150,halo,7,40 +150,halo,8,40 +151,risp,0,70 +151,risp,1,70 +151,risp,2,70 +151,risp,3,70 +151,risp,4,70 +151,risp,5,70 +151,risp,6,38 +151,risp,7,38 +151,risp,8,41 +152,halo,0,46 +152,halo,1,43 +152,halo,2,42 +152,halo,3,50 +152,halo,4,50 +152,halo,5,45 +152,halo,6,45 +152,halo,7,45 +152,halo,8,45 +153,risp,0,69 +153,risp,1,41 +153,risp,2,38 +153,risp,3,33 +153,risp,4,45 +153,risp,5,38 +153,risp,6,36 +153,risp,7,35 +153,risp,8,37 +154,halo,0,46 +154,halo,1,30 +154,halo,2,34 +154,halo,3,30 +154,halo,4,30 +154,halo,5,34 +154,halo,6,35 +154,halo,7,39 +154,halo,8,44 +155,halo,0,69 +155,halo,1,34 +155,halo,2,34 +155,halo,3,33 +155,halo,4,30 +155,halo,5,30 +155,halo,6,31 +155,halo,7,34 +155,halo,8,34 +156,risp,0,60 +156,risp,1,56 +156,risp,2,51 +156,risp,3,42 +156,risp,4,63 +156,risp,5,48 +156,risp,6,48 +156,risp,7,35 +156,risp,8,37 +157,halo,0,82 +157,halo,1,45 +157,halo,2,45 +157,halo,3,47 +157,halo,4,41 +157,halo,5,43 +157,halo,6,43 +157,halo,7,58 +157,halo,8,55 +158,risp,0,64 +158,risp,1,41 +158,risp,2,48 +158,risp,3,46 +158,risp,4,46 +158,risp,5,47 +158,risp,6,47 +158,risp,7,47 +158,risp,8,47 +159,risp,0,40 +159,risp,1,38 +159,risp,2,43 +159,risp,3,43 +159,risp,4,43 +159,risp,5,43 +159,risp,6,43 +159,risp,7,43 +159,risp,8,43 +160,risp,0,77 +160,risp,1,57 +160,risp,2,45 +160,risp,3,45 +160,risp,4,45 +160,risp,5,45 +160,risp,6,45 +160,risp,7,45 +160,risp,8,45 +161,halo,0,37 +161,halo,1,40 +161,halo,2,35 +161,halo,3,35 +161,halo,4,35 +161,halo,5,35 +161,halo,6,35 +161,halo,7,35 +161,halo,8,35 +162,risp,0,57 +162,risp,3,50 +162,risp,4,50 +162,risp,5,50 +162,risp,6,56 +162,risp,7,53 +162,risp,8,53 +163,risp,0,38 +163,risp,2,31 +163,risp,3,31 +163,risp,4,31 +163,risp,5,31 +163,risp,6,31 +163,risp,7,31 +163,risp,8,31 +164,risp,0,51 +164,risp,1,42 +164,risp,2,35 +164,risp,3,37 +164,risp,4,34 +164,risp,5,34 +164,risp,6,34 +164,risp,7,34 +164,risp,8,34 +165,risp,0,84 +165,risp,1,55 +165,risp,2,46 +165,risp,3,46 +165,risp,4,37 +165,risp,5,37 +165,risp,6,37 +165,risp,7,37 +165,risp,8,37 +166,halo,0,83 +166,halo,1,50 +166,halo,2,41 +166,halo,3,47 +166,halo,4,37 +166,halo,5,37 +166,halo,6,37 +166,halo,7,37 +166,halo,8,37 +167,risp,0,72 +167,risp,1,62 +167,risp,2,59 +167,risp,3,58 +167,risp,4,58 +167,risp,5,58 +167,risp,6,58 +167,risp,7,58 +167,risp,8,58 +168,halo,0,59 +168,halo,1,67 +168,halo,2,67 +168,halo,3,67 +168,halo,4,67 +168,halo,5,67 +168,halo,6,67 +168,halo,7,67 +168,halo,8,67 +169,halo,0,62 +169,halo,1,57 +169,halo,2,59 +169,halo,3,60 +169,halo,4,63 +169,halo,5,63 +169,halo,6,63 +169,halo,7,63 +169,halo,8,63 +170,halo,0,57 +170,halo,1,50 +170,halo,2,36 +170,halo,3,31 +170,halo,4,31 +170,halo,5,37 +170,halo,6,32 +170,halo,7,35 +170,halo,8,31 +171,risp,0,62 +171,risp,1,81 +171,risp,2,58 +171,risp,3,58 +171,risp,4,59 +171,risp,5,68 +171,risp,6,62 +171,risp,7,42 +171,risp,8,45 +172,halo,0,85 +172,halo,1,41 +172,halo,2,57 +172,halo,3,53 +172,halo,4,45 +172,halo,5,45 +172,halo,6,34 +172,halo,7,34 +172,halo,8,33 +173,risp,0,60 +173,risp,1,49 +173,risp,2,42 +173,risp,3,36 +173,risp,4,32 +173,risp,5,54 +173,risp,6,45 +173,risp,7,50 +173,risp,8,50 +174,halo,0,85 +174,halo,1,81 +174,halo,2,60 +174,halo,3,59 +174,halo,4,51 +174,halo,5,52 +174,halo,6,76 +174,halo,7,77 +174,halo,8,80 +175,halo,0,50 +175,halo,1,84 +175,halo,2,71 +175,halo,3,69 +175,halo,4,83 +175,halo,5,78 +175,halo,6,80 +175,halo,7,80 +175,halo,8,80 +176,risp,0,49 +176,risp,1,51 +176,risp,2,35 +176,risp,3,35 +176,risp,4,54 +176,risp,5,52 +176,risp,6,60 +176,risp,7,64 +176,risp,8,58 +177,risp,0,92 +177,risp,1,78 +177,risp,2,78 +177,risp,3,85 +177,risp,4,92 +177,risp,5,86 +177,risp,6,86 +177,risp,7,86 +177,risp,8,86 +178,halo,0,119 +178,halo,1,115 +178,halo,2,94 +178,halo,3,122 +178,halo,4,129 +178,halo,5,129 +178,halo,6,129 +178,halo,7,129 +178,halo,8,129 +179,risp,0,104 +179,risp,1,113 +179,risp,2,96 +179,risp,3,84 +179,risp,4,77 +179,risp,5,74 +179,risp,6,67 +179,risp,7,72 +179,risp,8,70 +180,risp,0,75 +180,risp,1,72 +180,risp,2,72 +180,risp,3,72 +180,risp,4,72 +180,risp,5,72 +180,risp,6,72 +180,risp,7,72 +180,risp,8,72 +181,halo,0,56 +181,halo,1,38 +181,halo,2,36 +181,halo,3,33 +181,halo,4,44 +181,halo,5,47 +181,halo,6,48 +181,halo,7,45 +181,halo,8,45 +182,risp,0,61 +182,risp,1,38 +182,risp,2,39 +182,risp,3,39 +182,risp,4,39 +182,risp,5,39 +182,risp,6,39 +182,risp,7,39 +182,risp,8,39 +183,risp,0,81 +183,risp,1,31 +183,risp,2,43 +183,risp,3,30 +183,risp,4,30 +183,risp,5,30 +183,risp,6,32 +183,risp,7,31 +183,risp,8,44 +184,risp,0,106 +184,risp,1,58 +184,risp,2,55 +184,risp,3,38 +184,risp,4,58 +184,risp,5,52 +184,risp,6,96 +184,risp,7,112 +184,risp,8,96 +185,halo,0,100 +185,halo,1,82 +185,halo,2,86 +185,halo,3,86 +185,halo,4,86 +185,halo,5,86 +185,halo,6,86 +185,halo,7,86 +185,halo,8,86 +186,risp,0,123 +186,risp,1,101 +186,risp,2,99 +186,risp,3,77 +186,risp,4,76 +186,risp,5,76 +186,risp,6,75 +186,risp,7,93 +186,risp,8,70 +187,risp,0,115 +187,risp,1,87 +187,risp,2,91 +187,risp,3,91 +187,risp,4,99 +187,risp,5,109 +187,risp,6,61 +187,risp,7,68 +187,risp,8,69 +188,risp,0,113 +188,risp,1,89 +188,risp,2,82 +188,risp,3,70 +188,risp,4,75 +188,risp,5,67 +188,risp,6,67 +188,risp,7,67 +188,risp,8,67 +189,halo,0,120 +189,halo,1,101 +189,halo,2,85 +189,halo,3,76 +189,halo,4,77 +189,halo,5,53 +189,halo,6,66 +189,halo,7,67 +189,halo,8,70 +190,halo,0,131 +190,halo,1,102 +190,halo,2,83 +190,halo,3,86 +190,halo,4,49 +190,halo,5,53 +190,halo,6,57 +190,halo,7,77 +190,halo,8,77 +191,halo,0,133 +191,halo,1,67 +191,halo,2,68 +191,halo,3,68 +191,halo,4,56 +191,halo,5,47 +191,halo,6,46 +191,halo,7,46 +191,halo,8,39 +192,risp,0,116 +192,risp,1,95 +192,risp,2,94 +192,risp,3,94 +192,risp,4,66 +192,risp,5,85 +192,risp,6,94 +192,risp,7,94 +192,risp,8,92 +193,halo,0,152 +193,halo,2,122 +193,halo,3,120 +193,halo,4,99 +193,halo,5,99 +193,halo,6,82 +193,halo,7,104 +193,halo,8,104 +194,halo,0,75 +194,halo,1,68 +194,halo,2,48 +194,halo,3,47 +194,halo,4,46 +194,halo,5,46 +194,halo,6,46 +194,halo,7,46 +194,halo,8,46 +195,risp,0,72 +195,risp,1,64 +195,risp,2,56 +195,risp,3,49 +195,risp,4,53 +195,risp,5,52 +195,risp,6,53 +195,risp,7,54 +195,risp,8,49 +196,halo,0,88 +196,halo,1,50 +196,halo,2,50 +196,halo,3,50 +196,halo,4,50 +196,halo,5,50 +196,halo,6,50 +196,halo,7,50 +196,halo,8,50 +197,halo,0,87 +197,halo,1,35 +197,halo,2,35 +197,halo,3,35 +197,halo,4,78 +197,halo,5,62 +197,halo,6,47 +197,halo,7,47 +197,halo,8,47 +198,risp,0,125 +198,risp,1,125 +198,risp,2,125 +198,risp,3,125 +198,risp,4,125 +198,risp,5,125 +198,risp,6,125 +198,risp,7,125 +198,risp,8,125 +199,halo,0,102 +199,halo,1,106 +199,halo,2,90 +199,halo,3,69 +199,halo,4,63 +199,halo,5,46 +199,halo,6,46 +199,halo,7,46 +199,halo,8,46 +200,risp,0,66 +200,risp,1,64 +200,risp,2,74 +200,risp,3,74 +200,risp,4,81 +200,risp,5,64 +200,risp,6,44 +200,risp,7,46 +200,risp,8,54 +201,halo,0,73 +201,halo,1,74 +201,halo,2,76 +201,halo,3,71 +201,halo,4,63 +201,halo,5,56 +201,halo,6,56 +201,halo,7,56 +201,halo,8,56 +202,halo,0,108 +202,halo,1,102 +202,halo,2,112 +202,halo,3,118 +202,halo,4,85 +202,halo,5,107 +202,halo,6,107 +202,halo,7,107 +202,halo,8,107 +203,halo,0,86 +203,halo,1,62 +203,halo,2,52 +203,halo,3,49 +203,halo,4,41 +203,halo,5,40 +203,halo,6,31 +203,halo,7,41 +203,halo,8,38 +204,halo,0,85 +204,halo,1,83 +204,halo,2,93 +204,halo,3,89 +204,halo,4,81 +204,halo,5,56 +204,halo,6,37 +204,halo,7,40 +204,halo,8,45 +205,halo,0,66 +205,halo,1,49 +205,halo,2,46 +205,halo,3,44 +205,halo,4,51 +205,halo,5,57 +205,halo,6,70 +205,halo,7,70 +205,halo,8,54 +206,halo,0,97 +206,halo,1,70 +206,halo,2,73 +206,halo,3,76 +206,halo,4,49 +206,halo,5,58 +206,halo,6,62 +206,halo,7,62 +206,halo,8,39 +207,risp,0,125 +207,risp,1,121 +207,risp,2,111 +207,risp,3,94 +207,risp,4,49 +207,risp,5,43 +207,risp,6,43 +207,risp,7,53 +207,risp,8,53 +208,risp,0,75 +208,risp,1,69 +208,risp,2,48 +208,risp,3,45 +208,risp,4,52 +208,risp,5,52 +208,risp,6,50 +208,risp,7,58 +208,risp,8,58 +209,risp,0,69 +209,risp,1,47 +209,risp,2,49 +209,risp,3,42 +209,risp,4,37 +209,risp,5,36 +209,risp,6,34 +209,risp,7,31 +209,risp,8,32 +210,halo,0,131 +210,halo,1,108 +210,halo,2,106 +210,halo,3,99 +210,halo,4,96 +210,halo,5,91 +210,halo,6,89 +210,halo,7,94 +210,halo,8,101 +211,risp,0,105 +211,risp,1,82 +211,risp,2,74 +211,risp,3,66 +211,risp,4,66 +211,risp,5,66 +211,risp,6,66 +211,risp,7,66 +211,risp,8,66 +212,halo,0,108 +212,halo,1,75 +212,halo,2,56 +212,halo,3,46 +212,halo,4,36 +212,halo,5,36 +212,halo,6,36 +212,halo,7,36 +212,halo,8,36 +213,risp,0,99 +213,risp,1,92 +213,risp,2,75 +213,risp,3,60 +213,risp,4,63 +213,risp,5,63 +213,risp,6,63 +213,risp,7,63 +213,risp,8,63 +214,halo,0,105 +214,halo,1,107 +214,halo,2,107 +214,halo,3,107 +214,halo,4,107 +214,halo,5,107 +214,halo,6,107 +214,halo,7,107 +214,halo,8,107 +215,risp,0,69 +215,risp,1,55 +215,risp,2,61 +215,risp,3,41 +215,risp,4,39 +215,risp,5,58 +215,risp,6,58 +215,risp,7,58 +215,risp,8,58 +216,risp,0,140 +216,risp,1,103 +216,risp,2,87 +216,risp,3,74 +216,risp,4,74 +216,risp,5,74 +216,risp,6,74 +216,risp,7,74 +216,risp,8,74 +217,risp,0,140 +217,risp,1,117 +217,risp,2,97 +217,risp,3,80 +217,risp,4,72 +217,risp,5,72 +217,risp,6,72 +217,risp,7,72 +217,risp,8,72 +218,halo,0,125 +218,halo,1,97 +218,halo,2,89 +218,halo,3,89 +218,halo,4,89 +218,halo,5,89 +218,halo,6,89 +218,halo,7,89 +218,halo,8,89 +219,halo,0,98 +219,halo,1,77 +219,halo,2,65 +219,halo,3,66 +219,halo,4,66 +219,halo,5,66 +219,halo,6,66 +219,halo,7,66 +219,halo,8,66 +220,halo,0,119 +220,halo,1,91 +220,halo,2,91 +220,halo,3,91 +220,halo,4,91 +220,halo,5,91 +220,halo,6,91 +220,halo,7,91 +220,halo,8,91 +221,halo,0,105 +221,halo,1,94 +221,halo,2,81 +221,halo,3,61 +221,halo,4,46 +221,halo,5,44 +221,halo,6,35 +221,halo,7,32 +221,halo,8,32 +222,risp,0,65 +222,risp,1,62 +222,risp,2,60 +222,risp,3,60 +222,risp,4,56 +222,risp,5,56 +222,risp,6,59 +222,risp,7,59 +222,risp,8,59 +223,halo,0,69 +223,halo,1,46 +223,halo,2,46 +223,halo,3,67 +223,halo,4,46 +223,halo,5,53 +223,halo,6,53 +223,halo,7,53 +223,halo,8,53 +224,risp,0,96 +224,risp,1,68 +224,risp,2,68 +224,risp,3,68 +224,risp,4,68 +224,risp,5,68 +224,risp,6,68 +224,risp,7,68 +224,risp,8,68 diff --git a/data/reisby.txt b/data/reisby.txt new file mode 100644 index 0000000..60f81d1 --- /dev/null +++ b/data/reisby.txt @@ -0,0 +1,399 @@ +## Data from Reisby et al., Psychopharmacology 54 (1977) 263-272 + id hamd week diag endweek +101 26 0 nonen 0 +101 22 1 nonen 0 +101 18 2 nonen 0 +101 7 3 nonen 0 +101 4 4 nonen 0 +101 3 5 nonen 0 +103 33 0 nonen 0 +103 24 1 nonen 0 +103 15 2 nonen 0 +103 24 3 nonen 0 +103 15 4 nonen 0 +103 13 5 nonen 0 +104 29 0 endog 0 +104 22 1 endog 1 +104 18 2 endog 2 +104 13 3 endog 3 +104 19 4 endog 4 +104 0 5 endog 5 +105 22 0 nonen 0 +105 12 1 nonen 0 +105 16 2 nonen 0 +105 16 3 nonen 0 +105 13 4 nonen 0 +105 9 5 nonen 0 +106 21 0 endog 0 +106 25 1 endog 1 +106 23 2 endog 2 +106 18 3 endog 3 +106 20 4 endog 4 +106 NA 5 endog 5 +107 21 0 endog 0 +107 21 1 endog 1 +107 16 2 endog 2 +107 19 3 endog 3 +107 NA 4 endog 4 +107 6 5 endog 5 +108 21 0 endog 0 +108 22 1 endog 1 +108 11 2 endog 2 +108 9 3 endog 3 +108 9 4 endog 4 +108 7 5 endog 5 +113 21 0 nonen 0 +113 23 1 nonen 0 +113 19 2 nonen 0 +113 23 3 nonen 0 +113 23 4 nonen 0 +113 NA 5 nonen 0 +114 NA 0 nonen 0 +114 17 1 nonen 0 +114 11 2 nonen 0 +114 13 3 nonen 0 +114 7 4 nonen 0 +114 7 5 nonen 0 +115 NA 0 endog 0 +115 16 1 endog 1 +115 16 2 endog 2 +115 16 3 endog 3 +115 16 4 endog 4 +115 11 5 endog 5 +117 19 0 endog 0 +117 16 1 endog 1 +117 13 2 endog 2 +117 12 3 endog 3 +117 7 4 endog 4 +117 6 5 endog 5 +118 NA 0 endog 0 +118 26 1 endog 1 +118 18 2 endog 2 +118 18 3 endog 3 +118 14 4 endog 4 +118 11 5 endog 5 +120 20 0 nonen 0 +120 19 1 nonen 0 +120 17 2 nonen 0 +120 18 3 nonen 0 +120 16 4 nonen 0 +120 17 5 nonen 0 +121 20 0 nonen 0 +121 22 1 nonen 0 +121 19 2 nonen 0 +121 19 3 nonen 0 +121 12 4 nonen 0 +121 14 5 nonen 0 +123 15 0 nonen 0 +123 15 1 nonen 0 +123 15 2 nonen 0 +123 13 3 nonen 0 +123 5 4 nonen 0 +123 5 5 nonen 0 +501 29 0 endog 0 +501 30 1 endog 1 +501 26 2 endog 2 +501 22 3 endog 3 +501 19 4 endog 4 +501 24 5 endog 5 +502 21 0 endog 0 +502 22 1 endog 1 +502 13 2 endog 2 +502 11 3 endog 3 +502 2 4 endog 4 +502 1 5 endog 5 +504 19 0 nonen 0 +504 17 1 nonen 0 +504 15 2 nonen 0 +504 16 3 nonen 0 +504 12 4 nonen 0 +504 12 5 nonen 0 +505 21 0 nonen 0 +505 11 1 nonen 0 +505 18 2 nonen 0 +505 0 3 nonen 0 +505 0 4 nonen 0 +505 4 5 nonen 0 +507 27 0 endog 0 +507 26 1 endog 1 +507 26 2 endog 2 +507 25 3 endog 3 +507 24 4 endog 4 +507 19 5 endog 5 +603 28 0 nonen 0 +603 22 1 nonen 0 +603 18 2 nonen 0 +603 20 3 nonen 0 +603 11 4 nonen 0 +603 13 5 nonen 0 +604 27 0 nonen 0 +604 27 1 nonen 0 +604 13 2 nonen 0 +604 5 3 nonen 0 +604 7 4 nonen 0 +604 NA 5 nonen 0 +606 19 0 endog 0 +606 33 1 endog 1 +606 12 2 endog 2 +606 12 3 endog 3 +606 3 4 endog 4 +606 1 5 endog 5 +607 30 0 endog 0 +607 39 1 endog 1 +607 30 2 endog 2 +607 27 3 endog 3 +607 20 4 endog 4 +607 4 5 endog 5 +608 24 0 nonen 0 +608 19 1 nonen 0 +608 14 2 nonen 0 +608 12 3 nonen 0 +608 3 4 nonen 0 +608 4 5 nonen 0 +609 NA 0 endog 1 +609 25 1 endog 1 +609 22 2 endog 2 +609 14 3 endog 3 +609 15 4 endog 4 +609 2 5 endog 5 +610 34 0 endog 0 +610 NA 1 endog 1 +610 33 2 endog 2 +610 23 3 endog 3 +610 NA 4 endog 4 +610 11 5 endog 5 +302 18 0 endog 0 +302 22 1 endog 1 +302 16 2 endog 2 +302 8 3 endog 3 +302 9 4 endog 4 +302 12 5 endog 5 +303 21 0 nonen 0 +303 21 1 nonen 0 +303 13 2 nonen 0 +303 14 3 nonen 0 +303 10 4 nonen 0 +303 5 5 nonen 0 +304 21 0 endog 0 +304 27 1 endog 1 +304 29 2 endog 2 +304 NA 3 endog 3 +304 12 4 endog 4 +304 24 5 endog 5 +305 19 0 nonen 0 +305 17 1 nonen 0 +305 15 2 nonen 0 +305 11 3 nonen 0 +305 5 4 nonen 0 +305 1 5 nonen 0 +308 22 0 nonen 0 +308 21 1 nonen 0 +308 18 2 nonen 0 +308 17 3 nonen 0 +308 12 4 nonen 0 +308 11 5 nonen 0 +309 22 0 nonen 0 +309 22 1 nonen 0 +309 16 2 nonen 0 +309 19 3 nonen 0 +309 20 4 nonen 0 +309 11 5 nonen 0 +310 24 0 endog 0 +310 19 1 endog 1 +310 11 2 endog 2 +310 7 3 endog 3 +310 6 4 endog 4 +310 NA 5 endog 5 +311 20 0 endog 0 +311 16 1 endog 1 +311 21 2 endog 2 +311 17 3 endog 3 +311 NA 4 endog 4 +311 15 5 endog 5 +312 17 0 endog 0 +312 NA 1 endog 1 +312 18 2 endog 2 +312 17 3 endog 3 +312 17 4 endog 4 +312 6 5 endog 5 +313 21 0 nonen 0 +313 19 1 nonen 0 +313 10 2 nonen 0 +313 11 3 nonen 0 +313 11 4 nonen 0 +313 8 5 nonen 0 +315 27 0 endog 0 +315 21 1 endog 1 +315 17 2 endog 2 +315 13 3 endog 3 +315 5 4 endog 4 +315 NA 5 endog 5 +316 32 0 endog 0 +316 26 1 endog 1 +316 23 2 endog 2 +316 26 3 endog 3 +316 23 4 endog 4 +316 24 5 endog 5 +318 17 0 endog 0 +318 18 1 endog 1 +318 19 2 endog 2 +318 21 3 endog 3 +318 17 4 endog 4 +318 11 5 endog 5 +319 24 0 endog 0 +319 18 1 endog 1 +319 10 2 endog 2 +319 14 3 endog 3 +319 13 4 endog 4 +319 12 5 endog 5 +322 28 0 endog 0 +322 21 1 endog 1 +322 25 2 endog 2 +322 32 3 endog 3 +322 34 4 endog 4 +322 NA 5 endog 5 +327 17 0 nonen 0 +327 18 1 nonen 0 +327 15 2 nonen 0 +327 8 3 nonen 0 +327 19 4 nonen 0 +327 17 5 nonen 0 +328 22 0 nonen 0 +328 24 1 nonen 0 +328 28 2 nonen 0 +328 26 3 nonen 0 +328 28 4 nonen 0 +328 29 5 nonen 0 +331 19 0 nonen 0 +331 21 1 nonen 0 +331 18 2 nonen 0 +331 16 3 nonen 0 +331 14 4 nonen 0 +331 10 5 nonen 0 +333 23 0 nonen 0 +333 20 1 nonen 0 +333 21 2 nonen 0 +333 20 3 nonen 0 +333 24 4 nonen 0 +333 14 5 nonen 0 +334 31 0 nonen 0 +334 25 1 nonen 0 +334 NA 2 nonen 0 +334 7 3 nonen 0 +334 8 4 nonen 0 +334 11 5 nonen 0 +335 21 0 nonen 0 +335 21 1 nonen 0 +335 18 2 nonen 0 +335 15 3 nonen 0 +335 12 4 nonen 0 +335 10 5 nonen 0 +337 27 0 nonen 0 +337 22 1 nonen 0 +337 23 2 nonen 0 +337 21 3 nonen 0 +337 12 4 nonen 0 +337 13 5 nonen 0 +338 22 0 nonen 0 +338 20 1 nonen 0 +338 22 2 nonen 0 +338 23 3 nonen 0 +338 19 4 nonen 0 +338 18 5 nonen 0 +339 27 0 endog 0 +339 NA 1 endog 1 +339 14 2 endog 2 +339 12 3 endog 3 +339 11 4 endog 4 +339 12 5 endog 5 +344 NA 0 endog 0 +344 21 1 endog 1 +344 12 2 endog 2 +344 13 3 endog 3 +344 13 4 endog 4 +344 18 5 endog 5 +345 29 0 nonen 0 +345 27 1 nonen 0 +345 27 2 nonen 0 +345 22 3 nonen 0 +345 22 4 nonen 0 +345 23 5 nonen 0 +346 25 0 endog 0 +346 24 1 endog 1 +346 19 2 endog 2 +346 23 3 endog 3 +346 14 4 endog 4 +346 21 5 endog 5 +347 18 0 endog 0 +347 15 1 endog 1 +347 14 2 endog 2 +347 10 3 endog 3 +347 8 4 endog 4 +347 NA 5 endog 5 +348 24 0 nonen 0 +348 21 1 nonen 0 +348 12 2 nonen 0 +348 13 3 nonen 0 +348 12 4 nonen 0 +348 5 5 nonen 0 +349 17 0 endog 0 +349 19 1 endog 1 +349 15 2 endog 2 +349 12 3 endog 3 +349 9 4 endog 4 +349 13 5 endog 5 +350 22 0 nonen 0 +350 25 1 nonen 0 +350 12 2 nonen 0 +350 16 3 nonen 0 +350 10 4 nonen 0 +350 16 5 nonen 0 +351 30 0 endog 0 +351 27 1 endog 1 +351 23 2 endog 2 +351 20 3 endog 3 +351 12 4 endog 4 +351 11 5 endog 5 +352 21 0 endog 0 +352 19 1 endog 1 +352 18 2 endog 2 +352 15 3 endog 3 +352 18 4 endog 4 +352 19 5 endog 5 +353 27 0 endog 0 +353 21 1 endog 1 +353 24 2 endog 2 +353 22 3 endog 3 +353 16 4 endog 4 +353 11 5 endog 5 +354 28 0 endog 0 +354 27 1 endog 1 +354 27 2 endog 2 +354 26 3 endog 3 +354 23 4 endog 4 +354 NA 5 endog 5 +355 22 0 endog 0 +355 26 1 endog 1 +355 20 2 endog 2 +355 13 3 endog 3 +355 10 4 endog 4 +355 7 5 endog 5 +357 27 0 endog 0 +357 22 1 endog 1 +357 24 2 endog 2 +357 25 3 endog 3 +357 19 4 endog 4 +357 19 5 endog 5 +360 21 0 endog 0 +360 28 1 endog 1 +360 27 2 endog 2 +360 29 3 endog 3 +360 28 4 endog 4 +360 33 5 endog 5 +361 30 0 endog 0 +361 22 1 endog 1 +361 11 2 endog 2 +361 8 3 endog 3 +361 7 4 endog 4 +361 19 5 endog 5 + diff --git a/exercises/schizo.md b/exercises/schizo.md new file mode 100644 index 0000000..75272e7 --- /dev/null +++ b/exercises/schizo.md @@ -0,0 +1,60 @@ +Exercise: Growth curve model with interaction +================ + +## Risperidone vs. haloperidol and schizophrenia (Möller et al. 2008) + +Analyze the data from [moeller.csv](../data/moeller.csv): + +- `pans`: Positive and Negative Symptom Scale for schizophrenia +- `treat`: medication group + - `risp`: atypical neuroleptic risperidone + - `halo`: conventional neuroleptic haloperidol + +``` r +dat <- read.table("../data/moeller.csv", header = TRUE, sep = ",") +dat$id <- factor(dat$id) +dat$treat <- factor(dat$treat, levels = c("risp", "halo")) + +lattice::xyplot(pans ~ week, data = dat, groups = treat, + type = c("g", "p", "a"), auto.key = TRUE) +``` + + + +- What is the sample size in each treatment group? + +- Estimate a random slope model + + - What are the estimates for the fixed effects and variance + components? + +- Fit a model with quadratic time trends for the population and + individual subjects + +- Add an interaction with treatment for the linear and the quadratic + effects for week and test them + +- Interpret your results: Which model would you choose? + +- Refit the model that you chose with a centered week variable: + + - Compare the estimates for your fixed effects and the covariance + components + - Which estimates change and why? + +### Reference + +
+ +
+ +Möller, H.-J., M. Riedel, M. Jäger, F. Wickelmaier, W. Maier, K.-U. +Kühn, G. Buchkremer, et al. 2008. “Short-Term Treatment with Risperidone +or Haloperidol in First-Episode Schizophrenia: 8-Week Results of a +Randomized Controlled Trial Within the German Research Network on +Schizophrenia.” *International Journal of Neuropsychopharmacology* 11 +(7): 985–97. . + +
+ +
diff --git a/exercises/schizo_files/figure-gfm/unnamed-chunk-2-1.png b/exercises/schizo_files/figure-gfm/unnamed-chunk-2-1.png new file mode 100644 index 0000000..a4cd55f Binary files /dev/null and b/exercises/schizo_files/figure-gfm/unnamed-chunk-2-1.png differ diff --git a/figures/cubic-gone-bad.pdf b/figures/cubic-gone-bad.pdf new file mode 100644 index 0000000..e5c6a03 Binary files /dev/null and b/figures/cubic-gone-bad.pdf differ diff --git a/figures/cubic.pdf b/figures/cubic.pdf new file mode 100644 index 0000000..de49f12 Binary files /dev/null and b/figures/cubic.pdf differ diff --git a/figures/hdrs-box.pdf b/figures/hdrs-box.pdf new file mode 100644 index 0000000..06e2450 Binary files /dev/null and b/figures/hdrs-box.pdf differ diff --git a/figures/hdrs-caterpillar.pdf b/figures/hdrs-caterpillar.pdf new file mode 100644 index 0000000..4545495 Binary files /dev/null and b/figures/hdrs-caterpillar.pdf differ diff --git a/figures/hdrs-ind_pred-quad.pdf b/figures/hdrs-ind_pred-quad.pdf new file mode 100644 index 0000000..9f17c72 Binary files /dev/null and b/figures/hdrs-ind_pred-quad.pdf differ diff --git a/figures/hdrs-ind_pred-randomslope.pdf b/figures/hdrs-ind_pred-randomslope.pdf new file mode 100644 index 0000000..bb38d1e Binary files /dev/null and b/figures/hdrs-ind_pred-randomslope.pdf differ diff --git a/figures/hdrs-plots.R b/figures/hdrs-plots.R new file mode 100644 index 0000000..5c48f2c --- /dev/null +++ b/figures/hdrs-plots.R @@ -0,0 +1,197 @@ +# hdrs-plots.R +# +# Example on growth curve models for Reisby et al. (1977) +# +# content: (1) Read data +# (2) Linear analysis +# (3) Fitting quadratic model +# (3) Centering of variables +# (4) Check random effects structure +# (5) Check model assumptions +# +# input: resiby.txt +# output: hdrs-ind_pred-randomslope.pdf +# hdrs-ind_pred-quad.pdf +# hdrs-caterpillar.pdf +# hdrs_shrinkage_int-week.pdf +# hdrs_shrinkage_int-weeksq.pdf +# hdrs_shrinkage_week-weeksq.pdf +# +# last mod: 2026-04-20, NW + +library("lme4") +library("lattice") + +#--------------- (1) Read data ------------------------------------------------ + +dat <- read.table("../data/reisby.txt", header = TRUE) +dat$id <- factor(dat$id) +dat$diag <- factor(dat$diag, levels = c("nonen", "endog")) +dat <- na.omit(dat) # drop missing values + +pdf("hdrs-box.pdf", height = 3.375, width = 3.375, pointsize = 10) + +par(mai = c(.6, .6, .1, .1), mgp = c(2, .4, 0)) +boxplot(hamd ~ week, dat, + col = "#3CB4DC", + ylab = "HDRS score", + xlab = "Time (week)") + +dev.off() + +#--------------- (2) Linear analysis ------------------------------------------ + +# random slope model +lme2 <- lmer(hamd ~ week + (week | id), dat, REML = FALSE) + +pdf("hdrs-ind_pred-randomslope.pdf", height = 6, width = 6) + +xyplot(hamd + predict(lme2) ~ week | id, data = dat, + type = c("p", "l", "g"), + #pch = 16, + ratio = "xy", + col = c("#3CB4DC", "#FF6900"), + distribute.type = TRUE, + layout = c(11, 6), + ylab = "HDRS score", + xlab = "Time (week)") + +dev.off() + +#--------------- (3) Fitting quadratic model ---------------------------------- + +# model with quadratic time trend +lme3 <- lmer(hamd ~ week + I(week^2) + (week + I(week^2) | id), dat, + REML = FALSE) + +pdf("hdrs-ind_pred-quad.pdf", height = 6, width = 6) + +xyplot(hamd + predict(lme3) ~ week | id, data = dat, + type = c("p", "l", "g"), + #pch = 16, + ratio = "xy", + col = c("#3CB4DC", "#FF6900"), + distribute.type = TRUE, + layout = c(11, 6), + ylab = "HDRS score", + xlab = "Time (week)") + +dev.off() + +#--------------- (4) Check random effects structure --------------- + +pdf("hdrs-caterpillar.pdf", height = 8, width = 20) + +dotplot(ranef(lme3), col = "#3CB4DC", + scales = list(x = list(relation = "free")))[[1]] + +dev.off() + +# Shrinkage plots + +df <- coef(lmList(hamd ~ week_c + I(week_c^2) | id, dat)) +cc1 <- as.data.frame(coef(lme3reml)$id) +names(cc1) <- c("A", "B", "C") + +df <- cbind(df, cc1) +ff <- fixef(lme3reml) + + +pdf("hdrs_shrinkage_int-week.pdf", height = 6, width = 6, pointsize = 10) + +with(df, + xyplot(`(Intercept)` ~ week_c, aspect = 1, + x1 = B, y1 = A, + panel = function(x, y, x1, y1, subscripts, ...) { + panel.grid(h = -1, v = -1) + x1 <- x1[subscripts] + y1 <- y1[subscripts] + larrows(x, y, x1, y1, type = "closed", length = 0.1, fill = "black", + angle = 15, ...) + lpoints(x, y, + pch = 16, + col = trellis.par.get("superpose.symbol")$col[2]) + lpoints(x1, y1, + pch = 16, + col = trellis.par.get("superpose.symbol")$col[1]) + lpoints(ff[2], ff[1], + pch = 16, + col = trellis.par.get("superpose.symbol")$col[3]) + }, + xlab = "week_c", + ylab = "(Intercept)", + key = list(space = "top", columns = 3, + text = list(c("Mixed model", "Within-subject", "Population")), + points = list(col = trellis.par.get("superpose.symbol")$col[1:3], + pch = 16)) + ) +) + +dev.off() + + +pdf("hdrs_shrinkage_int-weeksq.pdf", height = 6, width = 6, pointsize = 10) + +with(df, + xyplot(`(Intercept)` ~ `I(week_c^2)`, aspect = 1, + x1 = C, y1 = A, + panel = function(x, y, x1, y1, subscripts, ...) { + panel.grid(h = -1, v = -1) + x1 <- x1[subscripts] + y1 <- y1[subscripts] + larrows(x, y, x1, y1, type = "closed", length = 0.1, fill = "black", + angle = 15, ...) + lpoints(x, y, + pch = 16, + col = trellis.par.get("superpose.symbol")$col[2]) + lpoints(x1, y1, + pch = 16, + col = trellis.par.get("superpose.symbol")$col[1]) + lpoints(ff[3], ff[1], + pch = 16, + col = trellis.par.get("superpose.symbol")$col[3]) + }, + xlab = expression(week_c^2), + ylab = "(Intercept)", + key = list(space = "top", columns = 3, + text = list(c("Mixed model", "Within-subject", "Population")), + points = list(col = trellis.par.get("superpose.symbol")$col[1:3], + pch = 16)) + ) +) + +dev.off() + + +pdf("hdrs_shrinkage_week-weeksq.pdf", height = 6, width = 6, pointsize = 10) + +with(df, + xyplot(week_c ~ `I(week_c^2)`, aspect = 1, + x1 = C, y1 = B, + panel = function(x, y, x1, y1, subscripts, ...) { + panel.grid(h = -1, v = -1) + x1 <- x1[subscripts] + y1 <- y1[subscripts] + larrows(x, y, x1, y1, type = "closed", length = 0.1, fill = "black", + angle = 15, ...) + lpoints(x, y, + pch = 16, + col = trellis.par.get("superpose.symbol")$col[2]) + lpoints(x1, y1, + pch = 16, + col = trellis.par.get("superpose.symbol")$col[1]) + lpoints(ff[3], ff[2], + pch = 16, + col = trellis.par.get("superpose.symbol")$col[3]) + }, + xlab = expression(week_c^2), + ylab = "week_c", + key = list(space = "top", columns = 3, + text = list(c("Mixed model", "Within-subject", "Population")), + points = list(col = trellis.par.get("superpose.symbol")$col[1:3], + pch = 16)) + ) +) + +dev.off() + diff --git a/figures/hdrs-quad.pdf b/figures/hdrs-quad.pdf new file mode 100644 index 0000000..f469f68 Binary files /dev/null and b/figures/hdrs-quad.pdf differ diff --git a/figures/hdrs_shrinkage_int-week.pdf b/figures/hdrs_shrinkage_int-week.pdf new file mode 100644 index 0000000..32680b2 Binary files /dev/null and b/figures/hdrs_shrinkage_int-week.pdf differ diff --git a/figures/hdrs_shrinkage_int-weeksq.pdf b/figures/hdrs_shrinkage_int-weeksq.pdf new file mode 100644 index 0000000..bcacaa3 Binary files /dev/null and b/figures/hdrs_shrinkage_int-weeksq.pdf differ diff --git a/figures/hdrs_shrinkage_week-weeksq.pdf b/figures/hdrs_shrinkage_week-weeksq.pdf new file mode 100644 index 0000000..b7ee2cf Binary files /dev/null and b/figures/hdrs_shrinkage_week-weeksq.pdf differ diff --git a/figures/iwm_logo_rgb.pdf b/figures/iwm_logo_rgb.pdf new file mode 100644 index 0000000..d08238f Binary files /dev/null and b/figures/iwm_logo_rgb.pdf differ diff --git a/figures/quad.pdf b/figures/quad.pdf new file mode 100644 index 0000000..d3c4a7a Binary files /dev/null and b/figures/quad.pdf differ diff --git a/figures/sleep_box.pdf b/figures/sleep_box.pdf new file mode 100644 index 0000000..eb0a3a6 Binary files /dev/null and b/figures/sleep_box.pdf differ diff --git a/figures/sleep_random_intercept.pdf b/figures/sleep_random_intercept.pdf new file mode 100644 index 0000000..5e799a4 Binary files /dev/null and b/figures/sleep_random_intercept.pdf differ diff --git a/figures/sleep_random_slope.pdf b/figures/sleep_random_slope.pdf new file mode 100644 index 0000000..b0b9f73 Binary files /dev/null and b/figures/sleep_random_slope.pdf differ diff --git a/figures/sleep_shrinkfit.pdf b/figures/sleep_shrinkfit.pdf new file mode 100644 index 0000000..dad654f Binary files /dev/null and b/figures/sleep_shrinkfit.pdf differ diff --git a/figures/sleep_subjects.pdf b/figures/sleep_subjects.pdf new file mode 100644 index 0000000..063eda8 Binary files /dev/null and b/figures/sleep_subjects.pdf differ diff --git a/slides/lead_longitudinal.pdf b/slides/lead_longitudinal.pdf new file mode 100644 index 0000000..337b067 Binary files /dev/null and b/slides/lead_longitudinal.pdf differ diff --git a/slides/lead_longitudinal.tex b/slides/lead_longitudinal.tex new file mode 100644 index 0000000..ce186aa --- /dev/null +++ b/slides/lead_longitudinal.tex @@ -0,0 +1,789 @@ +\documentclass[aspectratio=169]{beamer} + +\usepackage{listings} +\usepackage[utf8,latin1]{inputenc} +\usepackage[style = apa, backend = biber, natbib = true]{biblatex} +\usepackage{tikz} +\addbibresource{lit.bib} + +\makeatletter \def\newblock{\beamer@newblock} \makeatother + +\beamertemplatenavigationsymbolsempty +\setbeamertemplate{itemize items}[circle] +\setbeamertemplate{section in toc}[circle] +\mode{\setbeamercolor{math text displayed}{fg=iwmgray}} +\setbeamercolor{block body}{bg=iwmorange!50!white} +\setbeamercolor{block title}{fg=white, bg=iwmorange} + +% Definitions for biblatex +\setbeamercolor{bibliography entry note}{fg=iwmgray} +\setbeamercolor{bibliography entry author}{fg=iwmgray} +\setbeamertemplate{bibliography item}{} + +\definecolor{iwmorange}{RGB}{255,105,0} +\definecolor{iwmgray}{RGB}{67,79,79} +\definecolor{iwmblue}{RGB}{60,180,220} +\definecolor{iwmgreen}{RGB}{145,200,110} +\definecolor{iwmpurple}{RGB}{120,0,75} + +\setbeamercolor{title}{fg=iwmorange} +\setbeamercolor{frametitle}{fg=iwmorange} +\setbeamercolor{structure}{fg=iwmorange} +\setbeamercolor{normal text}{fg=iwmgray} +\setbeamercolor{author}{fg=iwmgray} +\setbeamercolor{date}{fg=iwmgray} + +\lstset{language = R,% + basicstyle = \ttfamily\color{iwmgray}, + frame = single, + rulecolor = \color{iwmgray}, + commentstyle = \slshape\color{iwmgreen}, + keywordstyle = \bfseries\color{iwmgray}, + identifierstyle = \color{iwmpurple}, + stringstyle = \color{iwmblue}, + numbers = none,%left,numberstyle = \tiny, + basewidth = {.5em, .4em}, + showstringspaces = false, + emphstyle = \color{red!50!white}} + +\title{Introduction to mixed-effects models} +\subtitle{(for longitudinal data)} +\author{Nora Wickelmaier} +\institute{\includegraphics[scale=.2]{../figures/iwm_logo_rgb}} +\date{2026-04-28} + +\AtBeginSection[]{ + \frame{ + \tableofcontents[sectionstyle=show/hide, subsectionstyle=show/show/hide]}} + +\setbeamertemplate{headline}{ + \begin{beamercolorbox}{section in head} + \vskip5pt\insertsectionnavigationhorizontal{\paperwidth}{}{}\vskip2pt + \end{beamercolorbox} +} + +\setbeamertemplate{footline}{\vskip-2pt\hfill\insertframenumber$\;$\vskip2pt} + +\begin{document} + +\begin{frame}{} +\thispagestyle{empty} +\titlepage +\end{frame} + +\begin{frame}{Plan for today} + \begin{itemize}[<+->] + \item We will walk through an example for a longitudinal data set (6 time + points per subject) + \item I will explain the general concepts with the slides + \item We will switch to R and use the lme4 package to fit the models + \item You will use R to fit the model and try out different specifications + \item We will discuss the results + \item All the materials are here: + \url{https://gitea.iwm-tuebingen.de/nwickelmaier/lead_longitudinal} + \\~\\ + \item[$\to$] Try to code along in R! Ask as many questions as possible + \end{itemize} +\end{frame} + +\begin{frame}{Outline} +\tableofcontents +\end{frame} + +\section[Introduction]{Introduction to longitudinal data} + +\begin{frame}[fragile]{Longitudinal data} + \begin{itemize} + \item Consist of repeated measurements on the same subject taken over time + \item Are a frequent use case for mixed-effects models + \item Contain time as a predictor: + time trends within and between subjects are of interest + \end{itemize} + \begin{lstlisting} +library("lme4") +data("sleepstudy") +?sleepstudy +str(sleepstudy) +summary(sleepstudy) +head(sleepstudy) + \end{lstlisting} +\nocite{Bates10lme4, Bates15} +\end{frame} + +\begin{frame}{Sleep study} + \begin{itemize} + \item Average reaction time per day for subjects in a sleep deprivation + study + \item On day 0, the subjects had their normal amount of sleep + \item Starting that night they were restricted to 3 hours of sleep per + night + \item Observations represent the average reaction time on a series of + tests given each day to each subject + \end{itemize} + + \vfill + A data frame with 180 observations on the following 3 variables\\~\\ + + \begin{tabular}{ll} + \hline + \texttt{Reaction} & Average reaction time (ms) \\ + \texttt{Days} & Number of days of sleep deprivation \\ + \texttt{Subject} & Subject number on which the observation was made \\ + \hline + \end{tabular} +\end{frame} + + +\begin{frame}[fragile]{Visualization of data} + \begin{columns} + \begin{column}{.53\textwidth} + \includegraphics[scale=.8]{../figures/sleep_box} + \end{column} + \begin{column}{.55\textwidth} +\begin{lstlisting} +boxplot(Reaction ~ Days, sleepstudy) +\end{lstlisting} + \end{column} + \end{columns} +\end{frame} + +\begin{frame}[fragile]{Visualization of individual data} + \begin{columns} + \begin{column}{.53\textwidth} + \includegraphics[scale=.5]{../figures/sleep_subjects} + \end{column} + \begin{column}{.55\textwidth} +\begin{lstlisting} +library("lattice") + +xyplot(Reaction ~ Days | Subject, + data = sleepstudy, + type = c("g", "b"), + xlab = "Days of sleep deprivation", + ylab = "Average reaction time (ms)", + aspect = "xy") +\end{lstlisting} + \end{column} + \end{columns} +\end{frame} + +\begin{frame}{Random intercept model} + \begin{columns} + \begin{column}{.5\textwidth} + \begin{itemize} +\item The random intercept model adds a random intercept for each subject +\[ + y_{ij} = \beta_0 + \beta_1\,Days_{ij} + \upsilon_{0i} + \varepsilon_i +\] +with $\upsilon_{0i} \overset{iid}{\sim} N(0, \sigma^2_{\upsilon})$, + $\varepsilon_{ij} \overset{iid}{\sim} N(0, \sigma^2_\varepsilon)$ + \item The slope is identical for each subject (and the population) + \end{itemize} + \vspace{2.2cm} + \end{column} + \begin{column}{.5\textwidth} + \includegraphics[scale=.5]{../figures/sleep_random_intercept} + \end{column} + \end{columns} +\end{frame} + + +\begin{frame}{Fixed and random effects for random intercept model} +% m1 <- lmer(Reaction ~ Days + (1 | Subject), sleepstudy) +% tab <- cbind( +% sleepstudy[, c("Subject", "Days", "Reaction")], +% yhat = fixef(m1)["(Intercept)"] + +% fixef(m1)["Days"] * sleepstudy[["Days"]] + +% ranef(m1)[["Subject"]][["(Intercept)"]][sleepstudy[["Subject"]]] + +% resid(m1), +% Intercept = fixef(m1)["(Intercept)"], +% Slope = fixef(m1)["Days"] * sleepstudy[["Days"]], +% SubInt = ranef(m1)[["Subject"]][["(Intercept)"]][sleepstudy[["Subject"]]], +% Resid = resid(m1) +% ) +% tab <- cbind( +% sleepstudy[, c("Subject", "Days", "Reaction")], +% Intercept = fixef(m1)["(Intercept)"] |> unname(), +% Slope = fixef(m1)["Days"] |> unname(), +% SubInt = ranef(m1)[["Subject"]][["(Intercept)"]][sleepstudy[["Subject"]]], +% Resid = resid(m1) +% ) +% print( +% xtable(tab[tab$Subject %in% c("308", "309", "310"), ], +% digits = c(NA, 0, 0, 0, 1, 1, 1, 1)), +% include.rownames = FALSE, math.style.negative = TRUE) +% ) +% latex table generated in R 4.5.2 by xtable 1.8-8 package +% Thu Mar 5 11:18:32 2026 +\centering\footnotesize +\begin{tabular}{ccccccc} + \hline + Subject & Days & Reaction & Fixed & & Random & \\ + \cline{4-5} + \cline{6-7} + & & & Intercept & Slope & SubInt & Resid \\ + \hline + 308 & 0 & 250 & 251.4 & 10.5 & 40.8 & $-$42.6 \\ + 308 & 1 & 259 & 251.4 & 10.5 & 40.8 & $-$44.0 \\ + 308 & 2 & 251 & 251.4 & 10.5 & 40.8 & $-$62.3 \\ + \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots\\ + 308 & 8 & 431 & 251.4 & 10.5 & 40.8 & 54.7 \\ + 308 & 9 & 466 & 251.4 & 10.5 & 40.8 & 80.0 \\ + \hline + 309 & 0 & 223 & 251.4 & 10.5 & $-$77.8 & 49.2 \\ + 309 & 1 & 205 & 251.4 & 10.5 & $-$77.8 & 21.2 \\ + 309 & 2 & 203 & 251.4 & 10.5 & $-$77.8 & 8.5 \\ + \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots\\ + \hline + 310 & 0 & 199 & 251.4 & 10.5 & $-$63.1 & 10.8 \\ + 310 & 1 & 194 & 251.4 & 10.5 & $-$63.1 & $-$4.4 \\ + 310 & 2 & 234 & 251.4 & 10.5 & $-$63.1 & 25.1 \\ + \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots\\ + \hline + & & & & & $\sigma^2_{\upsilon}$ & $\sigma^2_\varepsilon$\\ +\end{tabular} +\end{frame} + +\begin{frame}[fragile]{Random slope model} + \begin{columns} + \begin{column}{.5\textwidth} + \begin{itemize} +\item The random slope model adds a random intercept and a random slope for each + subject +\[ + y_{ij} = \beta_0 + \beta_1\,Days_{ij} + \upsilon_{0i} + + \upsilon_{1i}\,Days_{ij} + \varepsilon_{ij} +\] +with +\begin{align*} + \begin{pmatrix} \upsilon_{0i}\\ \upsilon_{1i} \end{pmatrix} &\overset{iid}{\sim} + N \left(\begin{pmatrix} 0\\ 0 \end{pmatrix}, \, \mathbf{\Sigma}_\upsilon = + \begin{pmatrix} + \sigma^2_{\upsilon_0} & \sigma_{\upsilon_0 \upsilon_1} \\ + \sigma_{\upsilon_0 \upsilon_1} & \sigma^2_{\upsilon_1} \\ + \end{pmatrix} \right) + \\ + \varepsilon_{ij} & \overset{iid}{\sim} N(0, \sigma^2_\varepsilon) +\end{align*} + \vspace{-.5cm} + \item Individual slopes for each subject + \end{itemize} + \end{column} + \begin{column}{.5\textwidth} + \includegraphics[scale=.5]{../figures/sleep_random_slope} + \end{column} + \end{columns} +\end{frame} + + +\begin{frame}{Fixed and random effects for random slope model} +% m2 <- lmer(Reaction ~ Days + (Days | Subject), sleepstudy) +% tab <- cbind( +% sleepstudy[, c("Subject", "Days", "Reaction")], +% yhat = fixef(m2)["(Intercept)"] + +% fixef(m2)["Days"] * sleepstudy[["Days"]] + +% ranef(m2)[["Subject"]][["(Intercept)"]][sleepstudy[["Subject"]]] + +% ranef(m2)[["Subject"]][["Days"]][sleepstudy[["Subject"]]] * +% sleepstudy[["Days"]] + +% resid(m2), +% Intercept = fixef(m2)["(Intercept)"], +% Slope = fixef(m2)["Days"] * sleepstudy[["Days"]], +% SubInt = ranef(m2)[["Subject"]][["(Intercept)"]][sleepstudy[["Subject"]]], +% SubSlp = ranef(m2)[["Subject"]][["Days"]][sleepstudy[["Subject"]]] * +% sleepstudy[["Days"]], +% Resid = resid(m2) +% ) +% tab <- cbind( +% sleepstudy[, c("Subject", "Days", "Reaction")], +% Intercept = fixef(m2)["(Intercept)"] |> unname(), +% Slope = fixef(m2)["Days"] |> unname(), +% SubInt = ranef(m2)[["Subject"]][["(Intercept)"]][sleepstudy[["Subject"]]], +% SubSlp = ranef(m2)[["Subject"]][["Days"]][sleepstudy[["Subject"]]], +% Resid = resid(m2) +% ) +% print( +% xtable(tab[tab$Subject %in% c("308", "309", "310"), ], +% digits = c(NA, 0, 0, 0, 1, 1, 1, 1, 1)), +% include.rownames = FALSE, math.style.negative = TRUE) +% ) +% latex table generated in R 4.5.2 by xtable 1.8-8 package +% Thu Mar 5 11:18:32 2026 +\centering\footnotesize +\begin{tabular}{cccccccc} + \hline + Subject & Days & Reaction & Fixed & & Random & & \\ + \cline{4-5} + \cline{6-8} + & & & Intercept & Slope & SubInt & SubSlp & Resid \\ + \hline + 308 & 0 & 250 & 251.4 & 10.5 & 2.3 & 9.2 & $-$4.1 \\ + 308 & 1 & 259 & 251.4 & 10.5 & 2.3 & 9.2 & $-$14.6 \\ + 308 & 2 & 251 & 251.4 & 10.5 & 2.3 & 9.2 & $-$42.2 \\ + \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots \\ + 308 & 8 & 431 & 251.4 & 10.5 & 2.3 & 9.2 & 19.6 \\ + 308 & 9 & 466 & 251.4 & 10.5 & 2.3 & 9.2 & 35.7 \\ + \hline + 309 & 0 & 223 & 251.4 & 10.5 & $-$40.4 & $-$8.6 & 11.7 \\ + 309 & 1 & 205 & 251.4 & 10.5 & $-$40.4 & $-$8.6 & $-$7.6 \\ + 309 & 2 & 203 & 251.4 & 10.5 & $-$40.4 & $-$8.6 & $-$11.7 \\ + \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots \\ + \hline + 310 & 0 & 199 & 251.4 & 10.5 & $-$39.0 & $-$5.4 & $-$13.4 \\ + 310 & 1 & 194 & 251.4 & 10.5 & $-$39.0 & $-$5.4 & $-$23.1 \\ + 310 & 2 & 234 & 251.4 & 10.5 & $-$39.0 & $-$5.4 & 11.8 \\ + \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots \\ + \hline + &&&&& $\sigma^2_{\upsilon_0}$ & $\sigma^2_{\upsilon_1}$ & + $\sigma^2_\varepsilon$\\[-6pt] + &&&&& \multicolumn{2}{c}{$\sigma_{\upsilon_0\upsilon_1}$} & \\ +\end{tabular} +\end{frame} + +\begin{frame}{Partial pooling} + \begin{columns} + \begin{column}{.5\textwidth} + \includegraphics[scale=.5]{../figures/sleep_shrinkfit} + \end{column} + \begin{column}{.5\textwidth} + \begin{itemize} + \item Within-subject regression line shows regression line fitted to + data for each individual + \item Population regression line shows fixed effects for mixed-effects + model + \item Mixed model regression line shows individual regression lines as + predicted by mixed-effects models + \end{itemize} + \end{column} + \end{columns} +\end{frame} + +\section[Growth curve model]{Example: Depression and Imipramin} + + +\begin{frame}{Quadratic time trends} + \begin{columns} + \begin{column}{.35\textwidth} + \begin{itemize} + \item A lot of times the assumption of a linear time trend is too + simple + \item Change is not happening unbraked linearly but flattens out + \end{itemize} + \end{column} + \begin{column}{.65\textwidth} + \includegraphics[width=9cm]{../figures/quad} + \end{column} + \end{columns} +\end{frame} + +\begin{frame}{Quadratic time trends} + \begin{itemize} + \item Quadratic regression model +\begin{align*} + y_{ij} &= b_{0i} + b_{1i}\,t_{ij} + b_{2i}\,t^2_{ij} + \varepsilon_{ij}\\ + &= b_{0i} + (b_{1i} + b_{2i}\,t_{ij}) t_{ij} + \varepsilon_{ij} +\end{align*} + \item The linear change depends on time $t$ +\[ + \frac{\partial y}{\partial t} = b_{1i} + 2b_{2i} \, t +\] + \item The intercept $t = -b_{1i}/(2 b_{2i})$ is the point in time when a + positive (negative) trend becomes negative (postive) + \end{itemize} +\end{frame} + +\begin{frame}{Depression and Imipramin \citep{ReisbyGram77}} + \begin{itemize} + \item \citet{ReisbyGram77} studied the effect of Imipramin on 66 + inpatients treated for depression + \item Depression was measured with the Hamilton depression rating scale + (HDRS) + \item Additionally, the concentration of Imipramin and its metabolite + Desipramin was measured in their blood plasma + \item Patients were classified into endogenous and non-endogenous + depressed + \item Depression was measured weekly for 6 time points; the effect of + the antidepressant was observed starting at week 2 for four weeks + \end{itemize} + \vfill +\end{frame} + +\begin{frame}{Descriptive statistics} +\begin{columns} +\begin{column}{.55\textwidth} + \includegraphics[scale=.8]{../figures/hdrs-box} +\end{column} +\begin{column}{.6\textwidth} + HDRS score\\[1ex] + {\footnotesize +\begin{tabular}{rrrrrrr} + \hline + $t$ & W0 & W1 & W2 & W3 & W4 & W5 \\ + \hline + $M$ & 23.44 & 21.84 & 18.31 & 16.42 & 13.62 & 11.95 \\ + $SD$ & 4.53 & 4.70 & 5.49 & 6.42 & 6.97 & 7.22 \\ + $n$ & 61 & 63 & 65 & 65 & 63 & 58 \\ + \hline +\end{tabular} + } + + \vspace{.5cm} + Empirical correlation matrix of HDRS score\\[1ex] + {\footnotesize +\begin{tabular}{rrrrrrr} + \hline + & W0 & W1 & W2 & W3 & W4 & W5 \\ + \hline + Week 0 & 1 & .49 & .41 & .33 & .23 & .18 \\ + Week 1 & .49 & 1 & .49 & .41 & .31 & .22 \\ + Week 2 & .41 & .49 & 1 & .74 & .67 & .46 \\ + Week 3 & .33 & .41 & .74 & 1 & .82 & .57 \\ + Week 4 & .23 & .31 & .67 & .82 & 1 & .65 \\ + Week 5 & .18 & .22 & .46 & .57 & .65 & 1 \\ + \hline +\end{tabular} + } + + \vspace{1cm} +\end{column} +\end{columns} +\end{frame} + +\begin{frame}{Predictions random slope model} + \begin{columns} + \begin{column}{.6\textwidth} + \includegraphics[scale=.5]{../figures/hdrs-ind_pred-randomslope} + \end{column} + \begin{column}{.4\textwidth} + \[ + y_{ij} = \beta_0 + \beta_1 time + \upsilon_{0i} + \upsilon_{1i} time + + \varepsilon_{ij} + \] +with +\begin{align*} + \begin{pmatrix} \upsilon_{0i}\\ \upsilon_{1i} \end{pmatrix} &\overset{iid}\sim + N \left(\begin{pmatrix} 0\\ 0 \end{pmatrix}, \, \boldsymbol{\Sigma}_\upsilon = + \begin{pmatrix} + \sigma^2_{\upsilon_0} & \sigma_{\upsilon_0 \upsilon_1} \\ + \sigma_{\upsilon_0 \upsilon_1} & \sigma^2_{\upsilon_1} \\ + \end{pmatrix} \right)\\ + \boldsymbol{\varepsilon}_i &\overset{iid}\sim N(\mathbf{0}, \, \sigma^2 + \mathbf{I}_{n_i}) +\end{align*} + \end{column} +\end{columns} +\end{frame} + +\begin{frame}[fragile]{} + \begin{block}{Exercise} + \begin{itemize} + \item What would be a possible extension of this model?\pause + \item Write down the model equations + \begin{itemize} + \item What changes for the fixed effects? + \item How do the variance components for the random effects change? + \end{itemize}\pause + \item How can we interpret the random quadratic effects for this model?\pause + \item How do we add (random) quadratic effects to a random slope model + using \texttt{lme4::lmer()}?\pause + \item Fit this growth curve model in R\pause + \end{itemize} + \end{block} + \begin{lstlisting} +dat <- read.table("data/reisby.txt", header = TRUE) +dat$id <- factor(dat$id) +dat$diag <- factor(dat$diag, levels = c("nonen", "endog")) +dat <- na.omit(dat) # drop missing values + \end{lstlisting} +\end{frame} + +\begin{frame}[fragile]{Model with quadratic trend} + \begin{itemize} + \item Model with quadratic individual and quadratic group trend + \[ + y_{ij} = \beta_0 + \beta_1\,t_{ij} + \beta_2\,t^2_{ij} + \upsilon_{0i} + + \upsilon_1\,t_{ij} + \upsilon_2\,t^2_{ij} + \varepsilon_{ij} + \] +with +\begin{align*} + \begin{pmatrix} + \upsilon_{0i}\\ + \upsilon_{1i}\\ + \upsilon_{2i} + \end{pmatrix} &\overset{iid}{\sim} + N \left(\begin{pmatrix} + 0\\ 0\\ 0 + \end{pmatrix}, \, + \begin{pmatrix} + \sigma^2_{\upsilon_0} & \sigma_{\upsilon_0 \upsilon_1} & \sigma_{\upsilon_0 \upsilon_2}\\ + \sigma_{\upsilon_0 \upsilon_1} & \sigma^2_{\upsilon_1} & \sigma_{\upsilon_1 \upsilon_2}\\ + \sigma_{\upsilon_0 \upsilon_2} & \sigma_{\upsilon_1 \upsilon_2} & \sigma^2_{\upsilon_2}\\ + \end{pmatrix} \right) \\ + \boldsymbol{\varepsilon}_i &\overset{iid}{\sim} N(\mathbf{0}, \, \sigma^2 \mathbf{I}_{n_i}) +\end{align*} + \end{itemize} +\vspace{-.5cm} +\end{frame} + +\begin{frame}[fragile]{Model with quadratic trend} + \begin{lstlisting} +library("lme4") + +# random intercept model +lme1 <- lmer(hamd ~ week + (1 | id), data = dat, REML = FALSE) + +# random slope model +lme2 <- lmer(hamd ~ week + (week | id), data = dat, REML = FALSE) + +# model with quadratic time trend +lme3 <- lmer(hamd ~ week + I(week^2) + (week + I(week^2) | id), + data = dat, REML = FALSE) + \end{lstlisting} +\end{frame} + +\begin{frame}{Model predictions} +\begin{columns} +\begin{column}{6cm} +\includegraphics[width=6cm]{../figures/hdrs-quad} +\end{column} +% +\begin{column}{5cm} + \begin{itemize} + \item Averaged over persons an approximately linear trend is obtained, + $\hat{\beta}_1 = -2.63$, $\hat{\beta}_2 = 0.05$ + \item Some of the predicted individual trends are strongly nonlinear + \end{itemize} +\end{column} +\end{columns} + \begin{itemize} + \item Test against a model without individual quadratic trends\\[2ex] + +H$_0$: $\sigma^2_{\upsilon_2} = \sigma_{\upsilon_0 \upsilon_2} = +\sigma_{\upsilon_1 \upsilon_2} = 0$ \qquad +$G^2(3) = 10.98$, $p = .012$ + \end{itemize} +\end{frame} + +\begin{frame}[fragile]{Model predictions} + \begin{columns} + \begin{column}{8cm} + \includegraphics[scale=.5]{../figures/hdrs-ind_pred-quad} + \end{column} + \begin{column}{6cm} + \begin{lstlisting} +xyplot( + hamd + predict(lme3) + ~ week | id, + data = dat, + type = c("p", "l", "g"), + pch = 16, + distribute.type = TRUE, + ylab = "HDRS score", + xlab = "Time (Week)") + \end{lstlisting} + \end{column} + \end{columns} +\end{frame} + +% \begin{frame}[fragile]{Implied marginal covariance matrix} +% Predicted +% \begin{align*} +% \mathbf{Z}_i \boldsymbol{\hat\Sigma}_\upsilon \mathbf{Z}'_i + +% \hat\sigma^2 \mathbf{I}_{n_i} &= +% \begin{pmatrix} +% 20.96 & 9.41 & 8.16 & 6.68 & 4.98 & 3.06 \\ +% 9.41 & 23.86 & 15.57 & 16.08 & 14.88 & 11.97 \\ +% 8.16 & 15.57 & 31.07 & 23.11 & 23.26 & 20.98 \\ +% 6.68 & 16.08 & 23.11 & 38.31 & 30.12 & 30.09 \\ +% 4.98 & 14.88 & 23.26 & 30.12 & 45.98 & 39.29 \\ +% 3.06 & 11.97 & 20.98 & 30.09 & 39.29 & 59.11 +% \end{pmatrix}\\ +% % +% \intertext{Observed} +% % +% \widehat{Cov}(\mathbf{y}_i) &= +% \begin{pmatrix} +% 20.55 & 10.11 & 10.14 & 10.09 & 7.19 & 6.28 \\ +% 10.11 & 22.07 & 12.28 & 12.55 & 10.26 & 7.72 \\ +% 10.14 & 12.28 & 30.09 & 25.13 & 24.63 & 18.38 \\ +% 10.09 & 12.55 & 25.13 & 41.15 & 37.34 & 23.99 \\ +% 7.19 & 10.26 & 24.63 & 37.34 & 48.59 & 30.51 \\ +% 6.28 & 7.72 & 18.38 & 23.99 & 30.51 & 52.12 \\ +% \end{pmatrix}\\ +% \end{align*} +% % \vspace{-1cm} +% % \begin{lstlisting} +% % getVarCov(lme3, type="marginal") # predicted cov. matrix +% % +% % var.d <- crossprod(getME(lme3,"Lambdat")) +% % Zt <- getME(lme3, "Zt") +% % vr <- sigma(lme3)^2 +% % +% % var.b <- vr*(t(Zt) %*% var.d %*% Zt) +% % sI <- vr * Diagonal(nrow(dat2)) +% % var.y <- var.b + sI +% % +% % \end{lstlisting} +% \end{frame} + +\begin{frame}{Centering variables} + \begin{itemize} + \item If multiples of the time variables ($t$, $t^2$, $t^3$, etc.) are + entered into the regression equation, multicollinearity can become a + problem + \item For example, $t = 0, 1, 2, 3$ and $t^2 = 0, 1, 4, 9$ correlate almost + perfectly + \item By centering the variables, this problem can be diminished: $(t - + \bar{t}) = -1.5, -0.5, 0.5, 1.5$ and $(t - \bar{t})^2 = 2.25, 0.25, 0.25, + 2.25$ are uncorrelated + \item By centering variables the interpretation of the intercept in a linear + model changes: + \begin{itemize} + \item Uncentered intercepts represent the difference to the first time + point ($t = 0$) + \item Centered intercepts represent the difference after half of the + time + \end{itemize} + \end{itemize} +\end{frame} + +\begin{frame}{} + \begin{block}{Exercise} + \begin{itemize} + \item Center time\pause + \item Compare the correlation between time and time$^2$ with and without + centering\pause + \item Refit the model with quadratic effects when time is centered\pause + \item Which parameters change? + \begin{itemize} + \item What changes for the fixed effects? + \item How do the variance components for the random effects change? + \end{itemize}\pause + \item \emph{Why} do the variance components for the random effects change? + \end{itemize} + \end{block} +\end{frame} + +\begin{frame}[fragile]{Analysis with centered time variable} +\begin{lstlisting} +dat$week_c <- dat$week - mean(dat$week) +cor(dat$week, dat$week^2) # 0.96 +cor(dat$week_c, dat$week_c^2) # 0.01 + +# random slope model +lme2c <- lmer(hamd ~ week_c + (week_c | id), data = dat, REML = FALSE) + +# model with quadratic time trend +lme3c <- lmer(hamd ~ week_c + I(week_c^2) + (week_c + I(week_c^2)|id), + data = dat, REML = FALSE) +\end{lstlisting} +% \begin{itemize} +% \item When comparing the estimated parameters, it becomes obvious that not +% only the intercept changes but the estimates for the (co)variances do as +% well +% \item Why?\pause +% ~Be sure to make an informed choice when centering your +% variables! +% \end{itemize} + \nocite{Alday2025, Hedeker2006} +\end{frame} + +\begin{frame}[fragile]{Investigating random effects structure} + \begin{itemize} + \item In order to get a better understanding of the necessary random effects + it might be a good idea to take a closer look at them + \item Two plots often used are the so-called caterpillar and shrinkage plots + \item Play around with different models and compare how, e.\,g., the + caterpillar plots change with and without covariances in the model! + \end{itemize} + \vfill +\begin{lstlisting} +library("lattice") +dotplot(ranef(lme3), scales = list( x = list(relation = "free")))$id +\end{lstlisting} +\end{frame} + + +\begin{frame}[fragile]%{Investigating random effects structure} + {Caterpillar plot} + \begin{columns} + \begin{column}{1.13\textwidth} + \includegraphics[scale=.31]{../figures/hdrs-caterpillar} + \end{column} + \end{columns} +\end{frame} + +\begin{frame}[fragile]%{Investigating random effects structure} + {Shrinkage plots} + \begin{columns} + \begin{column}{.35\textwidth} + \includegraphics[scale=.35]{../figures/hdrs_shrinkage_int-week} + \end{column} + \begin{column}{.35\textwidth} + \includegraphics[scale=.35]{../figures/hdrs_shrinkage_int-weeksq} + \end{column} + \begin{column}{.35\textwidth} + \includegraphics[scale=.35]{../figures/hdrs_shrinkage_week-weeksq} + \end{column} + \end{columns} +\end{frame} + +% \begin{frame}{Higher-order polynomials} +% \begin{itemize} +% \item Nonlinear time trends can be modelled in a flexibel and parsimonious +% way by using higher-order polynomials +% \item For example, saddle or reversal points in a time trend can be +% described +% \item Polynomials have the advantage that the regression model stays linear +% in its parameters +% \item They have the disadvantage that extrapolated values can quickly be +% outside of a range that can still be interpreted in a meaningful way +% \end{itemize} +% \end{frame} +% +% \begin{frame}{Cubic time trends} +% \begin{center} +% \includegraphics[width=9cm]{../figures/cubic} +% \end{center} +% \end{frame} +% +% \begin{frame}{Polynomial regression: Extrapolation} +% \begin{center} +% \includegraphics[width=9cm]{../figures/cubic-gone-bad} +% \end{center} +% \end{frame} + +\begin{frame}{What we learned today\dots} + {\dots and how to go on} + \pause + \begin{enumerate}[<+->] + \item We learned + \begin{itemize} + \item The basic concept of random effects and why to include them in a + model + \item How to model data collected over several time points + \item How to compute mixed-effects model with quadratic time trends in R + using \texttt{lmer()} from the lme4 package + \item How to interpret parameters in mixed-effects model with quadratic + effects + \end{itemize} + \item Next steps + \begin{itemize} + \item Do this exercise + \url{https://gitea.iwm-tuebingen.de/nwickelmaier/lead_longitudinal/src/branch/master/exercises/schizo.md} + \item It has a very similar structure than the depression dataset and + this will help you to generalize the concepts we learned today + \item You can send questions to me and even make an appointment with me + to go over your solution + \end{itemize} + \end{enumerate} +\end{frame} + +\appendix +%\begin{frame}[allowframebreaks]{References} +\begin{frame}{References} + \printbibliography + \vfill +\end{frame} + +\end{document} + diff --git a/slides/lit.bib b/slides/lit.bib new file mode 100644 index 0000000..7522d6d --- /dev/null +++ b/slides/lit.bib @@ -0,0 +1,55 @@ +@book{Alday2025, + title = {Embrace Uncertainty -- {M}ixed-effects models with {J}ulia}, + author = {Phillip Alday and Reinhold Kliegl and Douglas Bates}, + year = {2025}, + url = {https://embraceuncertaintybook.com/} +} + +@misc{Bates10lme4, + title = {{lme4}: Mixed-effects modeling with {R} (book draft)}, + author = {Bates, D.}, + year = {2010}, + url = {https://lme4.r-forge.r-project.org/} +} + +@article{Bates15, + title = {Fitting Linear Mixed-Effects Models Using lme4}, + volume = {67}, + url = {https://CRAN.R-project.org/package=lme4/vignettes/lmer.pdf}, + doi = {10.18637/jss.v067.i01}, + number = {1}, + journal = {Journal of Statistical Software}, + author = {Bates, D. and M\"achler, M. and Bolker, B. and Walker, S.}, + year = {2015}, + pages = {1--48} +} + +@book{Hedeker2006, + author = {Hedeker, D. R. and Gibbons, R. D.}, + title = {Longitudinal data analysis}, + year = {2006}, + address = {Hoboken}, + publisher = {John Wiley} +} + +@article{Moeller08, + title = {Short-term treatment with risperidone or haloperidol in first-episode schizophrenia: 8-week results of a randomized controlled trial within the German Research Network on Schizophrenia}, + author = {M{\"o}ller, H.-J. and Riedel, M. and J{\"a}ger, M. and Wickelmaier, F. and Maier, W. and K{\"u}hn, K.-U. and Buchkremer, G. and Heuser, I. and Klosterk{\"o}tter, J. and Gastpar, M. and others}, + journal = {International Journal of Neuropsychopharmacology}, + volume = {11}, + number = {7}, + pages = {985--997}, + year = {2008}, + doi = {10.1017/S1461145708008791} +} + +@article{ReisbyGram77, + author = {Reisby, N. and Gram, L. F. and Bech, P. and Nagy, A. and Petersen, G. O. and Ortmann, J. and Ibsen, I. and Dencker, S. J. and Jacobsen, O. and Krautwald, O. and Sondergaard, I. and Christiansen, J.}, + title = {Imipramine: Clinical effects and pharmacokinetic variability}, + journal = {Psychopharmacology}, + year = {1977}, + volume = {54}, + pages = {263--272}, + doi = {10.1007/BF00426574} +} +