txt> print("Descriptive Statistics from Occasion-Level-1 Data")
[1] "Descriptive Statistics from Occasion-Level-1 Data"
txt> univSummary(data = Example5b, vars = c("baseage", "time", "dep", 
txt+ "PMgame", "WPgame"))

 Variable   N  Mean    SD Variance    Min    Max Median
 baseage  800 9.899 0.646    0.418  8.200 12.120  9.855
 time     800 1.500 1.121    1.258  0.000  3.310  1.540
 dep      800 9.581 3.611   13.042  0.000 20.000  9.000
 PMgame   800 3.104 1.061    1.125  0.750  6.250  3.000
 WPgame   800 0.000 0.869    0.755 -2.500  2.500  0.000

txt> LongPlot = rbind(data.frame(PersonID = Example5b$PersonID, time = Example5b$time, 
txt+ roundtime = Example5b$roundtime, Variable = "Depression", 
txt+ value = Example5b$dep), data.frame(PersonID = Example5b$PersonID, 
txt+ time = Example5b$time, roundtime = Example5b$roundtime, Variable = "Hours of Daily Gaming", 
txt+ value = Example5b$game))
txt> LongPlot$Variable = factor(LongPlot$Variable, levels = c("Depression", 
txt+ "Hours of Daily Gaming"))
txt> FigureSpaghetti = ggplot(LongPlot, aes(x = time, y = value)) + 
txt+ geom_line(aes(group = PersonID, color = factor(PersonID)), 
txt+ 
txt+ linewidth = 0.4, alpha = 0.7) + stat_summary(aes(x = roundtime, 
txt+ group = 1), fun = mean, geom = "line", color = "black", linewidth = 1.5) + 
txt+ stat_summary(aes(x = roundtime, group = 1), fun = mean, geom = "point", 
txt+ 
txt+ color = "black", size = 3) + facet_wrap(~Variable, nrow = 1, 
txt+ scales = "free_y") + guides(color = "none") + scale_x_continuous(breaks = 0:3, 
txt+ expand = expansion(mult = 0.05)) + labs(x = "Years Since Baseline", 
txt+ y = "Observed Value", title = expression(bold("Individual Trajectories of Depression and Gaming")), 
txt+ subtitle = expression(bolditalic("(thin colored lines = children; thick black line = mean at each occasion)"))) + 
txt+ theme_classic(base_size = 14) + theme(plot.title = element_text(face = "bold", 
txt+ hjust = 0.5, size = 16), plot.subtitle = element_text(face = "bold.italic", 
txt+ hjust = 0.5, size = 14), axis.title = element_text(face = "bold", 
txt+ size = 12), axis.text = element_text(size = 12, color = "black"), 
txt+ strip.text = element_text(face = "bold", size = 13), strip.background = element_rect(fill = "gray90", 
txt+ 
txt+ color = "black"), panel.border = element_rect(color = "black", 
txt+ 
txt+ fill = NA, linewidth = 0.8), panel.spacing = unit(1, 
txt+ 
txt+ "lines"))
txt> ggsave("Example5b_Spaghetti.png", FigureSpaghetti, width = 10, 
txt+ height = 5, dpi = 300)
txt> print(FigureSpaghetti)
txt> print("Model 1: Empty Means, Random Intercept for Person, to Get ICC for Depression")
[1] "Model 1: Empty Means, Random Intercept for Person, to Get ICC for Depression"
txt> Model1 = lmer(data = Example5b, REML = TRUE, formula = dep ~ 
txt+ 1 + (1 | PersonID))
txt> lmerSummary(Model1, ICC1 = TRUE, LRT = TRUE)

Model Fit for Predicting dep:
      AIC       BIC    logLik -2*logLik 
 4143.157  4157.211 -2068.579  4137.157 

Model for the Variance:
 Groups    Name      Variance Std.Dev
 PersonID  Intercept    5.905   2.430
 Residual               7.160   2.676
Number of obs: 800 , groups: PersonID, 200 

Fixed Effects:
 Term         Est    SE      DF t.value p.value
 Intercept  9.581 0.196 199.000  48.847   <.001

ICC1 intraclass correlation:
    Adjusted ICC: 0.452 
  Unadjusted ICC: 0.452 

LRT for random effect(s):
 Random Effect Dropped  Npar    LL -2LL  AIC  BIC     LRT DF p.value p.mix
 <none>                    3 -2069 4137 4143 4157                         
 Intercept | PersonID      2 -2163 4326 4330 4340 188.975  1   <.001 <.001

txt> print("Model 2: Saturated Means, UN Variances and Correlations by Rounded Occasion")
[1] "Model 2: Saturated Means, UN Variances and Correlations by Rounded Occasion"
txt> Model2 = gls(data = Example5b, method = "REML", model = dep ~ 
txt+ 1 + factor(occasion), correlation = corSymm(form = ~occasion | 
txt+ PersonID), weights = varIdent(form = ~1 | occasion))
txt> glsSummary(Model2, predMeans = TRUE)

Model Fit for Predicting dep:
      AIC      BIC    logLik -2*logLik
 3996.206 4061.721 -1984.103  3968.206

Residual Variance (sigma^2):
11.220 

Residual Correlation (corSymm):
 Structure Parameter Estimate
 Symm      cor(1,2)     0.606
 Symm      cor(1,3)     0.497
 Symm      cor(1,4)     0.409
 Symm      cor(2,3)     0.553
 Symm      cor(2,4)     0.486
 Symm      cor(3,4)     0.616

Residual Variance by Group (varIdent; SD.ratio is relative to the first group,
which is the single Residual Variance above):
 Group SD.ratio Residual.SD Residual.Var
 1        1.000       3.350       11.220
 2        0.971       3.254       10.587
 3        1.009       3.380       11.425
 4        1.143       3.830       14.668

Overall F-test for All Fixed Slopes:
 Model  NumDF   DenDF F.value p.value
 Model2     3 198.858  35.875   <.001

Fixed Effects:
 Term                 Est    SE      DF t.value p.value
 Intercept          8.225 0.237 206.548  34.725   <.001
 factor(occasion)2  0.855 0.207 199.035   4.126   <.001
 factor(occasion)3  1.725 0.239 199.128   7.230   <.001
 factor(occasion)4  2.845 0.277 198.292  10.257   <.001
Denominator df method: Satterthwaite (each effect has its own df)

Model-Predicted Means for factor(occasion):
 Level   Mean    SE      DF
 1      8.225 0.237 206.548
 2      9.080 0.230 194.243
 3      9.950 0.239 199.421
 4     11.070 0.271 200.661

txt> print("R matrix (variances and covariances) and correlations for first person")
[1] "R matrix (variances and covariances) and correlations for first person"
txt> getVarCov(Model2, individual = "1", type = "marginal")
Marginal variance covariance matrix
        [,1]    [,2]    [,3]    [,4]
[1,] 11.2200  6.6101  5.6294  5.2505
[2,]  6.6101 10.5870  6.0844  6.0597
[3,]  5.6294  6.0844 11.4250  7.9733
[4,]  5.2505  6.0597  7.9733 14.6680
  Standard Deviations: 3.3497 3.2537 3.38 3.8299 
txt> corMatrix(Model2$modelStruct$corStruct)[[1]]
          [,1]      [,2]      [,3]      [,4]
[1,] 1.0000000 0.6064888 0.4972048 0.4092631
[2,] 0.6064888 1.0000000 0.5532491 0.4862755
[3,] 0.4972048 0.5532491 1.0000000 0.6159253
[4,] 0.4092631 0.4862755 0.6159253 1.0000000
txt> print("Model 3a: Random Linear Change over Exact Time")
[1] "Model 3a: Random Linear Change over Exact Time"
txt> Model3a = lmer(data = Example5b, REML = TRUE, formula = dep ~ 
txt+ 1 + time + (1 + time | PersonID))
txt> lmerSummary(Model3a, LRT = TRUE, randCI = TRUE)

Model Fit for Predicting dep:
      AIC       BIC    logLik -2*logLik 
 3976.546  4004.653 -1982.273  3964.546 

Model for the Variance:
 Groups    Name      Variance Std.Dev   Corr
 PersonID  Intercept    7.045   2.654       
           time         0.793   0.890 -0.321
 Residual               4.341   2.084       
Number of obs: 800 , groups: PersonID, 200 

Overall F-test for All Fixed Slopes:
 Model   NumDF   DenDF F.value p.value
 Model3a     1 199.507 108.269   <.001

Fixed Effects:
 Term         Est    SE      DF t.value p.value
 Intercept  8.156 0.225 198.742  36.327   <.001
 time       0.948 0.091 199.507  10.405   <.001

LRT for random effect(s):
 Random Effect Dropped  Npar    LL -2LL  AIC  BIC    LRT DF p.value p.mix
 <none>                    6 -1982 3965 3977 4005                        
 time | PersonID           4 -1999 3998 4006 4025 33.241  2   <.001 <.001

Random Effect 95% CIs (Fixed Effect +/- 1.96*SQRT(random effect variance)):
 Term                 FixedEffect Variance    SD  Lower  Upper
 Intercept | PersonID       8.156    7.045 2.654  2.953 13.358
 time | PersonID            0.948    0.793 0.890 -0.797  2.693

txt> MeanM2 = predict(Model2, newdata = data.frame(occasion = 1:4))
txt> MeanM3a = predict(Model3a, newdata = data.frame(time = 0:3), 
txt+ re.form = NA)
txt> VarM2 = diag(as.matrix(getVarCov(Model2, individual = "1", type = "marginal")))
txt> V3a = marginalV(Model3a, times = 0:3)

Model-Implied Marginal Variance-Covariance Matrix (V) for Model3a:
       0      1      2      3
0 11.386  6.287  5.530  4.772
1  6.287 10.664  6.358  6.393
2  5.530  6.358 11.527  8.014
3  4.772  6.393  8.014 13.976

Model-Implied Marginal Correlation Matrix:
      0     1     2     3
0 1.000 0.571 0.483 0.378
1 0.571 1.000 0.573 0.524
2 0.483 0.573 1.000 0.631
3 0.378 0.524 0.631 1.000

(Rows and columns are values of time)

txt> VarM3a = diag(V3a$V)
txt> print("Predicted means and variances: Model 2 versus Model 3a")
[1] "Predicted means and variances: Model 2 versus Model 3a"
txt> print(data.frame(time = 0:3, MeanModel2 = MeanM2, MeanModel3a = MeanM3a, 
txt+ VarModel2 = VarM2, VarModel3a = VarM3a, row.names = NULL), 
txt+ digits = 5)
  time MeanModel2 MeanModel3a VarModel2 VarModel3a
1    0      8.225      8.1555    11.220     11.386
2    1      9.080      9.1034    10.587     10.664
3    2      9.950     10.0514    11.425     11.527
4    3     11.070     10.9993    14.668     13.976
txt> ModelLabels = c("Model 2 (Saturated Means, UN Variances)", "Model 3a (Random Linear Change)")
txt> PlotModels = data.frame(time = 0:3, Model = factor(rep(ModelLabels, 
txt+ each = 4), levels = ModelLabels), mean = unname(c(MeanM2, 
txt+ MeanM3a)), variance = unname(c(VarM2, VarM3a)))
txt> PlotStyle = list(scale_color_manual(values = c("black", "black")), 
txt+ scale_linetype_manual(values = c("solid", "dashed")), scale_x_continuous(breaks = 0:3, 
txt+ 
txt+ expand = expansion(mult = 0.05)), theme_classic(base_size = 14), 
txt+ theme(plot.title = element_text(face = "bold", hjust = 0.5, 
txt+ 
txt+ size = 16), axis.title = element_text(face = "bold", 
txt+ 
txt+ size = 12), axis.text = element_text(size = 12, color = "black"), 
txt+ 
txt+ panel.grid.major.y = element_line(color = "black", linewidth = 0.4), 
txt+ 
txt+ panel.border = element_rect(color = "black", fill = NA, 
txt+ 
txt+ 
txt+ linewidth = 0.8), legend.position = "top", legend.title = element_blank(), 
txt+ 
txt+ legend.key.width = unit(1.6, "cm"), legend.text = element_text(size = 11)))
txt> FigureMeans = ggplot(PlotModels, aes(x = time, y = mean, color = Model, 
txt+ linetype = Model)) + geom_line(linewidth = 1.2) + geom_point(data = PlotModels[1:4, 
txt+ ], size = 3.5) + PlotStyle + labs(x = "Years Since Baseline", 
txt+ y = "Predicted Mean Depression", title = "Predicted Means: Model 2 versus Model 3a Fixed Effects")
txt> ggsave("Example5b_Means.png", FigureMeans, width = 8.5, height = 5, 
txt+ dpi = 300)
txt> print(FigureMeans)
txt> FigureVar = ggplot(PlotModels, aes(x = time, y = variance, color = Model, 
txt+ linetype = Model)) + geom_line(linewidth = 1.2) + geom_point(data = PlotModels[1:4, 
txt+ ], size = 3.5) + PlotStyle + labs(x = "Years Since Baseline", 
txt+ y = "Predicted Variance of Depression", title = "Predicted Variances: Model 2 versus Model 3a")
txt> ggsave("Example5b_Variances.png", FigureVar, width = 8.5, height = 5, 
txt+ dpi = 300)
txt> print(FigureVar)
txt> print("Model 3b: Random Quadratic Change over Exact Time")
[1] "Model 3b: Random Quadratic Change over Exact Time"
txt> Model3b = lmer(data = Example5b, REML = TRUE, formula = dep ~ 
txt+ 1 + time + I(time^2) + (1 + time + I(time^2) | PersonID))
txt> lmerSummary(Model3b, LRT = TRUE, randCI = TRUE)

Model Fit for Predicting dep:
      AIC       BIC    logLik -2*logLik 
 3983.686  4030.533 -1981.843  3963.686 

Model for the Variance:
 Groups    Name      Variance Std.Dev   Corr       
 PersonID  Intercept    7.321   2.706              
           time         1.397   1.182 -0.328       
           I(time^2)    0.143   0.378  0.107 -0.687
 Residual               4.061   2.015              
Number of obs: 800 , groups: PersonID, 200 

Overall F-test for All Fixed Slopes:
 Model   NumDF   DenDF F.value p.value
 Model3b     2 200.950  54.172   <.001

Fixed Effects:
 Term         Est    SE      DF t.value p.value
 Intercept  8.233 0.236 198.946  34.846   <.001
 time       0.716 0.238 202.832   3.005   0.003
 I(time^2)  0.077 0.076 202.959   1.014   0.312

LRT for random effect(s):
 Random Effect Dropped  Npar    LL -2LL  AIC  BIC   LRT DF p.value p.mix
 <none>                   10 -1982 3964 3984 4031                       
 time | PersonID           7 -1983 3966 3980 4013 2.416  3   0.491 0.395
 I(time^2) | PersonID      7 -1983 3967 3981 4014 3.077  3   0.380 0.297

Random Effect 95% CIs (Fixed Effect +/- 1.96*SQRT(random effect variance)):
 Term                 FixedEffect Variance    SD  Lower  Upper
 Intercept | PersonID       8.233    7.321 2.706  2.930 13.536
 time | PersonID            0.716    1.397 1.182 -1.601  3.032
 I(time^2) | PersonID       0.077    0.143 0.378 -0.664  0.819

txt> print("Model 3b: Random Quadratic Change over Exact Time using lme")
[1] "Model 3b: Random Quadratic Change over Exact Time using lme"
txt> Model3b_lme = lme(data = Example5b, method = "REML", fixed = dep ~ 
txt+ 1 + time + I(time^2), random = ~1 + time + I(time^2) | PersonID)
txt> lmeSummary(Model3b_lme, LRT = TRUE, randCI = TRUE)

Model Fit for Predicting dep:
      AIC       BIC    logLik -2*logLik 
 3983.686  4030.495 -1981.843  3963.686 

Model for the Variance:
 Groups    Name      Variance Std.Dev   Corr       
 PersonID  Intercept    7.321   2.706              
           time         1.396   1.182 -0.328       
           I(time^2)    0.143   0.378  0.107 -0.687
 Residual               4.061   2.015              
Number of obs: 800 , groups: PersonID , 200 

Overall F-test for All Fixed Slopes:
 Model       NumDF   DenDF F.value p.value
 Model3b_lme     2 222.665  54.177   <.001

Fixed Effects:
 Term         Est    SE      DF t.value p.value
 Intercept  8.233 0.236 209.137  34.846   <.001
 time       0.716 0.238 246.167   3.005   0.003
 I(time^2)  0.077 0.076 236.200   1.014   0.312

LRT for each random effect and residual correlation, tested separately (full model vs. that piece dropped):
 Random Effect Dropped        Npar    LL -2LL  AIC  BIC    LRT DF p.value p.mix
 <none>                         10 -1982 3964 3984 4031                        
 time | PersonID                 7 -1983 3966 3980 4013  2.416  3   0.491 0.395
 I(time^2) | PersonID            7 -1983 3967 3981 4014  3.077  3   0.380 0.297
 time + I(time^2) | PersonID     5 -2000 4000 4010 4033 36.173  5   <.001      

Random Effect 95% CIs (Fixed Effect +/- 1.96*SQRT(random effect variance)):
 Term                 FixedEffect Variance    SD  Lower  Upper
 Intercept | PersonID       8.233    7.321 2.706  2.930 13.537
 time | PersonID            0.716    1.396 1.182 -1.600  3.032
 I(time^2) | PersonID       0.077    0.143 0.378 -0.664  0.819

txt> print("Model 3c: Add AR1 Residual Correlation to Random Linear Change")
[1] "Model 3c: Add AR1 Residual Correlation to Random Linear Change"
txt> Model3c = lme(data = Example5b, method = "REML", fixed = dep ~ 
txt+ 1 + time, random = ~1 + time | PersonID, correlation = corCAR1(form = ~time | 
txt+ PersonID))
txt> lmeSummary(Model3c, LRT = TRUE)

Model Fit for Predicting dep:
      AIC       BIC    logLik -2*logLik 
 3976.694  4009.469 -1981.347  3962.694 

Model for the Variance:
 Groups    Name      Variance Std.Dev   Corr
 PersonID  Intercept    6.345   2.519       
           time         0.632   0.795 -0.256
 Residual               4.883   2.210       
 Residual  CAR1 Phi                    0.104
Number of obs: 800 , groups: PersonID , 200 

Overall F-test for All Fixed Slopes:
 Model   NumDF   DenDF F.value p.value
 Model3c     1 202.087 108.640   <.001

Fixed Effects:
 Term         Est    SE      DF t.value p.value
 Intercept  8.163 0.224 201.943  36.391   <.001
 time       0.947 0.091 202.087  10.423   <.001

LRT for each random effect and residual correlation, tested separately (full model vs. that piece dropped):
 Random Effect Dropped        Npar    LL -2LL  AIC  BIC   LRT DF p.value p.mix
 <none>                          7 -1981 3963 3977 4009                       
 time | PersonID                 5 -1986 3972 3982 4005 9.052  2   0.011 0.007
 Residual correlation (CAR1)     6 -1982 3965 3977 4005 1.852  1   0.174      

txt> print("Model 4: Add Moderation of Linear Change by Baseline Age (centered at 10)")
[1] "Model 4: Add Moderation of Linear Change by Baseline Age (centered at 10)"
txt> Model4 = lmer(data = Example5b, REML = TRUE, formula = dep ~ 
txt+ 1 + time * baseage10 + (1 + time | PersonID))
txt> lmerSummary(Model4)

Model Fit for Predicting dep:
      AIC       BIC    logLik -2*logLik 
 3918.907  3956.384 -1951.453  3902.907 

Model for the Variance:
 Groups    Name      Variance Std.Dev   Corr
 PersonID  Intercept    6.270   2.504       
           time         0.342   0.585 -0.108
 Residual               4.341   2.083       
Number of obs: 800 , groups: PersonID, 200 

Overall F-test for All Fixed Slopes:
 Model  NumDF   DenDF F.value p.value
 Model4     3 199.465  74.976   <.001

Fixed Effects:
 Term               Est    SE      DF t.value p.value
 Intercept        8.018 0.218 197.359  36.727   <.001
 time             1.050 0.079 198.277  13.362   <.001
 baseage10       -1.400 0.334 196.782  -4.197   <.001
 time:baseage10   1.046 0.120 198.979   8.697   <.001

txt> PseudoR2(baseModel = Model3a, altModel1 = Model4)

Pseudo-R2
 Term       Model   Variance Pseudo.R2
 Intercept  Model3a    7.045          
 Intercept  Model4     6.270     0.110
 time       Model3a    0.793          
 time       Model4     0.342     0.568
 Residual   Model3a    4.341          
 Residual   Model4     4.341     0.000

txt> TotalR2(altModel1 = Model3a, altModel2 = Model4)

Total-R2 and Change in Total-R2
 Model   Total.R2 Change.R2
 Model3a    0.089          
 Model4     0.136     0.047

txt> print("Empty Means, Random Intercept to Get ICC for the Gaming Predictor")
[1] "Empty Means, Random Intercept to Get ICC for the Gaming Predictor"
txt> Gaming1 = lmer(data = Example5b, REML = TRUE, formula = game ~ 
txt+ 1 + (1 | PersonID))
txt> lmerSummary(Gaming1, ICC1 = TRUE, LRT = TRUE)

Model Fit for Predicting game:
      AIC       BIC    logLik -2*logLik 
 2583.508  2597.562 -1288.754  2577.508 

Model for the Variance:
 Groups    Name      Variance Std.Dev
 PersonID  Intercept    0.878   0.937
 Residual               1.005   1.003
Number of obs: 800 , groups: PersonID, 200 

Fixed Effects:
 Term         Est    SE      DF t.value p.value
 Intercept  3.104 0.075 199.000  41.300   <.001

ICC1 intraclass correlation:
    Adjusted ICC: 0.466 
  Unadjusted ICC: 0.466 

LRT for random effect(s):
 Random Effect Dropped  Npar    LL -2LL  AIC  BIC     LRT DF p.value p.mix
 <none>                    3 -1289 2578 2584 2598                         
 Intercept | PersonID      2 -1389 2779 2783 2792 201.169  1   <.001 <.001

txt> print("Gaming Model 2: Saturated Means, UN Variances and Correlations by Rounded Occasion")
[1] "Gaming Model 2: Saturated Means, UN Variances and Correlations by Rounded Occasion"
txt> Gaming2 = gls(data = Example5b, method = "REML", model = game ~ 
txt+ 1 + factor(occasion), correlation = corSymm(form = ~occasion | 
txt+ PersonID), weights = varIdent(form = ~1 | occasion))
txt> glsSummary(Gaming2, predMeans = TRUE)

Model Fit for Predicting game:
      AIC      BIC    logLik -2*logLik
 2597.187 2662.701 -1284.594  2569.187

Residual Variance (sigma^2):
1.751 

Residual Correlation (corSymm):
 Structure Parameter Estimate
 Symm      cor(1,2)     0.478
 Symm      cor(1,3)     0.410
 Symm      cor(1,4)     0.548
 Symm      cor(2,3)     0.512
 Symm      cor(2,4)     0.421
 Symm      cor(3,4)     0.441

Residual Variance by Group (varIdent; SD.ratio is relative to the first group,
which is the single Residual Variance above):
 Group SD.ratio Residual.SD Residual.Var
 1        1.000       1.323        1.751
 2        1.079       1.427        2.037
 3        1.035       1.369        1.874
 4        1.031       1.363        1.859

Overall F-test for All Fixed Slopes:
 Model   NumDF   DenDF F.value p.value
 Gaming2     3 200.430   1.784   0.151

Fixed Effects:
 Term                  Est    SE      DF t.value p.value
 Intercept           3.115 0.094 223.454  33.296   <.001
 factor(occasion)2  -0.055 0.100 199.951  -0.553   0.581
 factor(occasion)3  -0.110 0.103 198.878  -1.064   0.289
 factor(occasion)4   0.120 0.090 199.512   1.328   0.186
Denominator df method: Satterthwaite (each effect has its own df)

Model-Predicted Means for factor(occasion):
 Level  Mean    SE      DF
 1     3.115 0.094 223.454
 2     3.060 0.101 199.325
 3     3.005 0.097 200.509
 4     3.235 0.096 196.903

txt> print("R matrix (variances and covariances) and correlations for first person")
[1] "R matrix (variances and covariances) and correlations for first person"
txt> getVarCov(Gaming2, individual = "1", type = "marginal")
Marginal variance covariance matrix
        [,1]    [,2]    [,3]    [,4]
[1,] 1.75050 0.90261 0.74314 0.98791
[2,] 0.90261 2.03660 0.99970 0.82000
[3,] 0.74314 0.99970 1.87430 0.82294
[4,] 0.98791 0.82000 0.82294 1.85910
  Standard Deviations: 1.3231 1.4271 1.3691 1.3635 
txt> corMatrix(Gaming2$modelStruct$corStruct)[[1]]
          [,1]      [,2]      [,3]      [,4]
[1,] 1.0000000 0.4780411 0.4102613 0.5476291
[2,] 0.4780411 1.0000000 0.5116732 0.4214193
[3,] 0.4102613 0.5116732 1.0000000 0.4408529
[4,] 0.5476291 0.4214193 0.4408529 1.0000000
txt> print("Gaming Model 3a: Random Linear Change over Exact Time")
[1] "Gaming Model 3a: Random Linear Change over Exact Time"
txt> Gaming3a = lmer(data = Example5b, REML = TRUE, formula = game ~ 
txt+ 1 + time + (1 + time | PersonID))
txt> lmerSummary(Gaming3a, LRT = TRUE)

Model Fit for Predicting game:
      AIC       BIC    logLik -2*logLik 
 2593.563  2621.671 -1290.782  2581.563 

Model for the Variance:
 Groups    Name      Variance Std.Dev   Corr
 PersonID  Intercept    0.884   0.940       
           time         0.000   0.002 -1.000
 Residual               1.005   1.003       
Number of obs: 800 , groups: PersonID, 200 

Overall F-test for All Fixed Slopes:
 Model    NumDF   DenDF F.value p.value
 Gaming3a     1 599.362   1.010   0.315

Fixed Effects:
 Term         Est    SE      DF t.value p.value
 Intercept  3.056 0.089 214.106  34.321   <.001
 time       0.032 0.032 599.362   1.005   0.315

LRT for random effect(s):
 Random Effect Dropped  Npar    LL -2LL  AIC  BIC   LRT DF p.value p.mix
 <none>                    6 -1291 2582 2594 2622                       
 time | PersonID           4 -1291 2582 2590 2608 0.003  2   0.999 0.978

txt> marginalV(Gaming3a, times = 0:3)

Model-Implied Marginal Variance-Covariance Matrix (V) for Gaming3a:
      0     1     2     3
0 1.889 0.882 0.880 0.878
1 0.882 1.885 0.878 0.877
2 0.880 0.878 1.882 0.875
3 0.878 0.877 0.875 1.879

Model-Implied Marginal Correlation Matrix:
      0     1     2     3
0 1.000 0.467 0.467 0.466
1 0.467 1.000 0.466 0.466
2 0.467 0.466 1.000 0.465
3 0.466 0.466 0.465 1.000

(Rows and columns are values of time)

txt> print("Model 5: Separate Within-Person (WPgame) and Between-Person (PMgame3) Effects of Gaming")
[1] "Model 5: Separate Within-Person (WPgame) and Between-Person (PMgame3) Effects of Gaming"
txt> Model5 = lmer(data = Example5b, REML = TRUE, formula = dep ~ 
txt+ 1 + time * baseage10 + WPgame + time * PMgame3 + (1 + time | 
txt+ 
txt+ PersonID))
txt> lmerSummary(Model5)

Model Fit for Predicting dep:
      AIC       BIC    logLik -2*logLik 
 3671.165  3722.696 -1824.583  3649.165 

Model for the Variance:
 Groups    Name      Variance Std.Dev   Corr
 PersonID  Intercept    6.739   2.596       
           time         0.295   0.543 -0.101
 Residual               2.720   1.649       
Number of obs: 800 , groups: PersonID, 200 

Overall F-test for All Fixed Slopes:
 Model  NumDF   DenDF F.value p.value
 Model5     6 227.520 108.809   <.001

Fixed Effects:
 Term               Est    SE      DF t.value p.value
 Intercept        8.068 0.211 196.559  38.233   <.001
 time             1.006 0.066 198.477  15.313   <.001
 baseage10       -1.373 0.324 196.143  -4.231   <.001
 WPgame           1.283 0.071 561.181  18.111   <.001
 PMgame3          0.164 0.198 196.018   0.827   0.409
 time:baseage10   1.055 0.101 199.395  10.431   <.001
 time:PMgame3     0.016 0.062 199.802   0.252   0.801

txt> PseudoR2(baseModel = Model3a, altModel1 = Model4, altModel2 = Model5)

Pseudo-R2 and Change in Pseudo-R2
 Term       Model   Variance Pseudo.R2 Change.R2
 Intercept  Model3a    7.045                    
 Intercept  Model4     6.270     0.110          
 Intercept  Model5     6.739     0.043    -0.066
 time       Model3a    0.793                    
 time       Model4     0.342     0.568          
 time       Model5     0.295     0.628     0.060
 Residual   Model3a    4.341                    
 Residual   Model4     4.341     0.000          
 Residual   Model5     2.720     0.373     0.373

txt> TotalR2(altModel1 = Model3a, altModel2 = Model4, altModel3 = Model5)

Total-R2 and Change in Total-R2
 Model   Total.R2 Change.R2
 Model3a    0.089          
 Model4     0.136     0.047
 Model5     0.238     0.101

txt> Persons5b = unique(Example5b[, c("PersonID", "baseage10", "PMgame3")])
txt> Agesd = sd(Persons5b$baseage10)
txt> Agemean = mean(Persons5b$baseage10)
txt> Gamesd = sd(Persons5b$PMgame3)
txt> Gamemean = mean(Persons5b$PMgame3)
txt> DepLimits = c(5, 15)
txt> DepBreaks = 5:15
txt> print("Check predicted outcomes for baseline age by person mean gaming over time")
[1] "Check predicted outcomes for baseline age by person mean gaming over time"
txt> GridPlot = expand.grid(time = c(0, 1, 2, 3), baseage10 = Agemean + 
txt+ c(-Agesd, Agesd), PMgame3 = Gamemean + c(-Gamesd, Gamesd), 
txt+ WPgame = 0)
txt> GridPlot$dep = predict(Model5, newdata = GridPlot, re.form = NA)
txt> print(GridPlot)
   time baseage10    PMgame3 WPgame       dep
1     0  -0.74896 -0.9590629      0  8.939379
2     1  -0.74896 -0.9590629      0  9.140092
3     2  -0.74896 -0.9590629      0  9.340804
4     3  -0.74896 -0.9590629      0  9.541516
5     0   0.54626 -0.9590629      0  7.161576
6     1   0.54626 -0.9590629      0  8.728812
7     2   0.54626 -0.9590629      0 10.296049
8     3   0.54626 -0.9590629      0 11.863285
9     0  -0.74896  1.1665629      0  9.286939
10    1  -0.74896  1.1665629      0  9.520719
11    2  -0.74896  1.1665629      0  9.754499
12    3  -0.74896  1.1665629      0  9.988279
13    0   0.54626  1.1665629      0  7.509135
14    1   0.54626  1.1665629      0  9.109440
15    2   0.54626  1.1665629      0 10.709744
16    3   0.54626  1.1665629      0 12.310048
txt> GridPlot$Group = factor(paste(ifelse(GridPlot$baseage10 < Agemean, 
txt+ "Younger", "Older"), ifelse(GridPlot$PMgame3 < Gamemean, 
txt+ "Low", "High")), levels = c("Younger Low", "Younger High", 
txt+ "Older Low", "Older High"))
txt> legendLabels = c("Younger (-1 SD), Low Gaming (-1 SD)", "Younger (-1 SD), High Gaming (+1 SD)", 
txt+ "Older (+1 SD), Low Gaming (-1 SD)", "Older (+1 SD), High Gaming (+1 SD)")
txt> Figure1 = ggplot(GridPlot, aes(x = time, y = dep, color = Group, 
txt+ linetype = Group, shape = Group, fill = Group)) + geom_line(linewidth = 1.2) + 
txt+ geom_point(size = 3.5, stroke = 1) + scale_color_manual(values = c("blue", 
txt+ "blue", "red", "red"), labels = legendLabels) + scale_linetype_manual(values = c("dashed", 
txt+ "solid", "dashed", "solid"), labels = legendLabels) + scale_shape_manual(values = c(22, 
txt+ 22, 24, 24), labels = legendLabels) + scale_fill_manual(values = c("white", 
txt+ "blue", "white", "red"), labels = legendLabels) + scale_x_continuous(breaks = c(0, 
txt+ 1, 2, 3), expand = expansion(mult = 0.08)) + scale_y_continuous(limits = DepLimits, 
txt+ breaks = DepBreaks) + labs(x = "Years Since Baseline", y = "Depression", 
txt+ title = expression(bold("Depression Over Time by Baseline Age and Person Mean Gaming")), 
txt+ subtitle = expression(bolditalic("(at the person's own mean daily gaming)"))) + 
txt+ guides(color = guide_legend(nrow = 2, byrow = TRUE)) + theme_classic(base_size = 14) + 
txt+ theme(plot.title = element_text(face = "bold", hjust = 0.5, 
txt+ 
txt+ size = 16), plot.subtitle = element_text(face = "bold.italic", 
txt+ 
txt+ hjust = 0.5, size = 14), axis.title = element_text(face = "bold", 
txt+ 
txt+ size = 12), axis.text = element_text(size = 12, color = "black"), 
txt+ 
txt+ panel.grid.major.y = element_line(color = "black", linewidth = 0.4), 
txt+ 
txt+ panel.border = element_rect(color = "black", fill = NA, 
txt+ 
txt+ 
txt+ linewidth = 0.8), legend.position = "top", legend.title = element_blank(), 
txt+ 
txt+ legend.key.width = unit(1.6, "cm"), legend.text = element_text(size = 11), 
txt+ 
txt+ legend.box.margin = margin(0, 0, -6, 0), legend.margin = margin(0, 
txt+ 
txt+ 
txt+ 0, 0, 0), legend.key.height = unit(0.5, "cm"), plot.margin = margin(4, 
txt+ 
txt+ 
txt+ 10, 4, 5)) + theme(plot.title = element_text(margin = margin(0, 
txt+ 0, 0, 0)), plot.subtitle = element_text(margin = margin(0, 
txt+ 0, 0, 0)))
txt> ggsave("Example5b_Depression_by_BaseAge_Gaming.png", Figure1, 
txt+ width = 9.5, height = 5.5, dpi = 300)
txt> print(Figure1)
txt> WPsd = sd(Example5b$WPgame)
txt> BPsd = sd(Persons5b$PMgame3)
txt> BPmean = mean(Persons5b$PMgame3)
txt> print("Check predicted outcomes for the WP and BP gaming effects")
[1] "Check predicted outcomes for the WP and BP gaming effects"
txt> GridPlot2 = expand.grid(WPgame = c(-WPsd, 0, WPsd), PMgame3 = BPmean + 
txt+ c(-BPsd, 0, BPsd), time = 0, baseage10 = 0)
txt> GridPlot2$dep = predict(Model5, newdata = GridPlot2, re.form = NA)
txt> print(GridPlot2)
     WPgame    PMgame3 time baseage10      dep
1 -0.868911 -0.9590629    0         0 6.796698
2  0.000000 -0.9590629    0         0 7.911366
3  0.868911 -0.9590629    0         0 9.026034
4 -0.868911  0.1037500    0         0 6.970477
5  0.000000  0.1037500    0         0 8.085146
6  0.868911  0.1037500    0         0 9.199814
7 -0.868911  1.1665629    0         0 7.144257
8  0.000000  1.1665629    0         0 8.258925
9  0.868911  1.1665629    0         0 9.373594
txt> GridPlot2$Group = factor(rep(c("Low", "Average", "High"), each = 3), 
txt+ levels = c("Low", "Average", "High"))
txt> legendLabels2 = c("Low Mean Gaming (-1 SD)", "Average Mean Gaming", 
txt+ "High Mean Gaming (+1 SD)")
txt> Figure2 = ggplot(GridPlot2, aes(x = WPgame, y = dep, color = Group, 
txt+ shape = Group, fill = Group)) + geom_line(linewidth = 1.2) + 
txt+ geom_point(size = 3.5, stroke = 1) + scale_color_manual(values = c("blue", 
txt+ "green3", "red"), labels = legendLabels2) + scale_shape_manual(values = c(21, 
txt+ 24, 22), labels = legendLabels2) + scale_fill_manual(values = c("blue", 
txt+ "green3", "red"), labels = legendLabels2) + scale_x_continuous(breaks = c(-WPsd, 
txt+ 0, WPsd), labels = c("Lower than Usual", "Usual", "Higher than Usual"), 
txt+ expand = expansion(mult = 0.12)) + scale_y_continuous(limits = DepLimits, 
txt+ breaks = DepBreaks) + labs(x = "Within-Person Gaming", y = "Depression", 
txt+ title = expression(bold("Depression by Within-Person and Between-Person Gaming")), 
txt+ subtitle = expression(bolditalic("(at baseline for a 10-year-old, from the Model 5 fixed effects)"))) + 
txt+ theme_classic(base_size = 14) + theme(plot.title = element_text(face = "bold", 
txt+ hjust = 0.5, size = 16), plot.subtitle = element_text(face = "bold.italic", 
txt+ hjust = 0.5, size = 14), axis.title = element_text(face = "bold", 
txt+ size = 12), axis.text = element_text(size = 12, color = "black"), 
txt+ panel.grid.major.y = element_line(color = "black", linewidth = 0.4), 
txt+ panel.border = element_rect(color = "black", fill = NA, linewidth = 0.8), 
txt+ legend.position = "top", legend.title = element_blank(), 
txt+ legend.key.width = unit(1.6, "cm"), legend.text = element_text(size = 11), 
txt+ legend.box.margin = margin(0, 0, -6, 0), legend.margin = margin(0, 
txt+ 
txt+ 0, 0, 0), legend.key.height = unit(0.5, "cm"), plot.margin = margin(4, 
txt+ 
txt+ 10, 4, 5)) + theme(plot.title = element_text(margin = margin(0, 
txt+ 0, 0, 0)), plot.subtitle = element_text(margin = margin(0, 
txt+ 0, 0, 0)))
txt> ggsave("Example5b_Depression_by_WP_BP_Gaming.png", Figure2, width = 9.5, 
txt+ height = 5.5, dpi = 300)
txt> print(Figure2)
