### Week 07
### Longitudinal Mixed Effects
### First, set the working directory where R will find the data
setwd("H:/data/User/Classes/multiv_13/YTs/YTSplit/Week 07 wmv - 2017/Wk07.R")
getwd()## [1] "H:/data/User/Classes/multiv_13/YTs/YTSplit/Week 07 wmv - 2017/Wk07.R"### Read the data
class07<- read.csv("Wk07_Inclassdemo.csv")
head(class07)library(reshape)
library(plotrix)
library(plyr)##
## Attaching package: 'plyr'## The following objects are masked from 'package:reshape':
##
## rename, round_anyrequire(dplyr)## Loading required package: dplyr##
## Attaching package: 'dplyr'## The following objects are masked from 'package:plyr':
##
## arrange, count, desc, failwith, id, mutate, rename, summarise,
## summarize## The following object is masked from 'package:reshape':
##
## rename## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, unionrequire(ggplot2)## Loading required package: ggplot2require(graphics)
require(lme4)## Loading required package: lme4## Loading required package: Matrix##
## Attaching package: 'Matrix'## The following object is masked from 'package:reshape':
##
## expandrequire(nlme)## Loading required package: nlme##
## Attaching package: 'nlme'## The following object is masked from 'package:lme4':
##
## lmList## The following object is masked from 'package:dplyr':
##
## collapserequire(saemix)## Loading required package: saemix## Loading package saemix, version 2.3, December 2019
## please direct bugs, questions and feedback to emmanuelle.comets@inserm.frrequire(gee)## Loading required package: geerequire(matrixStats)## Loading required package: matrixStats##
## Attaching package: 'matrixStats'## The following object is masked from 'package:dplyr':
##
## count## The following object is masked from 'package:plyr':
##
## countrequire(mice)## Loading required package: mice##
## Attaching package: 'mice'## The following objects are masked from 'package:base':
##
## cbind, rbindrequire(nlmeU)## Loading required package: nlmeU##
## Attaching package: 'nlmeU'## The following object is masked from 'package:stats':
##
## sigmalibrary(Hmisc)## Loading required package: lattice## Loading required package: survival## Loading required package: Formula##
## Attaching package: 'Hmisc'## The following objects are masked from 'package:dplyr':
##
## src, summarize## The following objects are masked from 'package:plyr':
##
## is.discrete, summarize## The following objects are masked from 'package:base':
##
## format.pval, unitslibrary(MASS)##
## Attaching package: 'MASS'## The following object is masked from 'package:dplyr':
##
## selectlibrary(clubSandwich)## Registered S3 method overwritten by 'clubSandwich':
## method from
## bread.mlm sandwichlibrary(lmerTest)##
## Attaching package: 'lmerTest'## The following object is masked from 'package:lme4':
##
## lmer## The following object is masked from 'package:stats':
##
## steplibrary(RLRsim)
library(car)## Loading required package: carData## Registered S3 methods overwritten by 'car':
## method from
## influence.merMod lme4
## cooks.distance.influence.merMod lme4
## dfbeta.influence.merMod lme4
## dfbetas.influence.merMod lme4##
## Attaching package: 'car'## The following object is masked from 'package:dplyr':
##
## recode###-------------------------------------------------------------------------------------
### 0. Spaghetti plots
###-------------------------------------------------------------------------------------
#=================================================
# We will use the base plotting function plot( ) and
# the function ggplot( ) from the ggplot2 packdep.
# Both are valuable tools for exploring and creating
# publishable quality graphics. We will not be using
# the xyplot( ) function--although this can be used
# for the same purposes.
# Sometimes it may be difficult to see individuals'
# growth trajectories clearly when there are many
# observations in the data. In these scenarios it
# is wise to take a few random samples of a smaller
# size.
# Several individuals data connected by lines across
# time is called a "spaghetti plot" because it often
# resembles that famed pasta. To get the program to
# create this type of graph correctly, all NAs must
# be removed from the dataset before graphing. This
# is critical.
#==================================================
n.dat <- na.omit(class07)
head(n.dat, n=10)### Take a random sample of size = n
ids <- unique(class07$ID) ### This creates a vector of identification numbers in the data
ids## [1] 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18
## [19] 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36
## [37] 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54
## [55] 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72
## [73] 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90
## [91] 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108
## [109] 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126
## [127] 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144
## [145] 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162
## [163] 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180
## [181] 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198
## [199] 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216
## [217] 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234
## [235] 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252
## [253] 253 254 255 256 257 258 259 260 261 262 263 264 265 266 267 268 269 270
## [271] 271 272 273 274 275 276 277 278 279 280 281 282 283 284 285 286 287 288
## [289] 289 290 291 292 293 294 295 296 297 298 299 300 301 302 303 304 305 306
## [307] 307 308 309 310 311 312 313 314 315 316 317 318 319 320 321 322 323 324
## [325] 325 326 327 328 329 330 331 332 333 334 335 336 337 338 339 340 341 342
## [343] 343 344 345 346 347 348 349 350 351 352 353 354 355 356 357 358 359 360
## [361] 361 362 363 364 365 366 367 368 369 370 371 372 373 374 375 376 377 378
## [379] 379 380 381 382 383 384 385 386 387 388 389 390 391 392 393 394 395 396
## [397] 397 398 399 400subdata <- subset(class07,
ID %in% sample(ids,size=30))
### The second line randomly picks 30 identification numbers from "ids"
#subdata
### The subject identification number, ID, is used to associate subjects with their
### repeated measures. In the syntax below, the graph layers are saved in the objects:
### g1(first layer) and g2(second layer). Like layers of paint on a canvas.
### The object names for the layers are arbitrary.
### The g2 object is an accumulation of the layers and is printed in the last line of syntax.
g1 <- ggplot(data=na.omit(subdata), ### omit the missing values
aes(x=Occasion,y=posaff,group=ID))
print(g1)
### In the syntax for g1, the "data=" identifies the data frame. The aesthetic mapping, aes,
### identifies the variables for the axes. "x=Occasion" indicates the x-axis, "y=Security" indicates
### the y-axis, and "group=ID" indicates that the grouping variable is the subject identification
### number meaning a line will be drawn for each subject.
### Here, geom_line() indicates that a line will be drawn for each subject, geom_point()
### will be individuals' repeated measures and ggtitle() is used to
### manipulate the plot title.
### The DEFAULT graph has a gray background, and the axis labels indicate the variable names.
### In order to change the default options, we will add three components in two additional layers. That is,
### theme_bw() that changes the color theme to black-and-white, scale_x_continuous() that controls
### the scaling and labeling of the x-axis, and scale_y_continuous() that does the same for the
###y-axis. The graph is built in layers with g1 and g2 containing the same elements as previously, and
### the added options in g3, and g4 and g5.
g2 <- g1 + geom_line () + geom_point () + ggtitle("posaff Spaghetti Plot")
print(g2)
g3 <- g2 + theme_bw ()
print(g3)
lab.occ = c("0","1","2","3","4")
g4 <- g3 + scale_x_continuous (breaks = c(0,1,2,3,4), labels=lab.occ, name= "Occasion")
print(g4)
g5 <- g4 + scale_y_continuous (breaks = seq(-180,420,50), name= "posaff")
print(g5)
### Superimposing summary stats
#=================================================
# To get some idea of the functional form of the mean
# trajectory, you may want to first fit the empirical
# means of the data and connect them with a line. The
# curve of the OBSERVED MEANS at each time point can
# be drawn with stat_summary() using the options
# "fun.y=mean" and "geom="line". The "fun.y=mean" option
# computes the mean of the response variable
# at each time point and geom="line" draws the line.
# The means can be drawn as points by replacing
# geom="line" with geom="point".
# The syntax below draws the individual's observed
# data points and the observed mean curves.
#=================================================
g7 <- g5 + stat_summary(fun.y=mean,geom="point",size=2.5,shape=18,colour="orange")## Warning: `fun.y` is deprecated. Use `fun` instead.g8 <- g7 + stat_summary(fun.y=mean,geom="line",lwd=0.4,colour="blue")## Warning: `fun.y` is deprecated. Use `fun` instead.print(g8)
### The stat_smooth() function allows you to insert a particular functional form
### that you identify to see how well it fits the empirical means
g9 <- g8 + stat_smooth(method="lm",se=FALSE,lwd=0.2,formula=y~x,colour="red")
print(g9)
g10 <- g8 + stat_smooth(method="lm",se=FALSE,lwd=0.2,formula=y~poly(x,2),colour="green")
print(g10)
### In the syntax above, the option "lwd=" controls the
### line width and the options "size=" and "shape="
### control the size and type of the points, respectively.
###
### See the following webpdep for the possible line types and shape options.
### http://wiki.stdout.org/rcookbook/Graphs/Shapes%20and%20line%20types/
### Plot individual data
#=================================================
# Facet plots of individual curves using observed data
# points. To better examine the shape of individual
# trajectories, each subject's change curve can be
# depicted (not fitted) in its own graph. When several
# of these graphs appear in the same figure the graphs
# are known as panes, panels, or facets. Faceting for
# individual subjects is accomplished using facet_wrap().
# The syntax convention is to list the variable on which
# the faceting is based, preceded by a tilde(~). The
# following syntax produces a facet for each subject.
# We draw points for the observed values using geom_point()
# in addition to drawing lines. With numerous subjects a
# facet plot might have too many panels to be practical
# and it will be very busy. For easier comprehension,
# it is desirable to select a RANDOM subset of the sample for graphing.
# Like before, a random sample or subset of the larger
# number of unique cases can be chosen. This creates a
# vector of identification numbers in the data.
#===================================================
ids <- unique(class07$ID)
### Randomly select 6 observations from the data
subdata2 <- subset(class07,
ID %in% sample(ids,size=6)) ### The second line randomly picks 6 identification numbers from "ids"
#subdata2
### Plot the growth curves for each of the 6 individual
g11 <- ggplot(data=na.omit(subdata2),aes(x=Occasion,y=posaff,group=ID))
print(g11)
g12 <- g11 + geom_line () + geom_point ()
print(g12)
g13 <- g12 + facet_wrap(~ID,nrow=3,ncol=2) ### The number of rows or columns can be fixed using ncol= and nrow=
print(g13)
g14 <- g13 + theme_bw ()
print(g14)
g15 <- g14 + scale_x_continuous (breaks = c(0,1,2,3,4), labels=lab.occ, name= "Occasion")
print(g15)
g16 <- g15 + scale_y_continuous (breaks = seq(-180,420,50), name= "posaff")
print (g16)
#Conditioning on time-invariant covariates
### Spaghetti plot with lowess smoother by deprgroup to
## get a sense of the mean function by deprgroup status
sec.v=na.omit(class07)
par(mfrow=c(1,2),oma=c(0,0,0,0),adj=0.5,font.sub=3)
tmp1 = na.omit(cbind(sec.v$Occasion[sec.v$deprgroup=="0"],sec.v$posaff[sec.v$deprgroup=="0"]))
tmp2 = na.omit(cbind(sec.v$Occasion[sec.v$deprgroup=="1"],sec.v$posaff[sec.v$deprgroup=="1"]))
scatter.smooth(tmp1[,1],tmp1[,2],pch=42,cex=0.75,xlab="Occasion",
ylab="posaff")
title("deprgroup = 0")
scatter.smooth(tmp2[,1],tmp2[,2],pch=42,cex=0.75,xlab="Occasion",
ylab="posaff")
title("deprgroup = 1") 
#From last class, we can adopt techniques to look at scatter plots
#by deprgroup status
class07$dep2<-as.factor(class07$deprgroup)
#Here, we show occasion effect by dep group
qplot(y=posaff, x=Occasion, colour = dep2,
data = class07)
#Here, we superimpose best fitting line by und2
ggplot(class07) +
aes(x=class07$Occasion, y=class07$posaff, col=as.factor(class07$dep2)) + #
geom_point() +
stat_smooth(method="loess", se = FALSE)## Warning: Use of `class07$Occasion` is discouraged. Use `Occasion` instead.## Warning: Use of `class07$posaff` is discouraged. Use `posaff` instead.## Warning: Use of `class07$dep2` is discouraged. Use `dep2` instead.## Warning: Use of `class07$Occasion` is discouraged. Use `Occasion` instead.## Warning: Use of `class07$posaff` is discouraged. Use `posaff` instead.## Warning: Use of `class07$dep2` is discouraged. Use `dep2` instead.## `geom_smooth()` using formula 'y ~ x'
#Here, we vanish the data points in previous plot so we can see lines better; geom_point removed
ggplot(class07) +
aes(x=class07$Occasion, y=class07$posaff, col=as.factor(class07$dep2)) + #
stat_smooth(method="loess", se = FALSE)## Warning: Use of `class07$Occasion` is discouraged. Use `Occasion` instead.## Warning: Use of `class07$posaff` is discouraged. Use `posaff` instead.## Warning: Use of `class07$dep2` is discouraged. Use `dep2` instead.## `geom_smooth()` using formula 'y ~ x'
#removing grey background
ggplot(class07) +
aes(x=class07$Occasion, y=class07$posaff, col=as.factor(class07$dep2)) + #
stat_smooth(method="loess", se = FALSE) +
theme_bw ()## Warning: Use of `class07$Occasion` is discouraged. Use `Occasion` instead.## Warning: Use of `class07$posaff` is discouraged. Use `posaff` instead.## Warning: Use of `class07$dep2` is discouraged. Use `dep2` instead.## `geom_smooth()` using formula 'y ~ x'
#Faceting by dep group
g31 <- ggplot(data=na.omit(class07),aes(x=Occasion,y=posaff,group=ID)) +
geom_line ()+ theme_bw ()
g31
g32 <- g31 + stat_summary(data=na.omit(class07),aes(x=Occasion,y=posaff,group=dep2),fun.y=mean,geom="line",
lwd=1.2,colour="blue")## Warning: `fun.y` is deprecated. Use `fun` instead.g32
g33 <- g32 + facet_grid (.~dep2) ### Column faceting
print(g33)
g34 <- g32 + facet_grid (dep2~.) ### Row faceting
print(g34)
#=================================================
# #2
# FITTING RCM (ie., LGM) WITH
# lme() IN nlme PACKdep
# lmer() IN lme4 PACKdep
#=================================================
### fit an RCM using REML using the lme() function
###------0. Intercept only model (unconditional means/null) ----------- ###
d.1.out <- lme(fixed=posaff~1,
data=class07,
na.action=na.omit,
method="ML",
random=~1|ID)
summary(d.1.out)## Linear mixed-effects model fit by maximum likelihood
## Data: class07
## AIC BIC logLik
## 20765.79 20782.6 -10379.9
##
## Random effects:
## Formula: ~1 | ID
## (Intercept) Residual
## StdDev: 75.80368 30.76856
##
## Fixed effects: posaff ~ 1
## Value Std.Error DF t-value p-value
## (Intercept) 88.74896 3.853086 1600 23.03322 0
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -2.69132098 -0.52341602 -0.01714419 0.51002832 3.00352164
##
## Number of Observations: 2000
## Number of Groups: 400names(d.1.out)## [1] "modelStruct" "dims" "contrasts" "coefficients" "varFix"
## [6] "sigma" "apVar" "logLik" "numIter" "groups"
## [11] "call" "terms" "method" "fitted" "residuals"
## [16] "fixDF" "na.action" "data"###------1. Fixed and random linear time (unconditional growth) ----------- ###
d.2.out <- lme(fixed=posaff~1+Occasion,
data=class07,
na.action=na.omit,
method="ML",
#correlation=corCompSymm(), #options are corrAR1, corCAR1
random=~1+Occasion|ID)
summary(d.2.out)## Linear mixed-effects model fit by maximum likelihood
## Data: class07
## AIC BIC logLik
## 16913.06 16946.67 -8450.532
##
## Random effects:
## Formula: ~1 + Occasion | ID
## Structure: General positive-definite, Log-Cholesky parametrization
## StdDev Corr
## (Intercept) 70.547192 (Intr)
## Occasion 10.859323 0.156
## Residual 6.104389
##
## Fixed effects: posaff ~ 1 + Occasion
## Value Std.Error DF t-value p-value
## (Intercept) 57.38966 3.537043 1599 16.22532 0
## Occasion 15.67965 0.551754 1599 28.41782 0
## Correlation:
## (Intr)
## Occasion 0.144
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -2.659542809 -0.523680035 0.006297472 0.504120663 2.627962416
##
## Number of Observations: 2000
## Number of Groups: 400Anova(d.2.out)#?corClasses
###------2. Add fixed deprgroup (conditional intercept) ----------- ###
d.3.out <- lme(fixed=posaff~1+Occasion+deprgroup,
data=class07,
na.action=na.omit,
method="ML",
#correlation=corCompSymm(), #options are corrAR1, corCAR1
random=~1+Occasion|ID)
summary(d.3.out)## Linear mixed-effects model fit by maximum likelihood
## Data: class07
## AIC BIC logLik
## 16897.67 16936.87 -8441.834
##
## Random effects:
## Formula: ~1 + Occasion | ID
## Structure: General positive-definite, Log-Cholesky parametrization
## StdDev Corr
## (Intercept) 68.381874 (Intr)
## Occasion 10.859323 0.06
## Residual 6.104389
##
## Fixed effects: posaff ~ 1 + Occasion + deprgroup
## Value Std.Error DF t-value p-value
## (Intercept) 73.75605 4.909249 1599 15.023895 0
## Occasion 15.67965 0.551892 1599 28.410706 0
## deprgroup -31.93440 6.853473 398 -4.659595 0
## Correlation:
## (Intr) Occasn
## Occasion 0.034
## deprgroup -0.715 0.000
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -2.663896962 -0.526782654 0.006777205 0.502702976 2.626452784
##
## Number of Observations: 2000
## Number of Groups: 400Anova(d.3.out)###------3. Add fixed Level 2 mindful (mean) (conditional intercept) ----------- ###
d.4.out <- lme(fixed=posaff~1+Occasion+deprgroup+mindful_mean,
data=class07,
na.action=na.omit,
method="ML",
#correlation=corCompSymm(), #options are corrAR1, corCAR1
random=~1+Occasion|ID)
summary(d.4.out)## Linear mixed-effects model fit by maximum likelihood
## Data: class07
## AIC BIC logLik
## 16882.81 16927.62 -8433.404
##
## Random effects:
## Formula: ~1 + Occasion | ID
## Structure: General positive-definite, Log-Cholesky parametrization
## StdDev Corr
## (Intercept) 66.950706 (Intr)
## Occasion 10.859323 0.061
## Residual 6.104389
##
## Fixed effects: posaff ~ 1 + Occasion + deprgroup + mindful_mean
## Value Std.Error DF t-value p-value
## (Intercept) -6.836137 20.029839 1599 -0.341298 0.7329
## Occasion 15.679648 0.552030 1599 28.403591 0.0000
## deprgroup -28.773429 6.756210 397 -4.258812 0.0000
## mindful_mean 2.725558 0.657462 397 4.145578 0.0000
## Correlation:
## (Intr) Occasn dprgrp
## Occasion 0.008
## deprgroup -0.281 0.000
## mindful_mean -0.971 0.000 0.114
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -2.652342428 -0.524702532 0.006122746 0.499551177 2.627974326
##
## Number of Observations: 2000
## Number of Groups: 400Anova(d.4.out)###------4. Add fixed & random Level 1 mindful (centered) (time-varying) ----------- ###
d.5.out <- lme(fixed=posaff~1+Occasion+deprgroup+mindful_mean+mindful_center,
data=class07,
na.action=na.omit,
method="ML",
#correlation=corCompSymm(), #options are corrAR1, corCAR1
random=~1+Occasion+mindful_center|ID)
summary(d.5.out)## Linear mixed-effects model fit by maximum likelihood
## Data: class07
## AIC BIC logLik
## 16851.36 16918.57 -8413.68
##
## Random effects:
## Formula: ~1 + Occasion + mindful_center | ID
## Structure: General positive-definite, Log-Cholesky parametrization
## StdDev Corr
## (Intercept) 66.6809160 (Intr) Occasn
## Occasion 11.1362304 0.064
## mindful_center 0.7845291 -0.102 -0.301
## Residual 5.9355925
##
## Fixed effects: posaff ~ 1 + Occasion + deprgroup + mindful_mean + mindful_center
## Value Std.Error DF t-value p-value
## (Intercept) -5.047103 20.023269 1598 -0.252062 0.801
## Occasion 14.431765 0.602929 1598 23.936107 0.000
## deprgroup -28.715755 6.754063 397 -4.251627 0.000
## mindful_mean 2.748929 0.657277 397 4.182301 0.000
## mindful_center 0.834785 0.144579 1598 5.773917 0.000
## Correlation:
## (Intr) Occasn dprgrp mndfl_m
## Occasion 0.003
## deprgroup -0.281 -0.002
## mindful_mean -0.971 -0.002 0.114
## mindful_center 0.011 -0.409 0.004 0.004
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -2.752041519 -0.499103111 -0.001944696 0.493204914 2.659190914
##
## Number of Observations: 2000
## Number of Groups: 400Anova(d.5.out)###------5. Add fixed dep * occasion interaction (conditional slope) ----------- ###
d.6.out <- lme(fixed=posaff~1+Occasion+deprgroup+mindful_mean+mindful_center+OccDepr,
data=class07,
na.action=na.omit,
method="ML",
#correlation=corCompSymm(), #options are corrAR1, corCAR1
random=~1+Occasion+mindful_center|ID)
summary(d.6.out)## Linear mixed-effects model fit by maximum likelihood
## Data: class07
## AIC BIC logLik
## 16772.39 16845.2 -8373.194
##
## Random effects:
## Formula: ~1 + Occasion + mindful_center | ID
## Structure: General positive-definite, Log-Cholesky parametrization
## StdDev Corr
## (Intercept) 66.685327 (Intr) Occasn
## Occasion 10.045890 0.061
## mindful_center 0.789696 -0.091 -0.295
## Residual 5.934118
##
## Fixed effects: posaff ~ 1 + Occasion + deprgroup + mindful_mean + mindful_center + OccDepr
## Value Std.Error DF t-value p-value
## (Intercept) 6.03727 20.063189 1597 0.300913 0.7635
## Occasion 14.43178 0.552892 1597 26.102367 0.0000
## deprgroup -50.33435 7.132805 397 -7.056740 0.0000
## mindful_mean 2.74876 0.657460 397 4.180880 0.0000
## mindful_center 0.83505 0.144525 1597 5.777913 0.0000
## OccDepr -9.41587 0.995338 1597 -9.459965 0.0000
## Correlation:
## (Intr) Occasn dprgrp mndfl_m mndfl_c
## Occasion 0.002
## deprgroup -0.284 -0.002
## mindful_mean -0.969 -0.002 0.108
## mindful_center 0.011 -0.436 0.004 0.004
## OccDepr -0.058 -0.001 0.321 0.000 0.002
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -2.732249099 -0.497489705 0.001225253 0.499338011 2.664816653
##
## Number of Observations: 2000
## Number of Groups: 400Anova(d.6.out)anova(d.1.out,d.2.out,d.3.out,d.4.out,d.5.out,d.6.out)### lmer
sec.0 = lmer(posaff ~ 1 +
(1|ID), data=class07, REML=FALSE)
summary(sec.0)## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
## method [lmerModLmerTest]
## Formula: posaff ~ 1 + (1 | ID)
## Data: class07
##
## AIC BIC logLik deviance df.resid
## 20765.8 20782.6 -10379.9 20759.8 1997
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.69132 -0.52342 -0.01714 0.51003 3.00352
##
## Random effects:
## Groups Name Variance Std.Dev.
## ID (Intercept) 5746.2 75.80
## Residual 946.7 30.77
## Number of obs: 2000, groups: ID, 400
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 88.749 3.852 400.000 23.04 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1anova(sec.0)sec.1 = lmer(posaff ~ Occasion +
(1+Occasion|ID), data=class07, REML=FALSE)
summary(sec.1)## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
## method [lmerModLmerTest]
## Formula: posaff ~ Occasion + (1 + Occasion | ID)
## Data: class07
##
## AIC BIC logLik deviance df.resid
## 16913.1 16946.7 -8450.5 16901.1 1994
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.6595 -0.5237 0.0063 0.5041 2.6280
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## ID (Intercept) 4976.91 70.547
## Occasion 117.93 10.859 0.16
## Residual 37.26 6.104
## Number of obs: 2000, groups: ID, 400
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 57.3897 3.5353 399.9995 16.23 <2e-16 ***
## Occasion 15.6796 0.5515 399.9993 28.43 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## Occasion 0.144anova(sec.0,sec.1)sec.2 = lmer(posaff ~ Occasion + deprgroup +
(1+Occasion|ID), data=class07, REML=FALSE)
summary(sec.2)## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
## method [lmerModLmerTest]
## Formula: posaff ~ Occasion + deprgroup + (1 + Occasion | ID)
## Data: class07
##
## AIC BIC logLik deviance df.resid
## 16897.7 16936.9 -8441.8 16883.7 1993
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.66390 -0.52678 0.00678 0.50270 2.62645
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## ID (Intercept) 4676.08 68.382
## Occasion 117.92 10.859 0.06
## Residual 37.26 6.104
## Number of obs: 2000, groups: ID, 400
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 73.7560 4.9056 400.3820 15.035 < 2e-16 ***
## Occasion 15.6796 0.5515 400.0007 28.432 < 2e-16 ***
## deprgroup -31.9343 6.8483 399.9992 -4.663 4.25e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) Occasn
## Occasion 0.034
## deprgroup -0.715 0.000anova(sec.0,sec.1,sec.2)relgrad <- with(sec.2@optinfo$derivs,solve(Hessian,gradient))
max(abs(relgrad)) #note, less than .001, so probably okay## [1] 2.728362e-06sec.3 = lmer(posaff ~ Occasion + deprgroup +
mindful_mean +
(1+Occasion|ID), data=class07, REML=FALSE)
summary(sec.3)## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
## method [lmerModLmerTest]
## Formula: posaff ~ Occasion + deprgroup + mindful_mean + (1 + Occasion |
## ID)
## Data: class07
##
## AIC BIC logLik deviance df.resid
## 16882.8 16927.6 -8433.4 16866.8 1992
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.65235 -0.52470 0.00612 0.49955 2.62798
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## ID (Intercept) 4482.37 66.951
## Occasion 117.93 10.859 0.06
## Residual 37.26 6.104
## Number of obs: 2000, groups: ID, 400
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) -6.8365 20.0097 400.0545 -0.342 0.733
## Occasion 15.6796 0.5515 399.9877 28.432 < 2e-16 ***
## deprgroup -28.7728 6.7494 400.0038 -4.263 2.52e-05 ***
## mindful_mean 2.7256 0.6568 400.0008 4.150 4.07e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) Occasn dprgrp
## Occasion 0.008
## deprgroup -0.281 0.000
## mindful_men -0.971 0.000 0.114anova(sec.0,sec.1,sec.2,sec.3)relgrad <- with(sec.3@optinfo$derivs,solve(Hessian,gradient))
max(abs(relgrad)) #note, less than .001, so probably okay## [1] 2.351304e-05sec.4 = lmer(posaff ~ Occasion + deprgroup +
mindful_mean + mindful_center +
(1+Occasion+mindful_center||ID), data=class07, REML=FALSE)
summary(sec.4)## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
## method [lmerModLmerTest]
## Formula: posaff ~ Occasion + deprgroup + mindful_mean + mindful_center +
## (1 + Occasion + mindful_center || ID)
## Data: class07
##
## AIC BIC logLik deviance df.resid
## 16849.1 16899.5 -8415.6 16831.1 1991
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.73097 -0.50469 -0.00553 0.49777 2.62759
##
## Random effects:
## Groups Name Variance Std.Dev.
## ID (Intercept) 4482.5876 66.9521
## ID.1 Occasion 117.3917 10.8347
## ID.2 mindful_center 0.2718 0.5214
## Residual 35.7212 5.9767
## Number of obs: 2000, groups: ID, 400
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) -3.1005 20.0394 400.1849 -0.155 0.877
## Occasion 14.4339 0.5875 510.8217 24.567 < 2e-16 ***
## deprgroup -32.2185 6.7588 399.9742 -4.767 2.62e-06 ***
## mindful_mean 2.7436 0.6577 399.9939 4.171 3.72e-05 ***
## mindful_center 0.8308 0.1403 235.3303 5.923 1.11e-08 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) Occasn dprgrp mndfl_m
## Occasion -0.007
## deprgroup -0.281 -0.001
## mindful_men -0.971 -0.001 0.114
## mindfl_cntr 0.016 -0.345 0.001 0.004anova(sec.0,sec.1,sec.2,sec.3,sec.4)relgrad <- with(sec.4@optinfo$derivs,solve(Hessian,gradient))
max(abs(relgrad)) #note, less than .001, so probably okay## [1] 7.201962e-06sec.5 = lmer(posaff ~ Occasion + deprgroup +
mindful_mean + mindful_center + OccDepr +
(1+Occasion+mindful_center||ID), data=class07, REML=FALSE)
summary(sec.5)## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
## method [lmerModLmerTest]
## Formula: posaff ~ Occasion + deprgroup + mindful_mean + mindful_center +
## OccDepr + (1 + Occasion + mindful_center || ID)
## Data: class07
##
## AIC BIC logLik deviance df.resid
## 16770.2 16826.2 -8375.1 16750.2 1990
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.70692 -0.50583 -0.00037 0.50346 2.63628
##
## Random effects:
## Groups Name Variance Std.Dev.
## ID (Intercept) 4484.3602 66.9654
## ID.1 Occasion 95.1265 9.7533
## ID.2 mindful_center 0.2478 0.4978
## Residual 35.7576 5.9798
## Number of obs: 2000, groups: ID, 400
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 6.2773 20.0667 402.1190 0.313 0.755
## Occasion 14.4235 0.5379 535.4529 26.815 < 2e-16 ***
## deprgroup -50.4500 7.0245 462.8084 -7.182 2.76e-12 ***
## mindful_mean 2.7431 0.6578 400.0051 4.170 3.74e-05 ***
## mindful_center 0.8376 0.1397 245.5754 5.994 7.25e-09 ***
## OccDepr -9.4438 0.9968 399.0828 -9.474 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) Occasn dprgrp mndfl_m mndfl_c
## Occasion -0.008
## deprgroup -0.284 0.000
## mindful_men -0.969 -0.001 0.110
## mindfl_cntr 0.017 -0.377 0.000 0.004
## OccDepr -0.049 0.001 0.272 0.000 -0.003anova(sec.0,sec.1,sec.2,sec.3,sec.4,sec.5)relgrad <- with(sec.5@optinfo$derivs,solve(Hessian,gradient))
max(abs(relgrad)) #note, less than .001, so probably okay## [1] 3.600084e-06#######################################################################
#With R, using lme4, we do not have the option of specifying a covariance structure.
#If you do a Google search you will find that there are many people asking how to
#specify those structures in lmer, but no answers. The problem here is that mixed models,
#as done by R, does not let you specify a matrix form for these covariances, and that is
#that!
#######################################################################
### LME offers a fair number of options, if limited
### Many SAS/SPSS repeated covariance options are not available
### ?corClasses at bottom gives you all possible options
### And the single vs double line in random effects for lmer let's you
### alternate between variance components and unstructured.
########Comparing REPEATED covariance structures##########
#AR1
d.6.out_ar <- lme(fixed=posaff~1+Occasion+deprgroup+mindful_mean+mindful_center+OccDepr,
data=class07,
na.action=na.omit,
method="ML",
correlation=corAR1(), #options are corCompSymm, corCAR1
random=~1+Occasion+mindful_center|ID)
summary(d.6.out_ar)## Linear mixed-effects model fit by maximum likelihood
## Data: class07
## AIC BIC logLik
## 16774.1 16852.52 -8373.051
##
## Random effects:
## Formula: ~1 + Occasion + mindful_center | ID
## Structure: General positive-definite, Log-Cholesky parametrization
## StdDev Corr
## (Intercept) 66.6752198 (Intr) Occasn
## Occasion 10.0369467 0.062
## mindful_center 0.7858928 -0.093 -0.294
## Residual 6.0068860
##
## Correlation Structure: AR(1)
## Formula: ~1 | ID
## Parameter estimate(s):
## Phi
## 0.02958664
## Fixed effects: posaff ~ 1 + Occasion + deprgroup + mindful_mean + mindful_center + OccDepr
## Value Std.Error DF t-value p-value
## (Intercept) 6.00716 20.063218 1597 0.299412 0.7647
## Occasion 14.43169 0.552896 1597 26.101998 0.0000
## deprgroup -50.32776 7.132915 397 -7.055708 0.0000
## mindful_mean 2.74972 0.657461 397 4.182333 0.0000
## mindful_center 0.83390 0.144525 1597 5.769937 0.0000
## OccDepr -9.41790 0.995553 1597 -9.459972 0.0000
## Correlation:
## (Intr) Occasn dprgrp mndfl_m mndfl_c
## Occasion 0.002
## deprgroup -0.284 -0.002
## mindful_mean -0.969 -0.002 0.108
## mindful_center 0.011 -0.436 0.004 0.004
## OccDepr -0.058 -0.001 0.321 0.000 0.002
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -2.727425377 -0.495471472 0.002211132 0.487799342 2.653679311
##
## Number of Observations: 2000
## Number of Groups: 400Anova(d.6.out_ar)#CompSymp
d.6.out_cs <- lme(fixed=posaff~1+Occasion+deprgroup+mindful_mean+mindful_center+OccDepr,
data=class07,
na.action=na.omit,
method="ML",
correlation=corCompSymm(), #options are corAR1, corCAR1
random=~1+Occasion+mindful_center|ID)
summary(d.6.out_cs)## Linear mixed-effects model fit by maximum likelihood
## Data: class07
## AIC BIC logLik
## 16774.39 16852.8 -8373.194
##
## Random effects:
## Formula: ~1 + Occasion + mindful_center | ID
## Structure: General positive-definite, Log-Cholesky parametrization
## StdDev Corr
## (Intercept) 66.6913620 (Intr) Occasn
## Occasion 10.0458896 0.061
## mindful_center 0.7896939 -0.091 -0.295
## Residual 5.8660148
##
## Correlation Structure: Compound symmetry
## Formula: ~1 | ID
## Parameter estimate(s):
## Rho
## -0.02335467
## Fixed effects: posaff ~ 1 + Occasion + deprgroup + mindful_mean + mindful_center + OccDepr
## Value Std.Error DF t-value p-value
## (Intercept) 6.03727 20.063193 1597 0.300913 0.7635
## Occasion 14.43178 0.552892 1597 26.102369 0.0000
## deprgroup -50.33435 7.132806 397 -7.056739 0.0000
## mindful_mean 2.74876 0.657460 397 4.180880 0.0000
## mindful_center 0.83505 0.144525 1597 5.777914 0.0000
## OccDepr -9.41587 0.995338 1597 -9.459967 0.0000
## Correlation:
## (Intr) Occasn dprgrp mndfl_m mndfl_c
## Occasion 0.002
## deprgroup -0.284 -0.002
## mindful_mean -0.969 -0.002 0.108
## mindful_center 0.011 -0.436 0.004 0.004
## OccDepr -0.058 -0.001 0.321 0.000 0.002
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -2.762842e+00 -5.035448e-01 -5.831915e-07 5.047569e-01 2.696629e+00
##
## Number of Observations: 2000
## Number of Groups: 400Anova(d.6.out_cs)#Scaled identity
d.6.out_si <- lme(fixed=posaff~1+Occasion+deprgroup+mindful_mean+mindful_center+OccDepr,
data=class07,
na.action=na.omit,
method="ML",
#correlation=corCompSymm(), #options are corAR1, corCAR1
random=~1+Occasion+mindful_center|ID)
summary(d.6.out_si)## Linear mixed-effects model fit by maximum likelihood
## Data: class07
## AIC BIC logLik
## 16772.39 16845.2 -8373.194
##
## Random effects:
## Formula: ~1 + Occasion + mindful_center | ID
## Structure: General positive-definite, Log-Cholesky parametrization
## StdDev Corr
## (Intercept) 66.685327 (Intr) Occasn
## Occasion 10.045890 0.061
## mindful_center 0.789696 -0.091 -0.295
## Residual 5.934118
##
## Fixed effects: posaff ~ 1 + Occasion + deprgroup + mindful_mean + mindful_center + OccDepr
## Value Std.Error DF t-value p-value
## (Intercept) 6.03727 20.063189 1597 0.300913 0.7635
## Occasion 14.43178 0.552892 1597 26.102367 0.0000
## deprgroup -50.33435 7.132805 397 -7.056740 0.0000
## mindful_mean 2.74876 0.657460 397 4.180880 0.0000
## mindful_center 0.83505 0.144525 1597 5.777913 0.0000
## OccDepr -9.41587 0.995338 1597 -9.459965 0.0000
## Correlation:
## (Intr) Occasn dprgrp mndfl_m mndfl_c
## Occasion 0.002
## deprgroup -0.284 -0.002
## mindful_mean -0.969 -0.002 0.108
## mindful_center 0.011 -0.436 0.004 0.004
## OccDepr -0.058 -0.001 0.321 0.000 0.002
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -2.732249099 -0.497489705 0.001225253 0.499338011 2.664816653
##
## Number of Observations: 2000
## Number of Groups: 400Anova(d.6.out_si)anova(d.6.out_si,d.6.out_ar, d.6.out_cs)########Comparing RANDOM covariance structures##########
#Variance components
sec.5a = lmer(posaff ~ Occasion + deprgroup +
mindful_mean + mindful_center + OccDepr +
(1+Occasion+mindful_center||ID), data=class07, REML=FALSE)
summary(sec.5a)## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
## method [lmerModLmerTest]
## Formula: posaff ~ Occasion + deprgroup + mindful_mean + mindful_center +
## OccDepr + (1 + Occasion + mindful_center || ID)
## Data: class07
##
## AIC BIC logLik deviance df.resid
## 16770.2 16826.2 -8375.1 16750.2 1990
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.70692 -0.50583 -0.00037 0.50346 2.63628
##
## Random effects:
## Groups Name Variance Std.Dev.
## ID (Intercept) 4484.3602 66.9654
## ID.1 Occasion 95.1265 9.7533
## ID.2 mindful_center 0.2478 0.4978
## Residual 35.7576 5.9798
## Number of obs: 2000, groups: ID, 400
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 6.2773 20.0667 402.1190 0.313 0.755
## Occasion 14.4235 0.5379 535.4529 26.815 < 2e-16 ***
## deprgroup -50.4500 7.0245 462.8084 -7.182 2.76e-12 ***
## mindful_mean 2.7431 0.6578 400.0051 4.170 3.74e-05 ***
## mindful_center 0.8376 0.1397 245.5754 5.994 7.25e-09 ***
## OccDepr -9.4438 0.9968 399.0828 -9.474 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) Occasn dprgrp mndfl_m mndfl_c
## Occasion -0.008
## deprgroup -0.284 0.000
## mindful_men -0.969 -0.001 0.110
## mindfl_cntr 0.017 -0.377 0.000 0.004
## OccDepr -0.049 0.001 0.272 0.000 -0.003relgrad <- with(sec.5a@optinfo$derivs,solve(Hessian,gradient))
max(abs(relgrad)) #note, greater than .001, so ugh## [1] 3.600084e-06#Unstructured
sec.5b = lmer(posaff ~ Occasion + deprgroup +
mindful_mean + mindful_center + OccDepr +
(1+Occasion+mindful_center|ID), data=class07, REML=FALSE)## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.0115346 (tol = 0.002, component 1)summary(sec.5b)## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
## method [lmerModLmerTest]
## Formula: posaff ~ Occasion + deprgroup + mindful_mean + mindful_center +
## OccDepr + (1 + Occasion + mindful_center | ID)
## Data: class07
##
## AIC BIC logLik deviance df.resid
## 16772.4 16845.2 -8373.2 16746.4 1987
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.73202 -0.49755 0.00122 0.49933 2.66459
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## ID (Intercept) 4445.4494 66.6742
## Occasion 100.9123 10.0455 0.06
## mindful_center 0.6229 0.7893 -0.09 -0.30
## Residual 35.2197 5.9346
## Number of obs: 2000, groups: ID, 400
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 6.0374 20.0298 402.8625 0.301 0.763
## Occasion 14.4318 0.5520 380.7675 26.142 < 2e-16 ***
## deprgroup -50.3343 7.1208 401.3017 -7.069 6.98e-12 ***
## mindful_mean 2.7488 0.6564 400.2353 4.188 3.47e-05 ***
## mindful_center 0.8350 0.1443 201.9896 5.787 2.71e-08 ***
## OccDepr -9.4158 0.9938 399.6565 -9.475 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) Occasn dprgrp mndfl_m mndfl_c
## Occasion 0.002
## deprgroup -0.284 -0.002
## mindful_men -0.969 -0.002 0.108
## mindfl_cntr 0.011 -0.436 0.004 0.004
## OccDepr -0.058 -0.001 0.321 0.000 0.002
## convergence code: 0
## Model failed to converge with max|grad| = 0.0115346 (tol = 0.002, component 1)relgrad <- with(sec.5b@optinfo$derivs,solve(Hessian,gradient))
max(abs(relgrad)) #note, greater than .001, so ugh## [1] 0.002825362anova(sec.5a,sec.5b)#For lmer objects, we can use stargazer and tab_model
#Stargazer makes a nicely formated model summary table (albeit incomplete)
class(sec.0) <- "lmerMod"
class(sec.1) <- "lmerMod"
class(sec.2) <- "lmerMod"
class(sec.3) <- "lmerMod"
class(sec.4) <- "lmerMod"
class(sec.5) <- "lmerMod"
library(stargazer)##
## Please cite as:## Hlavac, Marek (2018). stargazer: Well-Formatted Regression and Summary Statistics Tables.## R package version 5.2.2. https://CRAN.R-project.org/package=stargazerstargazer(sec.0, type = "text",
digits = 3,
star.cutoffs = c(0.05, 0.01, 0.001),
digit.separator = "") ##
## =================================================
## Dependent variable:
## -----------------------------
## posaff
## -------------------------------------------------
## Constant 88.749***
## (3.852)
##
## -------------------------------------------------
## Observations 2000
## Log Likelihood -10379.900
## Akaike Inf. Crit. 20765.790
## Bayesian Inf. Crit. 20782.600
## =================================================
## Note: *p<0.05; **p<0.01; ***p<0.001###sjPlot::tab_model is even better, because it formats the table in APA style and adds
# fixed and random effects and p-values
library(sjPlot)## Install package "strengejacke" from GitHub (`devtools::install_github("strengejacke/strengejacke")`) to load all sj-packages at once!tab_model(sec.0,sec.1,sec.2,sec.3,sec.4,sec.5)| posaff | posaff | posaff | posaff | posaff | posaff | |||||||||||||
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| Predictors | Estimates | CI | p | Estimates | CI | p | Estimates | CI | p | Estimates | CI | p | Estimates | CI | p | Estimates | CI | p |
| (Intercept) | 88.75 | 81.20 – 96.30 | <0.001 | 57.39 | 50.46 – 64.32 | <0.001 | 73.76 | 64.14 – 83.37 | <0.001 | -6.84 | -46.05 – 32.38 | 0.733 | -3.10 | -42.38 – 36.18 | 0.877 | 6.28 | -33.05 – 45.61 | 0.754 |
| Occasion |
|
|
|
15.68 | 14.60 – 16.76 | <0.001 | 15.68 | 14.60 – 16.76 | <0.001 | 15.68 | 14.60 – 16.76 | <0.001 | 14.43 | 13.28 – 15.59 | <0.001 | 14.42 | 13.37 – 15.48 | <0.001 |
| deprgroup |
|
|
|
|
|
|
-31.93 | -45.36 – -18.51 | <0.001 | -28.77 | -42.00 – -15.54 | <0.001 | -32.22 | -45.47 – -18.97 | <0.001 | -50.45 | -64.22 – -36.68 | <0.001 |
| mindful_mean |
|
|
|
|
|
|
|
|
|
2.73 | 1.44 – 4.01 | <0.001 | 2.74 | 1.45 – 4.03 | <0.001 | 2.74 | 1.45 – 4.03 | <0.001 |
| mindful_center |
|
|
|
|
|
|
|
|
|
|
|
|
0.83 | 0.56 – 1.11 | <0.001 | 0.84 | 0.56 – 1.11 | <0.001 |
| OccDepr |
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
-9.44 | -11.40 – -7.49 | <0.001 |
| Random Effects | ||||||||||||||||||
| σ2 | 946.70 | 37.26 | 37.26 | 37.26 | 35.72 | 35.76 | ||||||||||||
| τ00 | 5746.20 ID | 4976.91 ID | 4676.08 ID | 4482.37 ID | 4482.59 ID | 4484.36 ID | ||||||||||||
|
|
117.39 ID.1 | 95.13 ID.1 | ||||||||||||||||
|
|
0.27 ID.2 | 0.25 ID.2 | ||||||||||||||||
| τ11 | 117.93 ID.Occasion | 117.92 ID.Occasion | 117.93 ID.Occasion | |||||||||||||||
| ρ01 | 0.16 ID | 0.06 ID | 0.06 ID | |||||||||||||||
| ICC | 0.86 | 0.99 | 0.99 | 0.99 | 0.99 | 0.99 | ||||||||||||
| N | 400 ID | 400 ID | 400 ID | 400 ID | 400 ID | 400 ID | ||||||||||||
| Observations | 2000 | 2000 | 2000 | 2000 | 2000 | 2000 | ||||||||||||
| Marginal R2 / Conditional R2 | 0.000 / 0.859 | 0.074 / 0.994 | 0.118 / 0.994 | 0.148 / 0.994 | 0.182 / 0.994 | 0.243 / 0.994 | ||||||||||||
#Plotting our model-based estimates; let's do it by under-represented group
#Save predicted values to data frame for plotting
pred_sec.5<-(fitted(sec.5))
pred5<-as.data.frame(pred_sec.5)
class07<-cbind(class07,pred5)
#Here, we show fixed and random occasion effect by dep group
qplot(y=pred_sec.5, x=Occasion, colour = deprgroup,
data = class07)
#Here, we superimpose best fitting line by dep group
ggplot(class07) +
aes(x=class07$Occasion, y=class07$pred_sec.5, col=as.factor(class07$dep2)) + #
geom_point(size=0.1) +
stat_smooth(method="loess", se = FALSE)## Warning: Use of `class07$Occasion` is discouraged. Use `Occasion` instead.## Warning: Use of `class07$pred_sec.5` is discouraged. Use `pred_sec.5` instead.## Warning: Use of `class07$dep2` is discouraged. Use `dep2` instead.## Warning: Use of `class07$Occasion` is discouraged. Use `Occasion` instead.## Warning: Use of `class07$pred_sec.5` is discouraged. Use `pred_sec.5` instead.## Warning: Use of `class07$dep2` is discouraged. Use `dep2` instead.## `geom_smooth()` using formula 'y ~ x'
#Here, we vanish the data points in previous plot so we can see lines better; geom_point removed
ggplot(class07) +
aes(x=class07$Occasion, y=class07$pred_sec.5, col=as.factor(class07$dep2)) + #
stat_smooth(method="loess", se = FALSE)## Warning: Use of `class07$Occasion` is discouraged. Use `Occasion` instead.## Warning: Use of `class07$pred_sec.5` is discouraged. Use `pred_sec.5` instead.## Warning: Use of `class07$dep2` is discouraged. Use `dep2` instead.## `geom_smooth()` using formula 'y ~ x'
#With white background
finalclass07<-ggplot(class07) +
aes(x=class07$Occasion, y=class07$pred_sec.5, col=as.factor(class07$dep2)) + #
stat_smooth(method="loess", se = FALSE) +
theme_bw () +
labs(y="Model-estimated posaff", x = "Occasion of measurement")
print(finalclass07 + labs(colour = "dep group"))## Warning: Use of `class07$Occasion` is discouraged. Use `Occasion` instead.## Warning: Use of `class07$pred_sec.5` is discouraged. Use `pred_sec.5` instead.## Warning: Use of `class07$dep2` is discouraged. Use `dep2` instead.## `geom_smooth()` using formula 'y ~ x'
End