library(arm)
## Loading required package: MASS
## Loading required package: Matrix
## Loading required package: lme4
## Loading required package: Rcpp
##
## arm (Version 1.7-03, built: 2014-4-27)
##
## Working directory is /Users/levibrackman/Google Drive/R/Emily Griffith
library(lmerTest)
## KernSmooth 2.23 loaded
## Copyright M. P. Wand 1997-2009
##
## Attaching package: 'lmerTest'
##
## The following object is masked from 'package:lme4':
##
## lmer
##
## The following object is masked from 'package:stats':
##
## step
library(psych)
##
## Attaching package: 'psych'
##
## The following objects are masked from 'package:arm':
##
## logit, rescale, sim
data <- read.csv("EmilyGriffith_all.csv")
data$ID <- data$Q1
data <- data[data$Time > 1, ]
# Create scale scores
data$meanCPS <- apply(data[, c("CPS_2", "CPS_3", "CPS_4", "CPS_10", "CPS_8",
"CPS_5", "CPS_6", "CPS_7", "CPS_1", "CPS_9")], 1, mean, na.rm = TRUE)
# Means or plotting
data$baseline <- ifelse(data$Time < 4, 0, 1)
pdata <- tapply(data[, "meanCPS"], data[, 3], mean, na.rm = TRUE)
plot(pdata, type = "l")
M0 <- lmer(meanCPS ~ 1 + (1 | ID), data = data)
fixef(M0)
## (Intercept)
## 6.293
confint(M0)
## Computing profile confidence intervals ...
## 2.5 % 97.5 %
## .sig01 0.0000 0.6695
## .sigma 0.5352 0.8368
## (Intercept) 6.0991 6.4819
M1 <- update(M0, . ~ . + Time, REML = FALSE)
fixef(M1)
## (Intercept) Time
## 5.8762 0.1625
confint(M1)
## Computing profile confidence intervals ...
## Warning: convergence code 3 from bobyqa: bobyqa -- a trust region step
## failed to reduce q
## 2.5 % 97.5 %
## .sig01 0.0000 0.6450
## .sigma 0.5351 0.8404
## (Intercept) 5.3065 6.4279
## Time -0.0394 0.3734
M2 <- update(M1, . ~ . + baseline)
fixef(M2)
## (Intercept) Time baseline
## 6.2129 0.0114 0.3869
confint(M2)
## Computing profile confidence intervals ...
## Warning: convergence code 3 from bobyqa: bobyqa -- a trust region step
## failed to reduce q
## 2.5 % 97.5 %
## .sig01 0.0000 0.6506
## .sigma 0.5265 0.8285
## (Intercept) 5.4051 6.9998
## Time -0.3119 0.3447
## baseline -0.2796 1.0517
M3 <- update(M2, . ~ . + I(Time^2))
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
fixef(M3)
## (Intercept) Time baseline
## 6.2129 0.0114 0.3869
confint(M3)
## Computing profile confidence intervals ...
## Warning: convergence code 3 from bobyqa: bobyqa -- a trust region step
## failed to reduce q
## 2.5 % 97.5 %
## .sig01 0.0000 0.6506
## .sigma 0.5265 0.8285
## (Intercept) 5.4051 6.9998
## Time -0.3119 0.3447
## baseline -0.2796 1.0517
Pdata <- tapply(data[, "meanCPS"], data[, 3], mean, na.rm = TRUE)
# Add random noise to time to better see the points of interest
data$TimeJIT <- data$Time + runif(126, min = -0.1, max = 0.1)
## Warning: longer object length is not a multiple of shorter object length
## Error: replacement has 126 rows, data has 89
with(data, plot(TimeJIT, meanCPS, col = "grey", pch = "*"))
## Error: object 'TimeJIT' not found
lines(pdata, col = "red", lwd = 2)