1 Data Management

pacman::p_load(knitr, simstudy, tidyverse, nlme, rms)

Read/Display Data

dta <- read.csv('Data/reading_g3.csv')

head(dta)
##   Group Post Pre  ID
## 1     C   41  32 S11
## 2     C   46  42 S12
## 3     C   48  45 S13
## 4     C   54  46 S14
## 5     C   52  49 S15
## 6     C   57  48 S16

Reshape data from wide to long

dtaL <- gather(data=dta, key=Time, value=Score, Pre:Post, factor_key=TRUE)  
kable(head(arrange(dtaL, ID), 10))
Group ID Time Score
C S11 Pre 32
C S11 Post 41
C S12 Pre 42
C S12 Post 46
C S13 Pre 45
C S13 Post 48
C S14 Pre 46
C S14 Post 54
C S15 Pre 49
C S15 Post 52

2 Visdualization

Generate a boxplot with Time (pre vs post) against Score(color different groups)

boxp <- ggplot(dtaL , 
               aes(x=Time, y=Score, col=Group))+ 
  geom_boxplot(outlier.size=0) + 
  geom_point(aes(fill=Group, col=NULL), 
             shape=21, alpha=0.5, size=2, 
             position=position_jitterdodge(jitter.width=0.2)) + 
  theme_bw() + 
  xlab("") 
boxp

Generate a regression line plot with pretest score against posttest score (color different groups)

ggplot(dta, 
       aes (x=Pre, y=Post, col=Group)) + 
  geom_point() + 
  geom_smooth(method="lm", se=FALSE) + 
  labs(x="Pre-test reading score",
       y="Post-test reading score")+
  theme_minimal()
## `geom_smooth()` using formula 'y ~ x'

# The constrained baseline model X is the model matrix comprises both main effect(time & group) and interaction effect(time x group)

X <- model.matrix(~ Time * Group, data = dtaL)
colnames(X)
## [1] "(Intercept)"      "TimePost"         "GroupT1"          "GroupT2"         
## [5] "TimePost:GroupT1" "TimePost:GroupT2"

Xalt is the model matrix excluding group information (GroupT1, GroupT2)

Xalt <- X[, c("TimePost", "TimePost:GroupT1", "TimePost:GroupT2")] 
colnames(Xalt)
## [1] "TimePost"         "TimePost:GroupT1" "TimePost:GroupT2"

3 Model specification

3.1 m0_gls

m0_gls is the model of Score ~ Xalt, with residual covariance matrix specificied as CS (weights= assumed equal variance at different levels of Time)

m0_gls <- gls(Score ~ Xalt, 
                weights = varIdent(form = ~ 1 | Time), 
                correlation=corSymm (form = ~ 1 | ID), 
                data = dtaL)
summary(m0_gls)
## Generalized least squares fit by REML
##   Model: Score ~ Xalt 
##   Data: dtaL 
##        AIC     BIC    logLik
##   350.2136 364.391 -168.1068
## 
## Correlation Structure: General
##  Formula: ~1 | ID 
##  Parameter estimate(s):
##  Correlation: 
##   1   
## 2 0.91
## Variance function:
##  Structure: Different standard deviations per stratum
##  Formula: ~1 | Time 
##  Parameter estimates:
##      Pre     Post 
## 1.000000 1.016413 
## 
## Coefficients:
##                         Value Std.Error  t-value p-value
## (Intercept)          50.50000 1.2296621 41.06819  0.0000
## XaltTimePost          6.38144 0.9027047  7.06924  0.0000
## XaltTimePost:GroupT1 -2.80337 1.2699097 -2.20753  0.0314
## XaltTimePost:GroupT2  7.15906 1.2699097  5.63745  0.0000
## 
##  Correlation: 
##                      (Intr) XltTmP XTP:GT1
## XaltTimePost         -0.102               
## XaltTimePost:GroupT1  0.000 -0.703        
## XaltTimePost:GroupT2  0.000 -0.703  0.500 
## 
## Standardized residuals:
##        Min         Q1        Med         Q3        Max 
## -2.7467890 -0.3836191  0.1374178  0.7475076  1.4493709 
## 
## Residual standard error: 6.735137 
## Degrees of freedom: 60 total; 56 residual

3.2 Generate prediction plot

Using Gls() function to plot Time against Estimated mean with 95 % CI

m0_Gls <- Gls(Score ~ Xalt, 
                weights = varIdent(form = ~ 1 | Time), 
                correlation=corSymm (form = ~ 1 | ID), 
                data = dtaL)
predictions <- cbind(dtaL, predict(m0_Gls, dtaL, conf.int=0.95))
predictions$Group2 <- factor(predictions$Group, levels=c("All", "C", "T1", "T2"))
predictions$Group2[predictions$Time=="Pre"] <- "All"
pd <- position_dodge(.08)
limits <- aes(ymax=upper , ymin=lower, shape=Group2)
pCI <- ggplot(predictions, 
              aes(y=linear.predictors, 
                  x=Time)) + 
  geom_errorbar(limits, 
                width=0.09, 
                position=pd) +
  geom_line(aes(group=Group, 
                y=linear.predictors), 
            position=pd)+
  geom_point(aes(fill=Group2),
             position=pd, 
             shape=21, 
             size=rel(4)) + 
  scale_fill_manual(values=c("white", "gray80", "gray20", "black")) + 
  theme_minimal() + 
  theme(panel.grid.major=element_blank(), 
        legend.title=element_blank(), 
        text=element_text(size=14), 
        legend.position="bottom") + 
  labs(x="Time", 
       y="Estimated mean with corresponding 95% CI")
## Warning: Ignoring unknown aesthetics: shape
pCI

The End