I am trying to make sure I understand the use of multilevel modelling to look at change (on questionnaire scores) across sessions of psychotherapy. This is widely done in my applied field. I am a researcher (and ex-therapist) not a statistician and have no direct availability of a statistician but I'm fairly savvy about stats. I am simulating, rather crudely, data which have both fixed and random intercepts and slopes with a quadratic as well as a linear effect of session number. I'm using lmer from lme4 overloaded by lmerTest in 64 bit R 3.4.1 on Windoze running in Rstudio.
The thing that's baffling me is that centering my predictor (session count) seems to have NO effect on my results, neither on the fixed nor the random effects. I would have expected the intercept to have changed at the very least given the effects I've simulated.
Details: session count can range from 1 to 51 (the number of sessions per client in my simulation can vary, realistically, from 4 to 51 so no client has fewer than four). I have modelled in a random linear and quadratic effect of session count but such that the overall effects are both negative. Code and session info are at the end of the post. Here is code and output from uncentred:
m4 <- lmer(DEP_TOTAL ~ poly(session,2)+(1+poly(session,2)|validPOC_ID),
data=longdat4)
summary(m4)
## Linear mixed model fit by REML t-tests use Satterthwaite approximations
## to degrees of freedom [lmerMod]
## Formula:
## DEP_TOTAL ~ poly(session, 2) + (1 + poly(session, 2) | validPOC_ID)
## Data: longdat4
##
## REML criterion at convergence: -29507
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.6274 -0.6583 0.0035 0.6420 4.0849
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## validPOC_ID (Intercept) 9.309e+00 3.05101
## poly(session, 2)1 7.390e+03 85.96686 0.20
## poly(session, 2)2 4.154e+02 20.38232 0.21 1.00
## Residual 4.973e-03 0.07052
## Number of obs: 15399, groups: validPOC_ID, 500
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 8.6192 0.1367 526.0000 63.07 <2e-16 ***
## poly(session, 2)1 -179.5508 3.9684 4712.0000 -45.25 <2e-16 ***
## poly(session, 2)2 -42.1333 0.9497 4029.0000 -44.37 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) p(,2)1
## ply(sss,2)1 0.210
## ply(sss,2)2 0.215 0.995
Here is the centered one (lengthCentre is 29.5).
longdat4 $ sessionC1 <- longdat4 $ session - lengthCentre summary ( longdat4 $ sessionC1 )
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -28.50 -21.50 -12.50 -10.09 -0.50 21.50
m5 <- lmer ( DEP_TOTAL ~ poly ( sessionC1 , 2 ) + ( 1 + poly ( sessionC1 , 2 ) | validPOC_ID ) , data = longdat4 ) summary ( m5 )
## Linear mixed model fit by REML t-tests use Satterthwaite approximations
## to degrees of freedom [lmerMod]
## Formula:
## DEP_TOTAL ~ poly(sessionC1, 2) + (1 + poly(sessionC1, 2) | validPOC_ID)
## Data: longdat4
##
## REML criterion at convergence: -29507
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.6274 -0.6583 0.0035 0.6420 4.0849
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## validPOC_ID (Intercept) 9.309e+00 3.05101
## poly(sessionC1, 2)1 7.390e+03 85.96686 0.20
## poly(sessionC1, 2)2 4.154e+02 20.38232 0.21 1.00
## Residual 4.973e-03 0.07052
## Number of obs: 15399, groups: validPOC_ID, 500
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 8.6192 0.1367 526.0000 63.07 <2e-16 ***
## poly(sessionC1, 2)1 -179.5508 3.9684 4712.0000 -45.25 <2e-16 ***
## poly(sessionC1, 2)2 -42.1333 0.9497 4029.0000 -44.37 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) p(C1,2)1
## ply(sC1,2)1 0.210
## ply(sC1,2)2 0.215 0.995
This makes no sense to me. I had expected that the very crude centering must shift the fixed intercept as the x=0 has shifted left by 29.5 and I have strong linear and quadratic effects in the model. I also thought (the whole reason I'm doing this) that it would change the correlation between the linear and quadratic terms as I should be moving from a strong positive correlation between them to a much smaller one.
I'm sure I'm doing something utterly stupid, or revealing shocking failure to comprehend the multilevel model, but I'm looking straight through the explanation and stuck. I hope someone will help me!
TIA,
Chris
Here is entire pertinent code.
```{r sim4,warning=FALSE,message=FALSE,cache=TRUE}
require(reshape2)
set.seed(12345)
### set up the parameters for the simulation
n <- 500
minNSessions <- 4
maxNSessions <- 55
linSlope <- -.02
errVar <- .005
### derived variables
lengthRange <- minNSessions:maxNSessions
lengthCentre <- (minNSessions+maxNSessions)/2
ID <- 1:n
### create a variable, length = n for the session counts
sessions <- sample(lengthRange,n,replace=TRUE)
### for now, give all starting scores of 10
startScores <- rnorm(n,mean=10,sd=3)
slopes <- rnorm(n,mean=linSlope,sd=.01) # random slopes
### set up matrix to store the simulated scores
scores <- matrix(rep(NA,n*maxNSessions),nrow=n)
dat4 <- as.data.frame(cbind(ID,sessions,startScores,slopes,scores))
for (i in 1:nrow(dat4)) {
nSessions <- dat4$sessions[i]
tmpScores <- rep(dat4$startScores[i],nSessions) # scores: no slopes
tmpScores <- tmpScores + c(0,slopes[i]*(2:nSessions)) # scores: add linear effect
tmpScores <- tmpScores + c(0,slopes[i]*((2:nSessions)^2))/10 # add quadratic effect
tmpScores <- tmpScores + rnorm(nSessions,mean=0,sd=sqrt(errVar)) # now add error
dat4[i,5:(4+nSessions)] <- tmpScores # put into data frame
}
### now melt to create a long data frame
longdat4 <- melt(dat4,id.vars="ID",measure.vars = 4:55)
colnames(longdat4) <- c("validPOC_ID","wks","DEP_TOTAL")
longdat4$session <- as.numeric(substr(longdat4$wks,2,3))-4 # minus 4 to get start at 1
# as data start in column 5
longdat4 <- longdat4[order(longdat4$validPOC_ID,longdat4$wks),]
longdat4 <- na.omit(longdat4) # ditch the empty data
```
So the data frame longdat4 has
`r nrow(longdat4)` rows and
`r ncol(longdat4)` columns and the column names are
`r colnames(longdat4)`. Here is a bit of the head of the wide data frame (dat4) and the head of the long one (longdat4).
```{r head4,warning=FALSE,message=FALSE,cache=TRUE}
require(knitr)
kable(head(dat4[,1:10]),digits=3)
kable(head(longdat4))
```
Elisa's lattice plot.
```{r spaghetti4,warning=FALSE,message=FALSE,cache=TRUE}
library(lattice)
samp <- sample(unique(longdat4$validPOC_ID), 35)
xyplot(DEP_TOTAL ~ wks|validPOC_ID, longdat4,
type = c("p","g","r"),
subset = (validPOC_ID %in% samp),
col="dark blue",
col.line="black",
strip = F,
xlab = "Weeks in treatment",
ylab = "Depression Scores")
```
Now Elisa's tangle plot. Just first 40 simulated cases.
```{r tangle4, warning=FALSE,message=FALSE,cache=TRUE}
library(ggplot2)
first40 <- longdat4[longdat4$validPOC_ID <=40,]
tanglePlot <- ggplot(data = longdat4, aes(x = wks, y = DEP_TOTAL)) +
geom_line(data = longdat4, alpha = 0.5, aes(group = validPOC_ID)) +
theme_bw()
tanglePlot
```
```{r model4, warning=FALSE, message=FALSE,cache=TRUE}
require(lmerTest)
### this next fails to converge and has perfect correlation of linear & quadratic effects
# m4 <- lmer(DEP_TOTAL ~ session+I(session^2)+(1+session+I(session^2)|validPOC_ID),
# data=longdat4)
### is this better:
summary(longdat4$session)
hist(longdat4$session,
main="Histogram of raw session counts",
xlab="Raw session n",
col="grey")
m4 <- lmer(DEP_TOTAL ~ poly(session,2)+(1+poly(session,2)|validPOC_ID),
data=longdat4)
summary(m4)
```
This next should test the random effect (and should find it highly significant I think).
```{r mode4rand,eval=FALSE}
rand(m4)
### fails with message:
### Error in terms.formula(tmp, simplify = TRUE) :
### invalid model formula in ExtractVars
```
And this summarises the random effects. These are deviations from the fixed effects so both are centred on zero.
```{r coeffs4,fig.width=11,fig.height=8,cache=TRUE}
apply(ranef(m4)$validPOC_ID,2,summary)
par(mfrow=c(1,3))
hist(ranef(m4)$validPOC_ID[,1],
xlab="Intercept deviation values",
main="Random intercepts",
col="grey")
abline(v=0)
hist(ranef(m4)$validPOC_ID[,2],
main="Random slopes",
xlab="Linear slope deviation values",
col="grey")
abline(v=0)
hist(ranef(m4)$validPOC_ID[,3],
main="Random slopes",
xlab="Quadratic slope deviation values",
col="grey")
abline(v=0)
```
# Simulation 5: random intercept and random slope with quadratic effect, CENTRED
OK. Random intercept and slope has linear and quadratic effects. Now going to centre the session count by subtracting lengthCentre
```{r model5, warning=FALSE, message=FALSE}
longdat4$sessionC1 <- longdat4$session - lengthCentre
summary(longdat4$sessionC1)
hist(longdat4$sessionC1,
main="Histogram of centred session counts",
xlab="Centred session n",
col="grey")
abline(v=0)
require(lmerTest)
m5 <- lmer(DEP_TOTAL ~ poly(sessionC1,2)+(1+poly(sessionC1,2)|validPOC_ID),
data=longdat4)
summary(m5)
```
```{r mode5rand,eval=FALSE}
rand(m5)
### fails with message:
### Error in terms.formula(tmp, simplify = TRUE) :
### invalid model formula in ExtractVars
```
```{r coeffs5,fig.width=11,fig.height=8}
apply(ranef(m5)$validPOC_ID,2,summary)
par(mfrow=c(1,3))
hist(ranef(m5)$validPOC_ID[,1],
xlab="Intercept deviation values",
main="Random intercepts",
col="grey")
abline(v=0)
hist(ranef(m5)$validPOC_ID[,2],
main="Random slopes",
xlab="Linear slope deviation values",
col="grey")
abline(v=0)
hist(ranef(m5)$validPOC_ID[,3],
main="Random slopes",
xlab="Quadratic slope deviation values",
col="grey")
abline(v=0)
```
And here is sessionInfo()