Loading the dataset

data.test4 <- read.csv("/Volumes/TOSHIBA EXT/Dropbox/ADULT STUDY/adult_study011615.csv")
# Load the psych package
library(psych)
items <- c("LET1", "LET2", "LET3", "LET4", "LET5", "LET6")
scaleKey <- c(-1,1,-1,1,-1,1)
data.test4$meanLET  <- scoreItems(scaleKey, items=data.test4[,items], delete=FALSE)$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 meanLET and ID Group and wave from data.test4 and create a new #dataset with only those variables.
data <- data.test4[,c("ID", "GROUP", "wave", "meanLET")]
#Use dcast to cnage from long-format data to wide format data
data <- dcast(data, ID + GROUP ~ wave, mean, value.var = "meanLET")
# 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 meanLET and wave 2 and 3 of meanLET. 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", "meanLET", "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
## 
## -- 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  9
## 
## -- Imputation 6 --
## 
##   1  2  3  4  5  6  7  8  9 10
## 
## -- Imputation 7 --
## 
##   1  2  3  4  5  6  7
## 
## -- Imputation 8 --
## 
##   1  2  3  4  5  6  7
## 
## -- Imputation 9 --
## 
##   1  2  3  4  5  6  7
## 
## -- Imputation 10 --
## 
##   1  2  3  4  5  6
## 
## -- Imputation 11 --
## 
##   1  2  3  4  5  6  7
## 
## -- Imputation 12 --
## 
##   1  2  3  4  5  6
## 
## -- 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  8  9
## 
## -- Imputation 16 --
## 
##   1  2  3  4  5  6  7  8
## 
## -- Imputation 17 --
## 
##   1  2  3  4  5  6  7  8
## 
## -- Imputation 18 --
## 
##   1  2  3  4  5  6  7  8  9
## 
## -- Imputation 19 --
## 
##   1  2  3  4  5  6  7
## 
## -- Imputation 20 --
## 
##   1  2  3  4  5  6  7  8  9
## 
## -- Imputation 21 --
## 
##   1  2  3  4  5  6  7  8
## 
## -- Imputation 22 --
## 
##   1  2  3  4  5  6  7  8  9
## 
## -- Imputation 23 --
## 
##   1  2  3  4  5  6  7  8  9 10
## 
## -- Imputation 24 --
## 
##   1  2  3  4  5  6  7  8  9 10
## 
## -- Imputation 25 --
## 
##   1  2  3  4
## 
## -- Imputation 26 --
## 
##   1  2  3  4  5  6  7
## 
## -- Imputation 27 --
## 
##   1  2  3  4  5  6  7  8
## 
## -- Imputation 28 --
## 
##   1  2  3  4  5  6  7  8  9
## 
## -- Imputation 29 --
## 
##   1  2  3  4  5  6  7
## 
## -- Imputation 30 --
## 
##   1  2  3  4  5  6  7  8  9
## 
## -- Imputation 31 --
## 
##   1  2  3  4  5  6  7  8  9
## 
## -- Imputation 32 --
## 
##   1  2  3  4  5  6  7  8
## 
## -- Imputation 33 --
## 
##   1  2  3  4  5  6  7
## 
## -- Imputation 34 --
## 
##   1  2  3  4  5  6  7
## 
## -- Imputation 35 --
## 
##   1  2  3  4  5  6  7  8
## 
## -- Imputation 36 --
## 
##   1  2  3  4  5  6  7  8  9 10 11
## 
## -- Imputation 37 --
## 
##   1  2  3  4  5  6  7  8  9
## 
## -- Imputation 38 --
## 
##   1  2  3  4  5  6  7  8
## 
## -- Imputation 39 --
## 
##   1  2  3  4  5  6  7  8  9 10 11 12
## 
## -- Imputation 40 --
## 
##   1  2  3  4  5  6  7  8  9
## 
## -- Imputation 41 --
## 
##   1  2  3  4  5  6  7  8  9 10 11 12
## 
## -- Imputation 42 --
## 
##   1  2  3  4  5  6  7  8  9
## 
## -- Imputation 43 --
## 
##   1  2  3  4  5  6  7
## 
## -- Imputation 44 --
## 
##   1  2  3  4  5  6  7  8
## 
## -- Imputation 45 --
## 
##   1  2  3  4  5  6  7  8  9 10
## 
## -- Imputation 46 --
## 
##   1  2  3  4  5  6
## 
## -- Imputation 47 --
## 
##   1  2  3  4  5  6
## 
## -- Imputation 48 --
## 
##   1  2  3  4  5  6  7  8  9 10
## 
## -- Imputation 49 --
## 
##   1  2  3  4  5  6  7  8  9 10 11
## 
## -- Imputation 50 --
## 
##   1  2  3  4  5  6  7  8  9 10

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(meanLET ~ 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.19461911 0.31930459  0.5668209 1.8224173     36 %
## GROUP1       0.01620297 0.26479893 -0.5037547 0.5361607     28 %
## WAVE        -0.08811118 0.12577137 -0.3354006 0.1591783     36 %
## BASELINE     0.75760294 0.07166653  0.6163538 0.8988521     48 %
## GROUP1:WAVE  0.21798581 0.16940359 -0.1148553 0.5508269     32 %

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(meanLET ~ 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.19024841 0.31956994  3.72453184 2.244588e-04
## GROUP1       0.01648232 0.27992142  0.05888195 9.530605e-01
## WAVE        -0.08811118 0.13365924 -0.65922255 5.100638e-01
## BASELINE     0.75874775 0.06960917 10.90011207 5.750087e-22
## GROUP1:WAVE  0.21796412 0.18070908  1.20616032 2.282019e-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(meanLET ~ 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.76548 -0.29103  0.00072  0.27887  1.40471 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  1.214531   0.245589   4.945 1.79e-06 ***
## GROUP1       0.009413   0.233995   0.040    0.968    
## WAVE        -0.088375   0.106156  -0.833    0.406    
## BASELINE     0.754303   0.046959  16.063  < 2e-16 ***
## GROUP1:WAVE  0.259281   0.147699   1.755    0.081 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.4922 on 173 degrees of freedom
## Multiple R-squared:  0.6124, Adjusted R-squared:  0.6034 
## F-statistic: 68.32 on 4 and 173 DF,  p-value: < 2.2e-16

Describe the meanLET 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 3.82 0.70   3.83    3.83 0.74 2.50 5.00  2.50 -0.23
## meanLET     2 86 3.96 0.78   4.03    4.01 0.94 1.83 5.08  3.25 -0.53
##          kurtosis   se
## BASELINE    -0.90 0.08
## meanLET     -0.46 0.08
## -------------------------------------------------------- 
## group: 1
##          vars  n mean   sd median trimmed  mad  min  max range  skew
## BASELINE    1 92 3.59 0.86   3.58    3.60 0.86 2.00 5.11  3.11 -0.07
## meanLET     2 92 4.19 0.77   4.17    4.23 0.85 2.17 5.76  3.59 -0.37
##          kurtosis   se
## BASELINE    -0.99 0.09
## meanLET     -0.41 0.08

Create a plot that visualizes meanLET 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(meanLET ~ 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$meanLET)) 
sel2 <- which(!is.na(data1$BASELINE))
data1$residual[intersect(sel1,sel2)] <- residual
qplot(GROUP, meanLET, 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 meanLET and the Residuals

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

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

Comparing Basline to Wave 2 and 3 by Group.

fullModeldata1 <- lme(meanLET ~ 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
##   256.6736 278.9461 -121.3368
## 
## Random effects:
##  Formula: ~1 | ID
##         (Intercept)  Residual
## StdDev:    0.235773 0.4242031
## 
## Fixed effects: meanLET ~ GROUP * WAVE + BASELINE 
##                  Value  Std.Error DF   t-value p-value
## (Intercept)  1.2519196 0.24897313 87  5.028332  0.0000
## GROUP1       0.0061428 0.21101349 87  0.029111  0.9768
## WAVE        -0.0883748 0.09279862 86 -0.952329  0.3436
## BASELINE     0.7445098 0.05181198 86 14.369451  0.0000
## GROUP1:WAVE  0.2600041 0.12913630 86  2.013408  0.0472
##  Correlation: 
##             (Intr) GROUP1 WAVE   BASELI
## GROUP1      -0.500                     
## WAVE        -0.559  0.660              
## BASELINE    -0.795  0.082  0.000       
## GROUP1:WAVE  0.425 -0.920 -0.719 -0.030
## 
## Standardized Within-Group Residuals:
##          Min           Q1          Med           Q3          Max 
## -3.423233630 -0.511156998  0.005912427  0.517379535  2.119213021 
## 
## Number of Observations: 178
## Number of Groups: 89

Another random imputation

data10=MI$imputations[[10]]
library("Zelig")
zelig.fitdata10 <- zelig(meanLET ~ 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.76515 -0.26932  0.05965  0.28325  1.21746 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  1.24347    0.25614   4.855 2.68e-06 ***
## GROUP1       0.12442    0.24197   0.514    0.608    
## WAVE        -0.10665    0.10992  -0.970    0.333    
## BASELINE     0.75592    0.04928  15.338  < 2e-16 ***
## GROUP1:WAVE  0.11032    0.15289   0.722    0.472    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.5097 on 173 degrees of freedom
## Multiple R-squared:  0.5802, Adjusted R-squared:  0.5705 
## F-statistic: 59.77 on 4 and 173 DF,  p-value: < 2.2e-16

Describe the meanLET 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 3.82 0.70   3.83    3.83 0.74 2.50 5.00  2.50 -0.23
## meanLET     2 86 3.97 0.77   4.05    4.02 0.91 1.83 5.37  3.54 -0.55
##          kurtosis   se
## BASELINE    -0.90 0.08
## meanLET     -0.27 0.08
## -------------------------------------------------------- 
## group: 1
##          vars  n mean   sd median trimmed  mad  min  max range  skew
## BASELINE    1 92 3.61 0.84   3.67    3.62 0.74 2.00 5.00  3.00 -0.15
## meanLET     2 92 4.10 0.78   4.17    4.16 0.82 1.91 5.99  4.08 -0.53
##          kurtosis   se
## BASELINE    -0.93 0.09
## meanLET      0.00 0.08

Create a plot that visualizes meanLET variable by the GROUP variable

library(ggplot2)
library(influence.ME)

Take a look at the residuals

residual <- lm(meanLET ~ 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$meanLET)) 
sel2 <- which(!is.na(data10$BASELINE))
data10$residual[intersect(sel1,sel2)] <- residual
qplot(GROUP, meanLET, 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 meanLET and the Residuals

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

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

Comparing Basline to Wave 2 and 3 by Group.

fullModeldata10 <- lme(meanLET ~ 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
##   272.142 294.4145 -129.071
## 
## Random effects:
##  Formula: ~1 | ID
##         (Intercept)  Residual
## StdDev:   0.1935451 0.4636793
## 
## Fixed effects: meanLET ~ GROUP * WAVE + BASELINE 
##                  Value  Std.Error DF   t-value p-value
## (Intercept)  1.2444061 0.25935628 87  4.798057  0.0000
## GROUP1       0.1243663 0.22721572 87  0.547349  0.5855
## WAVE        -0.1066509 0.10143443 86 -1.051427  0.2960
## BASELINE     0.7556760 0.05280775 86 14.309946  0.0000
## GROUP1:WAVE  0.1103208 0.14109171 86  0.781908  0.4364
##  Correlation: 
##             (Intr) GROUP1 WAVE   BASELI
## GROUP1      -0.490                     
## WAVE        -0.587  0.670              
## BASELINE    -0.777  0.049  0.000       
## GROUP1:WAVE  0.422 -0.931 -0.719  0.000
## 
## Standardized Within-Group Residuals:
##         Min          Q1         Med          Q3         Max 
## -3.35041203 -0.51409409  0.07754312  0.52205476  2.27745853 
## 
## Number of Observations: 178
## Number of Groups: 89

Another random imputation

data15=MI$imputations[[15]]
library("Zelig")
zelig.fitdata15 <- zelig(meanLET ~ 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 
## -1.76579 -0.29131  0.01747  0.26648  1.28830 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  0.87291    0.25084   3.480 0.000635 ***
## GROUP1       0.09590    0.23723   0.404 0.686520    
## WAVE        -0.12646    0.10774  -1.174 0.242095    
## BASELINE     0.84722    0.04823  17.568  < 2e-16 ***
## GROUP1:WAVE  0.19442    0.14987   1.297 0.196259    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.4996 on 173 degrees of freedom
## Multiple R-squared:  0.6478, Adjusted R-squared:  0.6396 
## F-statistic: 79.53 on 4 and 173 DF,  p-value: < 2.2e-16

Describe the meanLET 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 3.82 0.70   3.83    3.83 0.74 2.50 5.00  2.50 -0.23
## meanLET     2 86 3.92 0.83   4.00    3.97 0.99 1.83 5.29  3.46 -0.44
##          kurtosis   se
## BASELINE    -0.90 0.08
## meanLET     -0.75 0.09
## -------------------------------------------------------- 
## group: 1
##          vars  n mean   sd median trimmed  mad  min  max range  skew
## BASELINE    1 92 3.60 0.85   3.67    3.62 0.74 2.00 5.00  3.00 -0.13
## meanLET     2 92 4.12 0.83   4.17    4.18 0.85 1.59 5.72  4.12 -0.61
##          kurtosis   se
## BASELINE    -0.94 0.09
## meanLET      0.03 0.09

Create a plot that visualizes meanLET variable by the GROUP variable

library(ggplot2)
library(influence.ME)

Take a look at the residuals

residual <- lm(meanLET ~ 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$meanLET)) 
sel2 <- which(!is.na(data15$BASELINE))
data15$residual[intersect(sel1,sel2)] <- residual
qplot(GROUP, meanLET, 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 meanLET and the Residuals

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

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

Comparing Basline to Wave 2 and 3 by Group.

fullModeldata15 <- lme(meanLET ~ 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
##   265.7446 288.0171 -125.8723
## 
## Random effects:
##  Formula: ~1 | ID
##         (Intercept)  Residual
## StdDev:   0.1696558 0.4623641
## 
## Fixed effects: meanLET ~ GROUP * WAVE + BASELINE 
##                  Value  Std.Error DF   t-value p-value
## (Intercept)  0.8728462 0.25332195 87  3.445601  0.0009
## GROUP1       0.0959081 0.22575507 87  0.424833  0.6720
## WAVE        -0.1264644 0.10114672 86 -1.250307  0.2146
## BASELINE     0.8472334 0.05099659 86 16.613530  0.0000
## GROUP1:WAVE  0.1944190 0.14069376 86  1.381860  0.1706
##  Correlation: 
##             (Intr) GROUP1 WAVE   BASELI
## GROUP1      -0.501                     
## WAVE        -0.599  0.672              
## BASELINE    -0.769  0.054  0.000       
## GROUP1:WAVE  0.435 -0.935 -0.719 -0.006
## 
## Standardized Within-Group Residuals:
##          Min           Q1          Med           Q3          Max 
## -3.438678698 -0.538050860 -0.001059709  0.525879724  2.735928431 
## 
## Number of Observations: 178
## Number of Groups: 89

Another random imputation

data25=MI$imputations[[25]]

library("Zelig")
zelig.fitdata25 <- zelig(meanLET ~ 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.7648 -0.2759  0.0393  0.2795  1.4578 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  1.51947    0.26606   5.711  4.8e-08 ***
## GROUP1      -0.08916    0.25230  -0.353    0.724    
## WAVE        -0.11009    0.11454  -0.961    0.338    
## BASELINE     0.69603    0.05106  13.633  < 2e-16 ***
## GROUP1:WAVE  0.22488    0.15933   1.411    0.160    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.5311 on 173 degrees of freedom
## Multiple R-squared:  0.5236, Adjusted R-squared:  0.5126 
## F-statistic: 47.53 on 4 and 173 DF,  p-value: < 2.2e-16

Describe the meanLET 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 3.82 0.70   3.83    3.83 0.74 2.50 5.00  2.50 -0.23
## meanLET     2 86 4.01 0.77   4.17    4.07 0.74 1.83 5.42  3.59 -0.59
##          kurtosis   se
## BASELINE    -0.90 0.08
## meanLET     -0.21 0.08
## -------------------------------------------------------- 
## group: 1
##          vars  n mean   sd median trimmed  mad  min  max range  skew
## BASELINE    1 92 3.61 0.85   3.67    3.62 0.86 2.00 5.00  3.00 -0.13
## meanLET     2 92 4.11 0.75   4.17    4.18 0.82 2.17 5.47  3.31 -0.60
##          kurtosis   se
## BASELINE    -0.98 0.09
## meanLET     -0.36 0.08

Create a plot that visualizes meanLET variable by the GROUP variable

library(ggplot2)
library(influence.ME)

Take a look at the residuals

residual <- lm(meanLET ~ 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$meanLET)) 
sel2 <- which(!is.na(data25$BASELINE))
data25$residual[intersect(sel1,sel2)] <- residual
qplot(GROUP, meanLET, 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 meanLET and the Residuals

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

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

Comparing Basline to Wave 2 and 3 by Group.

fullModeldata25 <- lme(meanLET ~ 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
##   288.1751 310.4476 -137.0876
## 
## Random effects:
##  Formula: ~1 | ID
##         (Intercept)  Residual
## StdDev:   0.1501555 0.5015752
## 
## Fixed effects: meanLET ~ GROUP * WAVE + BASELINE 
##                  Value  Std.Error DF   t-value p-value
## (Intercept)  1.5222982 0.26773950 87  5.685744  0.0000
## GROUP1      -0.0893649 0.24391532 87 -0.366377  0.7150
## WAVE        -0.1100928 0.10972454 86 -1.003356  0.3185
## BASELINE     0.6952895 0.05306648 86 13.102237  0.0000
## GROUP1:WAVE  0.2249167 0.15264099 86  1.473501  0.1443
##  Correlation: 
##             (Intr) GROUP1 WAVE   BASELI
## GROUP1      -0.515                     
## WAVE        -0.615  0.675              
## BASELINE    -0.757  0.060  0.000       
## GROUP1:WAVE  0.454 -0.939 -0.719 -0.015
## 
## Standardized Within-Group Residuals:
##         Min          Q1         Med          Q3         Max 
## -3.26913656 -0.50623027  0.06696432  0.46405913  2.86479759 
## 
## Number of Observations: 178
## Number of Groups: 89

Check assumptions on model without any imputations

Describe the meanLET 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
## BASELINE    1 86 3.82 0.70   3.83    3.83 0.74 2.50   5  2.50 -0.23
## meanLET     2 59 4.05 0.81   4.17    4.14 0.74 1.83   5  3.17 -0.86
##          kurtosis   se
## BASELINE    -0.90 0.08
## meanLET     -0.03 0.11
## -------------------------------------------------------- 
## group: 1
##          vars  n mean   sd median trimmed  mad  min max range  skew
## BASELINE    1 88 3.59 0.86   3.58    3.60 0.86 2.00   5  3.00 -0.09
## meanLET     2 54 4.20 0.72   4.25    4.29 0.74 2.17   5  2.83 -0.98
##          kurtosis   se
## BASELINE    -0.98 0.09
## meanLET      0.49 0.10

Create a plot that visualizes meanLET variable by the GROUP variable

library(ggplot2)
library(influence.ME)

Take a look at the residuals

residual <- lm(meanLET ~ 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$meanLET)) 
sel2 <- which(!is.na(data2$BASELINE))
data2$residual[intersect(sel1,sel2)] <- residual
qplot(GROUP, meanLET, 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 meanLET and the Residuals

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

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

Comparing Basline to Wave 2 and 3 by Group.

fullModel <- lme(meanLET ~ 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
##   168.0512 186.8907 -77.02562
## 
## Random effects:
##  Formula: ~1 | ID
##         (Intercept)  Residual
## StdDev:   0.3759615 0.3593626
## 
## Fixed effects: meanLET ~ GROUP * WAVE + BASELINE 
##                  Value Std.Error DF   t-value p-value
## (Intercept)  1.2053278 0.3373420 66  3.573014  0.0007
## GROUP1      -0.0243520 0.2443397 66 -0.099665  0.9209
## WAVE        -0.0945574 0.1054288 38 -0.896884  0.3754
## BASELINE     0.7535852 0.0772543 66  9.754605  0.0000
## GROUP1:WAVE  0.2524914 0.1540213 38  1.639328  0.1094
##  Correlation: 
##             (Intr) GROUP1 WAVE   BASELI
## GROUP1      -0.384                     
## WAVE        -0.387  0.590              
## BASELINE    -0.871  0.060 -0.048       
## GROUP1:WAVE  0.273 -0.871 -0.684  0.024
## 
## Standardized Within-Group Residuals:
##         Min          Q1         Med          Q3         Max 
## -3.33870466 -0.39782905  0.03127719  0.43784183  1.96357889 
## 
## Number of Observations: 109
## Number of Groups: 69
Table with P-value

|             |       Value|  Std.Error|  DF|     t-value|    p-value|
|:------------|-----------:|----------:|---:|-----------:|----------:|
|(Intercept)  |   1.2053278|  0.3373420|  66|   3.5730142|  0.0006652|
|GROUP1       |  -0.0243520|  0.2443397|  66|  -0.0996646|  0.9209127|
|WAVE         |  -0.0945574|  0.1054288|  38|  -0.8968841|  0.3754279|
|BASELINE     |   0.7535852|  0.0772543|  66|   9.7546051|  0.0000000|
|GROUP1:WAVE  |   0.2524914|  0.1540213|  38|   1.6393282|  0.1093991|

Table with confidence intervals

est. lower upper
(Intercept) 1.2053278 0.5474319 1.8632237
GROUP1 -0.0243520 -0.5008716 0.4521676
WAVE -0.0945574 -0.3030341 0.1139194
BASELINE 0.7535852 0.6029212 0.9042492
GROUP1:WAVE 0.2524914 -0.0520730 0.5570558