Loading the dataset

data.test4 <- read.csv("~/Final_Adult_Study_R_Docs/adult_study011615.csv")
# Load the psych package
library(psych)
items <- grep("LOTR[0-9]", names(data.test4), value=TRUE)
scaleKey <- c(1, 1, -1, 1, 1, 1, -1, 1, -1, 1)
data.test4[,items] <- apply(data.test4[,items], 2, as.numeric)
data.test4$meanLOTR <- scoreItems(scaleKey, items = data.test4[, items])$score
library(reshape2); library(car); library(Amelia);library(mitools);library(nlme);library(predictmeans)
## 
## Attaching package: 'car'
## 
## The following object is masked from 'package:psych':
## 
##     logit
## 
## Loading required package: Rcpp
## ## 
## ## Amelia II: Multiple Imputation
## ## (Version 1.7.3, built: 2014-11-14)
## ## Copyright (C) 2005-2015 James Honaker, Gary King and Matthew Blackwell
## ## Refer to http://gking.harvard.edu/amelia/ for more information
## ## 
## Loading required package: lme4
## Loading required package: Matrix
## 
## Attaching package: 'lme4'
## 
## The following object is masked from 'package:nlme':
## 
##     lmList
#Remove the meanLOTR and ID Group and wave from data.test4 and create a new #dataset with only those variables.
data <- data.test4[,c("ID", "GROUP", "wave", "meanLOTR")]
#Use dcast to cnage from long-format data to wide format data
data <- dcast(data, ID + GROUP ~ wave, mean, value.var = "meanLOTR")
# Changing all NaNs to NA
data[,3:5] <- apply(data[,3:5],2,function(x) recode(x, "NaN = NA") )

Unsing the mapply function we create a new data set with ID Group baseline meanLOTR and wave 2 and 3 of meanLOTR. So we have a Baseline, which is Time 1 (placed in column 3 one on top of the other) to compare to both Time 2 and 3 (placed in column 4 one on top of the other). In the next line of code we then create a separate column called “wave” which calls the first 89 (which compares Time 2 to Baseline) “wave 1” and then the second 89 we call “wave 2” which compares Time 3 to Baseline. In the third line of code we add names to the new columns of the dataset.

data2 <- as.data.frame(mapply(c,data[,1:4], data[,c(1:3,5)]))
data2$wave <- rep(1:2, each=89)
names(data2) <- c("ID", "GROUP", "BASELINE", "meanLOTR", "WAVE")

Intention to treat model (ITT) where we keep the cases who dropped out and did not complete the study (http://en.wikipedia.org/wiki/Intention-to-treat_analysis).

data2[which(data2$GROUP ==2), "GROUP"] <- 1

Make GROUP and ID a factor

data2$GROUP <-as.factor(data2$GROUP)
data2$ID <-as.factor(data2$ID)

Imputing missing data

MI <- amelia(data2, 50, idvars = c("ID"), ords = "GROUP")
## -- Imputation 1 --
## 
##   1  2  3  4  5  6  7
## 
## -- Imputation 2 --
## 
##   1  2  3  4  5  6  7  8  9 10 11 12
## 
## -- Imputation 3 --
## 
##   1  2  3  4  5  6  7  8
## 
## -- Imputation 4 --
## 
##   1  2  3  4  5  6  7  8  9
## 
## -- Imputation 5 --
## 
##   1  2  3  4  5  6  7  8
## 
## -- Imputation 6 --
## 
##   1  2  3  4  5  6  7  8  9 10 11
## 
## -- Imputation 7 --
## 
##   1  2  3  4  5  6  7  8  9 10 11 12
## 
## -- Imputation 8 --
## 
##   1  2  3  4  5  6  7  8  9 10
## 
## -- Imputation 9 --
## 
##   1  2  3  4  5  6  7  8  9 10
## 
## -- Imputation 10 --
## 
##   1  2  3  4  5  6  7  8  9 10 11
## 
## -- Imputation 11 --
## 
##   1  2  3  4  5  6  7
## 
## -- Imputation 12 --
## 
##   1  2  3  4  5  6  7  8
## 
## -- Imputation 13 --
## 
##   1  2  3  4  5  6  7  8  9
## 
## -- Imputation 14 --
## 
##   1  2  3  4  5  6  7  8  9
## 
## -- Imputation 15 --
## 
##   1  2  3  4  5  6  7
## 
## -- Imputation 16 --
## 
##   1  2  3  4  5  6  7  8
## 
## -- Imputation 17 --
## 
##   1  2  3  4  5  6  7  8  9
## 
## -- Imputation 18 --
## 
##   1  2  3  4  5  6
## 
## -- Imputation 19 --
## 
##   1  2  3  4  5  6  7  8
## 
## -- Imputation 20 --
## 
##   1  2  3  4  5  6  7  8
## 
## -- Imputation 21 --
## 
##   1  2  3  4  5  6  7  8  9 10 11 12
## 
## -- Imputation 22 --
## 
##   1  2  3  4  5  6  7  8  9 10 11 12
## 
## -- Imputation 23 --
## 
##   1  2  3  4  5  6  7  8  9 10 11
## 
## -- Imputation 24 --
## 
##   1  2  3  4  5  6  7  8  9 10
## 
## -- Imputation 25 --
## 
##   1  2  3  4  5  6  7  8  9
## 
## -- Imputation 26 --
## 
##   1  2  3  4  5  6
## 
## -- Imputation 27 --
## 
##   1  2  3  4  5  6  7  8  9 10 11
## 
## -- Imputation 28 --
## 
##   1  2  3  4  5  6  7  8
## 
## -- Imputation 29 --
## 
##   1  2  3  4  5  6  7  8  9 10 11
## 
## -- Imputation 30 --
## 
##   1  2  3  4  5  6
## 
## -- Imputation 31 --
## 
##   1  2  3  4  5  6  7  8  9 10 11
## 
## -- Imputation 32 --
## 
##   1  2  3  4  5  6  7  8  9
## 
## -- Imputation 33 --
## 
##   1  2  3  4  5  6  7  8
## 
## -- Imputation 34 --
## 
##   1  2  3  4  5  6
## 
## -- Imputation 35 --
## 
##   1  2  3  4  5  6  7  8  9
## 
## -- Imputation 36 --
## 
##   1  2  3  4  5  6  7  8  9 10 11
## 
## -- Imputation 37 --
## 
##   1  2  3  4  5  6
## 
## -- Imputation 38 --
## 
##   1  2  3  4  5  6  7  8
## 
## -- Imputation 39 --
## 
##   1  2  3  4  5  6  7  8  9 10
## 
## -- Imputation 40 --
## 
##   1  2  3  4  5  6  7  8  9 10
## 
## -- Imputation 41 --
## 
##   1  2  3  4  5  6  7  8  9 10
## 
## -- Imputation 42 --
## 
##   1  2  3  4  5  6  7
## 
## -- Imputation 43 --
## 
##   1  2  3  4  5
## 
## -- Imputation 44 --
## 
##   1  2  3  4  5  6  7  8  9
## 
## -- Imputation 45 --
## 
##   1  2  3  4  5  6  7  8
## 
## -- Imputation 46 --
## 
##   1  2  3  4  5  6  7  8  9
## 
## -- Imputation 47 --
## 
##   1  2  3  4  5  6  7  8
## 
## -- Imputation 48 --
## 
##   1  2  3  4  5  6  7
## 
## -- Imputation 49 --
## 
##   1  2  3  4  5  6
## 
## -- Imputation 50 --
## 
##   1  2  3  4  5  6  7

Creating new dataset with missing data imputed

data(MI$imputations)
## Warning in data(MI$imputations): data set 'MI$imputations' not found
allimplogreg<-lapply(MI$imputations,function(X) {lme(meanLOTR ~ GROUP * WAVE + BASELINE, random = ~1 | ID, data = X, method = "ML", na.action = "na.omit")})
betas<-MIextract(allimplogreg, fun=fixef)
vars<-MIextract(allimplogreg, fun=vcov)
summary(MIcombine(betas,vars))
## Multiple imputation results:
##       MIcombine.default(betas, vars)
##                results         se     (lower    upper) missInfo
## (Intercept)  1.1703617 0.31597247  0.5492315 1.7914919     35 %
## GROUP1      -0.0131131 0.31319407 -0.6278620 0.6016358     25 %
## WAVE         0.0610420 0.15473346 -0.2433546 0.3654386     39 %
## BASELINE     0.4036809 0.08772351  0.2310290 0.5763328     41 %
## GROUP1:WAVE -0.1402872 0.19949995 -0.5320922 0.2515178     29 %

Check results with Imputations using Zelig

library("Zelig")
## Loading required package: boot
## 
## Attaching package: 'boot'
## 
## The following object is masked from 'package:car':
## 
##     logit
## 
## The following object is masked from 'package:psych':
## 
##     logit
## 
## Loading required package: MASS
## Loading required package: sandwich
## ZELIG (Versions 4.2-1, built: 2013-09-12)
## 
## +----------------------------------------------------------------+
## |  Please refer to http://gking.harvard.edu/zelig for full       |
## |  documentation or help.zelig() for help with commands and      |
## |  models support by Zelig.                                      |
## |                                                                |
## |  Zelig project citations:                                      |
## |    Kosuke Imai, Gary King, and Olivia Lau.  (2009).            |
## |    ``Zelig: Everyone's Statistical Software,''                 |
## |    http://gking.harvard.edu/zelig                              |
## |   and                                                          |
## |    Kosuke Imai, Gary King, and Olivia Lau. (2008).             |
## |    ``Toward A Common Framework for Statistical Analysis        |
## |    and Development,'' Journal of Computational and             |
## |    Graphical Statistics, Vol. 17, No. 4 (December)             |
## |    pp. 892-913.                                                |
## |                                                                |
## |   To cite individual Zelig models, please use the citation     |
## |   format printed with each model run and in the documentation. |
## +----------------------------------------------------------------+
## 
## 
## 
## Attaching package: 'Zelig'
## 
## The following objects are masked from 'package:psych':
## 
##     alpha, describe, sim
## 
## The following object is masked from 'package:utils':
## 
##     cite
zelig.fit <- zelig(meanLOTR ~ GROUP * WAVE + BASELINE, random = ~1 | ID, data = MI$imputations,  model = "ls", cite = FALSE)
summary(zelig.fit)
## 
##   Model: ls
##   Number of multiply imputed data sets: 50 
## 
## Combined results:
## 
## Call:
## lm(formula = formula, weights = weights, model = F, data = data)
## 
## Coefficients:
##                  Value Std. Error      t-stat      p-value
## (Intercept)  1.1674750 0.32250563  3.62001435 3.279572e-04
## GROUP1      -0.0130177 0.33645831 -0.03869038 9.691443e-01
## WAVE         0.0610420 0.16649774  0.36662358 7.140764e-01
## BASELINE     0.4048653 0.08433436  4.80071610 2.718064e-06
## GROUP1:WAVE -0.1402987 0.21706050 -0.64635768 5.182247e-01
## 
## For combined results from datasets i to j, use summary(x, subset = i:j).
## For separate results, use print(summary(x), subset = i:j).

Check assumptions with Random Computations

data1=MI$imputations[[1]]
library("Zelig")
zelig.fitdata1 <- zelig(meanLOTR ~ GROUP * WAVE + BASELINE, random = ~1 | ID, data = data1,  model = "ls", cite = FALSE)
summary(zelig.fitdata1)
## 
## Call:
## lm(formula = formula, weights = weights, model = F, data = data)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.44274 -0.42456 -0.01929  0.31641  1.95271 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  1.11354    0.28082   3.965 0.000107 ***
## GROUP1      -0.03583    0.31775  -0.113 0.910354    
## WAVE         0.06619    0.14446   0.458 0.647380    
## BASELINE     0.45150    0.06702   6.737 2.31e-10 ***
## GROUP1:WAVE -0.13629    0.20094  -0.678 0.498528    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.6698 on 173 degrees of freedom
## Multiple R-squared:  0.2361, Adjusted R-squared:  0.2184 
## F-statistic: 13.37 on 4 and 173 DF,  p-value: 1.641e-09

Describe the meanLOTR variable by the GROUP variable

describeBy(data1[,3:4], group = data1$GROUP)
## group: 0
##          vars  n mean   sd median trimmed  mad  min max range skew
## BASELINE    1 86 2.44 0.79    2.5    2.39 0.89 1.10 4.3  3.20 0.43
## meanLOTR    2 86 2.31 0.86    2.2    2.26 0.85 0.65 4.2  3.55 0.51
##          kurtosis   se
## BASELINE    -0.47 0.09
## meanLOTR    -0.61 0.09
## -------------------------------------------------------- 
## group: 1
##          vars  n mean   sd median trimmed  mad  min max range skew
## BASELINE    1 92 2.36 0.71   2.30    2.33 0.82 1.10   4  2.90 0.40
## meanLOTR    2 92 2.04 0.63   2.06    2.01 0.69 0.71   4  3.29 0.41
##          kurtosis   se
## BASELINE    -0.64 0.07
## meanLOTR     0.09 0.07

Create a plot that visualizes meanLOTR variable by the GROUP variable

library(ggplot2)
## 
## Attaching package: 'ggplot2'
## 
## The following object is masked from 'package:psych':
## 
##     %+%
library(influence.ME)
## 
## Attaching package: 'influence.ME'
## 
## The following object is masked from 'package:stats':
## 
##     influence

Take a look at the residuals

residual <- lm(meanLOTR ~ BASELINE, data=data1)$residual

Plot the residuals to see that they are random

plot(density(residual))# A density plot

qqnorm(residual) # A quantile normal plot to checking normality
qqline(residual)

Checking the different between intervention and control groups residuals. This allows us to control for individual unsystematic differences.

data2$residual <- NA
sel1 <- which(!is.na(data1$meanLOTR)) 
sel2 <- which(!is.na(data1$BASELINE))
data1$residual[intersect(sel1,sel2)] <- residual
qplot(GROUP, meanLOTR, data=data1, geom="boxplot")

Plot of the difference between intervention and control groups.

qplot(GROUP, residual, data=data1, geom="boxplot")

Two way repeated measures ======================================================== Graphing the Two-Way Interaction. Both meanLOTR and the Residuals

# Load the nlme package
library(nlme)
with(data1, boxplot(meanLOTR ~ WAVE + GROUP))

with(data1, boxplot(residual ~ WAVE + GROUP))
Linear Mixed-Effects Model

Comparing Basline to Wave 2 and 3 by Group.

fullModeldata1 <- lme(meanLOTR ~ GROUP * WAVE + BASELINE, random = ~1 | ID, data = data1, method = "ML", na.action = "na.omit")

Cooks Distence

CookD(fullModeldata1)

Plot Cook’s distance:

plot(fullModeldata1, which="cook")
Check results on this random Imputation model
Results

Explanation of significance:

We asses the significance of our models by comparing them from the baseline model using the anova() function.
(Intercept): Where everything is 0
GROUP1: Is there a difference between group. If it is significant than there is a difference and the treatment had an effect.
WAVE: Asseses whether the effects gets bigger beteen time 2 and 3 (does not have to be significant)
BASELINE: Should not be significant. If it is then it shows that there is a difference between groups before the treatment.
GROUP1:WAVE: If this is significant then it means that the effect was either fleeting or it happened after the treatment i.e. between time 2 and 3.

summary(fullModeldata1)
## Linear mixed-effects model fit by maximum likelihood
##  Data: data1 
##        AIC      BIC   logLik
##   364.3479 386.6204 -175.174
## 
## Random effects:
##  Formula: ~1 | ID
##         (Intercept)  Residual
## StdDev:   0.3471795 0.5617483
## 
## Fixed effects: meanLOTR ~ GROUP * WAVE + BASELINE 
##                  Value  Std.Error DF   t-value p-value
## (Intercept)  1.1151834 0.27297711 87  4.085263  0.0001
## GROUP1      -0.0358719 0.28044457 87 -0.127911  0.8985
## WAVE         0.0661940 0.12288799 86  0.538653  0.5915
## BASELINE     0.4508313 0.07552153 86  5.969573  0.0000
## GROUP1:WAVE -0.1362951 0.17093472 86 -0.797352  0.4274
##  Correlation: 
##             (Intr) GROUP1 WAVE   BASELI
## GROUP1      -0.543                     
## WAVE        -0.675  0.657              
## BASELINE    -0.674  0.017  0.000       
## GROUP1:WAVE  0.482 -0.914 -0.719  0.005
## 
## Standardized Within-Group Residuals:
##         Min          Q1         Med          Q3         Max 
## -2.05293020 -0.50440966 -0.06844674  0.43782803  2.57390169 
## 
## Number of Observations: 178
## Number of Groups: 89

Another random imputation

data10=MI$imputations[[10]]
library("Zelig")
zelig.fitdata10 <- zelig(meanLOTR ~ GROUP * WAVE + BASELINE, random = ~1 | ID, data = data10,  model = "ls", cite = FALSE)
summary(zelig.fitdata10)
## 
## Call:
## lm(formula = formula, weights = weights, model = F, data = data)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.53479 -0.41756 -0.01204  0.34162  1.92129 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  1.38326    0.26500   5.220 5.09e-07 ***
## GROUP1      -0.12221    0.29955  -0.408    0.684    
## WAVE        -0.04923    0.13620  -0.361    0.718    
## BASELINE     0.37922    0.06337   5.984 1.22e-08 ***
## GROUP1:WAVE -0.13364    0.18945  -0.705    0.482    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.6315 on 173 degrees of freedom
## Multiple R-squared:  0.2272, Adjusted R-squared:  0.2093 
## F-statistic: 12.72 on 4 and 173 DF,  p-value: 4.29e-09

Describe the meanLOTR variable by the GROUP variable

describeBy(data10[,3:4], group = data10$GROUP)
## group: 0
##          vars  n mean   sd median trimmed  mad  min max range skew
## BASELINE    1 86 2.44 0.79   2.50    2.39 0.89 1.10 4.3  3.20 0.43
## meanLOTR    2 86 2.23 0.76   2.18    2.18 0.82 0.55 4.2  3.65 0.49
##          kurtosis   se
## BASELINE    -0.47 0.09
## meanLOTR    -0.27 0.08
## -------------------------------------------------------- 
## group: 1
##          vars  n mean   sd median trimmed  mad  min max range skew
## BASELINE    1 92 2.39 0.71    2.3    2.37 0.74 1.10   4  2.90 0.35
## meanLOTR    2 92 1.89 0.62    1.8    1.86 0.59 0.59   4  3.41 0.63
##          kurtosis   se
## BASELINE    -0.68 0.07
## meanLOTR     0.51 0.06

Create a plot that visualizes meanLOTR variable by the GROUP variable

library(ggplot2)
library(influence.ME)

Take a look at the residuals

residual <- lm(meanLOTR ~ BASELINE, data=data10)$residual

Plot the residuals to see that they are random

plot(density(residual))# A density plot

qqnorm(residual) # A quantile normal plot to checking normality
qqline(residual)

Checking the different between intervention and control groups residuals. This allows us to control for individual unsystematic differences.

data2$residual <- NA
sel1 <- which(!is.na(data10$meanLOTR)) 
sel2 <- which(!is.na(data10$BASELINE))
data10$residual[intersect(sel1,sel2)] <- residual
qplot(GROUP, meanLOTR, data=data10, geom="boxplot")

Plot of the difference between intervention and control groups.

qplot(GROUP, residual, data=data10, geom="boxplot")

Two way repeated measures ======================================================== Graphing the Two-Way Interaction. Both meanLOTR and the Residuals

# Load the nlme package
library(nlme)
with(data10, boxplot(meanLOTR ~ WAVE + GROUP))

with(data10, boxplot(residual ~ WAVE + GROUP))
Linear Mixed-Effects Model

Comparing Basline to Wave 2 and 3 by Group.

fullModeldata10 <- lme(meanLOTR ~ GROUP * WAVE + BASELINE, random = ~1 | ID, data = data10, method = "ML", na.action = "na.omit")

Cooks Distence

CookD(fullModeldata10)

Plot Cook’s distance:

plot(fullModeldata10, which="cook")
Check results on this random Imputation model
Results

Explanation of significance:

We asses the significance of our models by comparing them from the baseline model using the anova() function.
(Intercept): Where everything is 0
GROUP1: Is there a difference between group. If it is significant than there is a difference and the treatment had an effect.
WAVE: Asseses whether the effects gets bigger beteen time 2 and 3 (does not have to be significant)
BASELINE: Should not be significant. If it is then it shows that there is a difference between groups before the treatment.
GROUP1:WAVE: If this is significant then it means that the effect was either fleeting or it happened after the treatment i.e. between time 2 and 3.

summary(fullModeldata10)
## Linear mixed-effects model fit by maximum likelihood
##  Data: data10 
##        AIC      BIC    logLik
##   336.3607 358.6332 -161.1804
## 
## Random effects:
##  Formula: ~1 | ID
##         (Intercept) Residual
## StdDev:   0.3852014 0.489153
## 
## Fixed effects: meanLOTR ~ GROUP * WAVE + BASELINE 
##                  Value  Std.Error DF   t-value p-value
## (Intercept)  1.4004169 0.25479555 87  5.496238  0.0000
## GROUP1      -0.1224068 0.24951791 87 -0.490573  0.6250
## WAVE        -0.0492319 0.10700707 86 -0.460081  0.6466
## BASELINE     0.3721840 0.07424606 86  5.012845  0.0000
## GROUP1:WAVE -0.1337246 0.14884578 86 -0.898411  0.3715
##  Correlation: 
##             (Intr) GROUP1 WAVE   BASELI
## GROUP1      -0.512                     
## WAVE        -0.630  0.643              
## BASELINE    -0.710  0.008  0.000       
## GROUP1:WAVE  0.449 -0.895 -0.719  0.006
## 
## Standardized Within-Group Residuals:
##          Min           Q1          Med           Q3          Max 
## -1.942870672 -0.511480363 -0.008999573  0.426305211  2.736067156 
## 
## Number of Observations: 178
## Number of Groups: 89

Another random imputation

data15=MI$imputations[[15]]
library("Zelig")
zelig.fitdata15 <- zelig(meanLOTR ~ GROUP * WAVE + BASELINE, random = ~1 | ID, data = data15,  model = "ls", cite = FALSE)
summary(zelig.fitdata15)
## 
## Call:
## lm(formula = formula, weights = weights, model = F, data = data)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -2.02630 -0.38073 -0.04305  0.35135  1.99687 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  1.06953    0.26930   3.972 0.000105 ***
## GROUP1       0.05325    0.30498   0.175 0.861608    
## WAVE         0.14994    0.13862   1.082 0.280904    
## BASELINE     0.40216    0.06420   6.264 2.88e-09 ***
## GROUP1:WAVE -0.23544    0.19284  -1.221 0.223778    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.6427 on 173 degrees of freedom
## Multiple R-squared:  0.231,  Adjusted R-squared:  0.2132 
## F-statistic: 12.99 on 4 and 173 DF,  p-value: 2.845e-09

Describe the meanLOTR variable by the GROUP variable

describeBy(data15[,3:4], group = data15$GROUP)
## group: 0
##          vars  n mean   sd median trimmed  mad  min max range skew
## BASELINE    1 86 2.44 0.79    2.5    2.39 0.89 1.10 4.3  3.20 0.43
## meanLOTR    2 86 2.27 0.79    2.1    2.22 0.81 0.56 4.2  3.64 0.51
##          kurtosis   se
## BASELINE    -0.47 0.09
## meanLOTR    -0.49 0.09
## -------------------------------------------------------- 
## group: 1
##          vars  n mean   sd median trimmed  mad  min max range skew
## BASELINE    1 92 2.38 0.72   2.30    2.35 0.81 1.10   4  2.90 0.40
## meanLOTR    2 92 1.95 0.62   1.95    1.94 0.67 0.28   4  3.72 0.23
##          kurtosis   se
## BASELINE    -0.71 0.07
## meanLOTR     0.68 0.06

Create a plot that visualizes meanLOTR variable by the GROUP variable

library(ggplot2)
library(influence.ME)

Take a look at the residuals

residual <- lm(meanLOTR ~ BASELINE, data=data15)$residual

Plot the residuals to see that they are random

plot(density(residual))# A density plot

qqnorm(residual) # A quantile normal plot to checking normality
qqline(residual)

Checking the different between intervention and control groups residuals. This allows us to control for individual unsystematic differences.

data2$residual <- NA
sel1 <- which(!is.na(data15$meanLOTR)) 
sel2 <- which(!is.na(data15$BASELINE))
data15$residual[intersect(sel1,sel2)] <- residual
qplot(GROUP, meanLOTR, data=data15, geom="boxplot")

Plot of the difference between intervention and control groups.

qplot(GROUP, residual, data=data15, geom="boxplot")

Two way repeated measures ======================================================== Graphing the Two-Way Interaction. Both meanLOTR and the Residuals

# Load the nlme package
library(nlme)
with(data15, boxplot(meanLOTR ~ WAVE + GROUP))

with(data15, boxplot(residual ~ WAVE + GROUP))
Linear Mixed-Effects Model

Comparing Basline to Wave 2 and 3 by Group.

fullModeldata15 <- lme(meanLOTR ~ GROUP * WAVE + BASELINE, random = ~1 | ID, data = data15, method = "ML", na.action = "na.omit")

Cooks Distence

CookD(fullModeldata15)

Plot Cook’s distance:

plot(fullModeldata15, which="cook")
Check results on this random Imputation model
Results

Explanation of significance:

We asses the significance of our models by comparing them from the baseline model using the anova() function.
(Intercept): Where everything is 0
GROUP1: Is there a difference between group. If it is significant than there is a difference and the treatment had an effect.
WAVE: Asseses whether the effects gets bigger beteen time 2 and 3 (does not have to be significant)
BASELINE: Should not be significant. If it is then it shows that there is a difference between groups before the treatment.
GROUP1:WAVE: If this is significant then it means that the effect was either fleeting or it happened after the treatment i.e. between time 2 and 3.

summary(fullModeldata15)
## Linear mixed-effects model fit by maximum likelihood
##  Data: data15 
##        AIC      BIC    logLik
##   354.5097 376.7822 -170.2549
## 
## Random effects:
##  Formula: ~1 | ID
##         (Intercept)  Residual
## StdDev:   0.2506896 0.5819571
## 
## Fixed effects: meanLOTR ~ GROUP * WAVE + BASELINE 
##                  Value  Std.Error DF   t-value p-value
## (Intercept)  1.0711314 0.26489784 87  4.043564  0.0001
## GROUP1       0.0531588 0.28528678 87  0.186335  0.8526
## WAVE         0.1499395 0.12730888 86  1.177761  0.2421
## BASELINE     0.4014978 0.06883871 86  5.832442  0.0000
## GROUP1:WAVE -0.2354082 0.17711315 86 -1.329140  0.1873
##  Correlation: 
##             (Intr) GROUP1 WAVE   BASELI
## GROUP1      -0.576                     
## WAVE        -0.721  0.669              
## BASELINE    -0.633  0.032  0.000       
## GROUP1:WAVE  0.530 -0.932 -0.719 -0.019
## 
## Standardized Within-Group Residuals:
##         Min          Q1         Med          Q3         Max 
## -3.09433524 -0.52220563 -0.04242576  0.50662510  2.88563194 
## 
## Number of Observations: 178
## Number of Groups: 89

Another random imputation

data25=MI$imputations[[25]]

library("Zelig")
zelig.fitdata25 <- zelig(meanLOTR ~ GROUP * WAVE + BASELINE, random = ~1 | ID, data = data25,  model = "ls", cite = FALSE)
summary(zelig.fitdata25)
## 
## Call:
## lm(formula = formula, weights = weights, model = F, data = data)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.84632 -0.31827 -0.00072  0.35353  2.12655 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  0.84111    0.26379   3.189   0.0017 ** 
## GROUP1       0.20804    0.29851   0.697   0.4868    
## WAVE         0.22088    0.13569   1.628   0.1054    
## BASELINE     0.41851    0.06297   6.646 3.78e-10 ***
## GROUP1:WAVE -0.27710    0.18875  -1.468   0.1439    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.6292 on 173 degrees of freedom
## Multiple R-squared:  0.236,  Adjusted R-squared:  0.2183 
## F-statistic: 13.36 on 4 and 173 DF,  p-value: 1.657e-09

Describe the meanLOTR variable by the GROUP variable

describeBy(data25[,3:4], group = data25$GROUP)
## group: 0
##          vars  n mean   sd median trimmed  mad  min max range skew
## BASELINE    1 86 2.44 0.79   2.50    2.39 0.89 1.10 4.3  3.20 0.43
## meanLOTR    2 86 2.19 0.78   2.12    2.15 0.77 0.26 4.2  3.94 0.42
##          kurtosis   se
## BASELINE    -0.47 0.09
## meanLOTR    -0.02 0.08
## -------------------------------------------------------- 
## group: 1
##          vars  n mean   sd median trimmed  mad  min max range skew
## BASELINE    1 92 2.36 0.71   2.30    2.33 0.82 1.10   4  2.90 0.40
## meanLOTR    2 92 1.95 0.62   1.96    1.91 0.65 0.78   4  3.22 0.66
##          kurtosis   se
## BASELINE    -0.67 0.07
## meanLOTR     0.57 0.07

Create a plot that visualizes meanLOTR variable by the GROUP variable

library(ggplot2)
library(influence.ME)

Take a look at the residuals

residual <- lm(meanLOTR ~ BASELINE, data=data25)$residual

Plot the residuals to see that they are random

plot(density(residual))# A density plot

qqnorm(residual) # A quantile normal plot to checking normality
qqline(residual)

Checking the different between intervention and control groups residuals. This allows us to control for individual unsystematic differences.

data2$residual <- NA
sel1 <- which(!is.na(data25$meanLOTR)) 
sel2 <- which(!is.na(data25$BASELINE))
data25$residual[intersect(sel1,sel2)] <- residual
qplot(GROUP, meanLOTR, data=data25, geom="boxplot")

Plot of the difference between intervention and control groups.

qplot(GROUP, residual, data=data25, geom="boxplot")

Two way repeated measures ======================================================== Graphing the Two-Way Interaction. Both meanLOTR and the Residuals

# Load the nlme package
library(nlme)
with(data25, boxplot(meanLOTR ~ WAVE + GROUP))

with(data25, boxplot(residual ~ WAVE + GROUP))
Linear Mixed-Effects Model

Comparing Basline to Wave 2 and 3 by Group.

fullModeldata25 <- lme(meanLOTR ~ GROUP * WAVE + BASELINE, random = ~1 | ID, data = data25, method = "ML", na.action = "na.omit")

Cooks Distence

CookD(fullModeldata25)

Plot Cook’s distance:

plot(fullModeldata25, which="cook")
Check results on this random Imputation model
Results

Explanation of significance:

We asses the significance of our models by comparing them from the baseline model using the anova() function.
(Intercept): Where everything is 0
GROUP1: Is there a difference between group. If it is significant than there is a difference and the treatment had an effect.
WAVE: Asseses whether the effects gets bigger beteen time 2 and 3 (does not have to be significant)
BASELINE: Should not be significant. If it is then it shows that there is a difference between groups before the treatment.
GROUP1:WAVE: If this is significant then it means that the effect was either fleeting or it happened after the treatment i.e. between time 2 and 3.

summary(fullModeldata25)
## Linear mixed-effects model fit by maximum likelihood
##  Data: data25 
##        AIC      BIC    logLik
##   346.9598 369.2323 -166.4799
## 
## Random effects:
##  Formula: ~1 | ID
##         (Intercept)  Residual
## StdDev:   0.2439879 0.5702676
## 
## Fixed effects: meanLOTR ~ GROUP * WAVE + BASELINE 
##                  Value  Std.Error DF   t-value p-value
## (Intercept)  0.8408463 0.25972310 87  3.237472  0.0017
## GROUP1       0.2080476 0.27944525 87  0.744502  0.4586
## WAVE         0.2208819 0.12475167 86  1.770572  0.0802
## BASELINE     0.4186142 0.06757494 86  6.194814  0.0000
## GROUP1:WAVE -0.2771045 0.17353138 86 -1.596855  0.1140
##  Correlation: 
##             (Intr) GROUP1 WAVE   BASELI
## GROUP1      -0.573                     
## WAVE        -0.720  0.670              
## BASELINE    -0.634  0.027  0.000       
## GROUP1:WAVE  0.523 -0.932 -0.719 -0.008
## 
## Standardized Within-Group Residuals:
##           Min            Q1           Med            Q3           Max 
## -2.8948006423 -0.5172560659  0.0003487264  0.4847578463  3.1341393623 
## 
## Number of Observations: 178
## Number of Groups: 89

Check assumptions on model without any imputations

Describe the meanLOTR variable by the GROUP variable

describeBy(data2[,3:4], group = data2$GROUP)
## group: 0
##          vars  n mean   sd median trimmed  mad min max range skew kurtosis
## BASELINE    1 86 2.44 0.79    2.5    2.39 0.89 1.1 4.3   3.2 0.43    -0.47
## meanLOTR    2 59 2.19 0.80    2.1    2.11 0.74 1.1 4.2   3.1 0.86    -0.17
##            se
## BASELINE 0.09
## meanLOTR 0.10
## -------------------------------------------------------- 
## group: 1
##          vars  n mean   sd median trimmed  mad min max range skew kurtosis
## BASELINE    1 88 2.38 0.72   2.30    2.35 0.82 1.1   4   2.9 0.37    -0.69
## meanLOTR    2 54 1.94 0.57   1.85    1.88 0.52 1.1   4   2.9 1.22     1.97
##            se
## BASELINE 0.08
## meanLOTR 0.08

Create a plot that visualizes meanLOTR variable by the GROUP variable

library(ggplot2)
library(influence.ME)

Take a look at the residuals

residual <- lm(meanLOTR ~ BASELINE, data=data2)$residual

Plot the residuals to see that they are random

plot(density(residual))# A density plot

qqnorm(residual) # A quantile normal plot to checking normality
qqline(residual)

Checking the different between intervention and control groups residuals. This allows us to control for individual unsystematic differences.

data2$residual <- NA
sel1 <- which(!is.na(data2$meanLOTR)) 
sel2 <- which(!is.na(data2$BASELINE))
data2$residual[intersect(sel1,sel2)] <- residual
qplot(GROUP, meanLOTR, data=data2, geom="boxplot")
## Warning: Removed 65 rows containing non-finite values (stat_boxplot).

Plot of the difference between intervention and control groups.

qplot(GROUP, residual, data=data2, geom="boxplot")
## Warning: Removed 69 rows containing non-finite values (stat_boxplot).

Two way repeated measures ======================================================== Graphing the Two-Way Interaction. Both meanLOTR and the Residuals

# Load the nlme package
library(nlme)
with(data2, boxplot(meanLOTR ~ WAVE + GROUP))

with(data2, boxplot(residual ~ WAVE + GROUP))
Linear Mixed-Effects Model

Comparing Basline to Wave 2 and 3 by Group.

fullModel <- lme(meanLOTR ~ GROUP * WAVE + BASELINE, random = ~1 | ID, data = data2, method = "ML", na.action = "na.omit")

Cooks Distence

CookD(fullModel)

Plot Cook’s distance:

plot(fullModel, which="cook")
Results on Model with data that contains no imputations
Results

Explanation of significance:

We asses the significance of our models by comparing them from the baseline model using the anova() function.
(Intercept): Where everything is 0
GROUP1: Is there a difference between group. If it is significant than there is a difference and the treatment had an effect.
WAVE: Asseses whether the effects gets bigger beteen time 2 and 3 (does not have to be significant)
BASELINE: Should not be significant. If it is then it shows that there is a difference between groups before the treatment.
GROUP1:WAVE: If this is significant then it means that the effect was either fleeting or it happened after the treatment i.e. between time 2 and 3.

summary(fullModel)
## Linear mixed-effects model fit by maximum likelihood
##  Data: data2 
##        AIC      BIC    logLik
##   217.4728 236.3123 -101.7364
## 
## Random effects:
##  Formula: ~1 | ID
##         (Intercept)  Residual
## StdDev:   0.3868753 0.5005022
## 
## Fixed effects: meanLOTR ~ GROUP * WAVE + BASELINE 
##                  Value Std.Error DF   t-value p-value
## (Intercept)  1.2002102 0.3122824 66  3.843348  0.0003
## GROUP1       0.0476590 0.3244938 66  0.146872  0.8837
## WAVE         0.0636678 0.1435847 38  0.443416  0.6600
## BASELINE     0.3966105 0.0910289 66  4.356972  0.0000
## GROUP1:WAVE -0.1889549 0.2102791 38 -0.898591  0.3745
##  Correlation: 
##             (Intr) GROUP1 WAVE   BASELI
## GROUP1      -0.513                     
## WAVE        -0.668  0.614              
## BASELINE    -0.711  0.052  0.047       
## GROUP1:WAVE  0.457 -0.903 -0.683 -0.032
## 
## Standardized Within-Group Residuals:
##        Min         Q1        Med         Q3        Max 
## -1.9149812 -0.4879502 -0.0951565  0.3848770  2.6604162 
## 
## Number of Observations: 109
## Number of Groups: 69
Table with P-value

|             |       Value|  Std.Error|  DF|     t-value|    p-value|
|:------------|-----------:|----------:|---:|-----------:|----------:|
|(Intercept)  |   1.2002102|  0.3122824|  66|   3.8433481|  0.0002754|
|GROUP1       |   0.0476590|  0.3244938|  66|   0.1468719|  0.8836808|
|WAVE         |   0.0636678|  0.1435847|  38|   0.4434164|  0.6599790|
|BASELINE     |   0.3966105|  0.0910289|  66|   4.3569722|  0.0000471|
|GROUP1:WAVE  |  -0.1889549|  0.2102791|  38|  -0.8985911|  0.3745293|

Table with confidence intervals

est. lower upper
(Intercept) 1.2002102 0.5911863 1.8092340
GROUP1 0.0476590 -0.5851798 0.6804978
WAVE 0.0636678 -0.2202592 0.3475949
BASELINE 0.3966105 0.2190828 0.5741382
GROUP1:WAVE -0.1889549 -0.6047645 0.2268547