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
data2 <- data2[-c(which(data2$GROUP ==2)),]

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  8
## 
## -- Imputation 2 --
## 
##   1  2  3  4  5  6
## 
## -- Imputation 3 --
## 
##   1  2  3  4  5  6  7  8
## 
## -- Imputation 4 --
## 
##   1  2  3  4  5  6  7
## 
## -- Imputation 5 --
## 
##   1  2  3  4  5  6  7
## 
## -- Imputation 6 --
## 
##   1  2  3  4  5  6  7  8  9
## 
## -- Imputation 7 --
## 
##   1  2  3  4  5  6  7
## 
## -- Imputation 8 --
## 
##   1  2  3  4  5  6  7  8  9
## 
## -- Imputation 9 --
## 
##   1  2  3  4  5  6  7
## 
## -- Imputation 10 --
## 
##   1  2  3  4  5
## 
## -- Imputation 11 --
## 
##   1  2  3  4  5  6
## 
## -- Imputation 12 --
## 
##   1  2  3  4  5  6
## 
## -- Imputation 13 --
## 
##   1  2  3  4  5  6  7  8
## 
## -- Imputation 14 --
## 
##   1  2  3  4  5  6  7  8  9
## 
## -- Imputation 15 --
## 
##   1  2  3  4  5  6
## 
## -- Imputation 16 --
## 
##   1  2  3  4  5  6  7  8  9
## 
## -- Imputation 17 --
## 
##   1  2  3  4  5  6  7  8  9 10 11 12 13
## 
## -- Imputation 18 --
## 
##   1  2  3  4  5  6  7  8
## 
## -- Imputation 19 --
## 
##   1  2  3  4  5  6  7  8  9
## 
## -- Imputation 20 --
## 
##   1  2  3  4  5  6  7  8  9
## 
## -- Imputation 21 --
## 
##   1  2  3  4  5  6  7  8  9
## 
## -- Imputation 22 --
## 
##   1  2  3  4  5  6  7
## 
## -- Imputation 23 --
## 
##   1  2  3  4  5  6  7  8  9
## 
## -- Imputation 24 --
## 
##   1  2  3  4  5  6  7
## 
## -- Imputation 25 --
## 
##   1  2  3  4  5  6  7  8
## 
## -- 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
## 
## -- Imputation 29 --
## 
##   1  2  3  4  5  6  7  8
## 
## -- Imputation 30 --
## 
##   1  2  3  4  5  6  7  8  9
## 
## -- Imputation 31 --
## 
##   1  2  3  4  5  6  7
## 
## -- Imputation 32 --
## 
##   1  2  3  4  5
## 
## -- Imputation 33 --
## 
##   1  2  3  4  5  6  7
## 
## -- Imputation 34 --
## 
##   1  2  3  4  5  6  7  8
## 
## -- Imputation 35 --
## 
##   1  2  3  4  5  6  7
## 
## -- Imputation 36 --
## 
##   1  2  3  4  5  6
## 
## -- Imputation 37 --
## 
##   1  2  3  4  5  6  7
## 
## -- 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
## 
## -- Imputation 41 --
## 
##   1  2  3  4  5  6  7  8
## 
## -- Imputation 42 --
## 
##   1  2  3  4  5  6  7
## 
## -- Imputation 43 --
## 
##   1  2  3  4  5  6  7  8  9 10 11 12 13
## 
## -- Imputation 44 --
## 
##   1  2  3  4  5  6  7
## 
## -- Imputation 45 --
## 
##   1  2  3  4  5  6  7  8  9
## 
## -- Imputation 46 --
## 
##   1  2  3  4  5  6  7  8
## 
## -- Imputation 47 --
## 
##   1  2  3  4  5  6  7  8  9 10
## 
## -- Imputation 48 --
## 
##   1  2  3  4  5  6
## 
## -- Imputation 49 --
## 
##   1  2  3  4  5  6  7
## 
## -- Imputation 50 --
## 
##   1  2  3  4  5  6  7  8

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.08052002 0.31342997  0.4646095 1.6964306     33 %
## GROUP1      -0.01167188 0.27549694 -0.5529791 0.5296354     32 %
## WAVE        -0.10147979 0.13313300 -0.3636762 0.1607166     45 %
## BASELINE     0.79292524 0.07005706  0.6549879 0.9308626     43 %
## GROUP1:WAVE  0.21250897 0.18740936 -0.1564856 0.5815036     43 %

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.07734069 0.31305021  3.44143099 6.294522e-04
## GROUP1      -0.01149875 0.29009342 -0.03963809 9.683949e-01
## WAVE        -0.10147979 0.14028034 -0.72340710 4.699747e-01
## BASELINE     0.79375800 0.06797235 11.67766067 2.917485e-25
## GROUP1:WAVE  0.21248612 0.19768104  1.07489376 2.832094e-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.73460 -0.26834  0.07787  0.29745  1.28131 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  0.98154    0.25953   3.782 0.000217 ***
## GROUP1      -0.09882    0.24877  -0.397 0.691718    
## WAVE        -0.13926    0.11048  -1.261 0.209255    
## BASELINE     0.82183    0.05028  16.346  < 2e-16 ***
## GROUP1:WAVE  0.30409    0.15719   1.935 0.054750 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.5123 on 165 degrees of freedom
## Multiple R-squared:  0.6296, Adjusted R-squared:  0.6206 
## F-statistic: 70.12 on 4 and 165 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.5 -0.23
## meanLET     2 86 3.91 0.83   4.10    3.97 0.87 1.83 5.04   3.2 -0.53
##          kurtosis   se
## BASELINE    -0.90 0.08
## meanLET     -0.71 0.09
## -------------------------------------------------------- 
## group: 1
##          vars  n mean   sd median trimmed  mad  min  max range  skew
## BASELINE    1 84 3.66 0.86   3.67    3.66 0.77 2.00 5.75  3.75  0.00
## meanLET     2 84 4.14 0.83   4.17    4.20 0.78 2.12 5.69  3.57 -0.58
##          kurtosis   se
## BASELINE    -0.79 0.09
## meanLET     -0.26 0.09

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
##   263.6734 285.624 -124.8367
## 
## Random effects:
##  Formula: ~1 | ID
##         (Intercept)  Residual
## StdDev:   0.1193405 0.4903714
## 
## Fixed effects: meanLET ~ GROUP * WAVE + BASELINE 
##                  Value  Std.Error DF   t-value p-value
## (Intercept)  0.9927666 0.26057001 83  3.809980  0.0003
## GROUP1      -0.0994962 0.24317517 83 -0.409154  0.6835
## WAVE        -0.1392626 0.10734661 82 -1.297317  0.1982
## BASELINE     0.8188922 0.05155893 82 15.882647  0.0000
## GROUP1:WAVE  0.3042340 0.15273210 82  1.991946  0.0497
##  Correlation: 
##             (Intr) GROUP1 WAVE   BASELI
## GROUP1      -0.497                     
## WAVE        -0.618  0.662              
## BASELINE    -0.755  0.049  0.000       
## GROUP1:WAVE  0.447 -0.943 -0.703 -0.016
## 
## Standardized Within-Group Residuals:
##        Min         Q1        Med         Q3        Max 
## -3.3599567 -0.5582208  0.1231809  0.5494301  2.5464484 
## 
## Number of Observations: 170
## Number of Groups: 85

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.79268 -0.29738  0.06134  0.31003  1.19753 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  1.28749    0.26409   4.875 2.54e-06 ***
## GROUP1      -0.05807    0.25032  -0.232    0.817    
## WAVE        -0.09640    0.11122  -0.867    0.387    
## BASELINE     0.74770    0.05161  14.489  < 2e-16 ***
## GROUP1:WAVE  0.20861    0.15823   1.318    0.189    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.5157 on 165 degrees of freedom
## Multiple R-squared:  0.5652, Adjusted R-squared:  0.5547 
## F-statistic: 53.62 on 4 and 165 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 4.00 0.79   4.17    4.05 0.87 1.83 5.26  3.43 -0.58
##          kurtosis   se
## BASELINE    -0.90 0.08
## meanLET     -0.44 0.09
## -------------------------------------------------------- 
## group: 1
##          vars  n mean   sd median trimmed  mad  min  max range  skew
## BASELINE    1 84 3.65 0.83   3.67    3.67 0.74 2.00 5.00  3.00 -0.16
## meanLET     2 84 4.13 0.76   4.22    4.20 0.66 2.17 5.33  3.16 -0.78
##          kurtosis   se
## BASELINE    -0.94 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
##   260.6133 282.5639 -123.3066
## 
## Random effects:
##  Formula: ~1 | ID
##         (Intercept)  Residual
## StdDev:   0.2553815 0.4392334
## 
## Fixed effects: meanLET ~ GROUP * WAVE + BASELINE 
##                  Value  Std.Error DF   t-value p-value
## (Intercept)  1.2907011 0.27066287 83  4.768667  0.0000
## GROUP1      -0.0582032 0.22366738 83 -0.260222  0.7953
## WAVE        -0.0963982 0.09615206 82 -1.002560  0.3190
## BASELINE     0.7468597 0.05773299 82 12.936447  0.0000
## GROUP1:WAVE  0.2086098 0.13678655 82  1.525076  0.1311
##  Correlation: 
##             (Intr) GROUP1 WAVE   BASELI
## GROUP1      -0.442                     
## WAVE        -0.533  0.645              
## BASELINE    -0.814  0.042  0.000       
## GROUP1:WAVE  0.374 -0.917 -0.703  0.000
## 
## Standardized Within-Group Residuals:
##         Min          Q1         Med          Q3         Max 
## -3.30836296 -0.46498123  0.06892531  0.50577161  2.08794787 
## 
## Number of Observations: 170
## Number of Groups: 85

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.80508 -0.26728  0.07447  0.33266  1.24076 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 1.281253   0.266851   4.801 3.51e-06 ***
## GROUP1      0.233793   0.253221   0.923    0.357    
## WAVE        0.006305   0.112432   0.056    0.955    
## BASELINE    0.706196   0.052128  13.547  < 2e-16 ***
## GROUP1:WAVE 0.022925   0.159961   0.143    0.886    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.5213 on 165 degrees of freedom
## Multiple R-squared:  0.5309, Adjusted R-squared:  0.5196 
## F-statistic: 46.69 on 4 and 165 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.7   3.83    3.83 0.74 2.50 5.00  2.50 -0.23
## meanLET     2 86 3.99 0.8   4.15    4.05 0.77 1.83 5.49  3.66 -0.61
##          kurtosis   se
## BASELINE    -0.90 0.08
## meanLET     -0.25 0.09
## -------------------------------------------------------- 
## group: 1
##          vars  n mean   sd median trimmed  mad  min  max range  skew
## BASELINE    1 84 3.63 0.83   3.67    3.65 0.74 2.00 5.00  3.00 -0.11
## meanLET     2 84 4.12 0.69   4.17    4.19 0.74 2.17 5.21  3.04 -0.64
##          kurtosis   se
## BASELINE    -0.96 0.09
## meanLET     -0.03 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=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
##   268.938 290.8886 -127.469
## 
## Random effects:
##  Formula: ~1 | ID
##         (Intercept)  Residual
## StdDev:   0.1671064 0.4856585
## 
## Fixed effects: meanLET ~ GROUP * WAVE + BASELINE 
##                 Value  Std.Error DF   t-value p-value
## (Intercept) 1.2842849 0.26951171 83  4.765228  0.0000
## GROUP1      0.2335978 0.24232561 83  0.963983  0.3379
## WAVE        0.0063047 0.10631493 82  0.059302  0.9529
## BASELINE    0.7054020 0.05476151 82 12.881345  0.0000
## GROUP1:WAVE 0.0229572 0.15126072 82  0.151773  0.8797
##  Correlation: 
##             (Intr) GROUP1 WAVE   BASELI
## GROUP1      -0.486                     
## WAVE        -0.592  0.658              
## BASELINE    -0.776  0.055  0.000       
## GROUP1:WAVE  0.427 -0.937 -0.703 -0.015
## 
## Standardized Within-Group Residuals:
##        Min         Q1        Med         Q3        Max 
## -3.3997342 -0.4956551  0.1343276  0.5805018  2.4447755 
## 
## Number of Observations: 170
## Number of Groups: 85

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.71695 -0.29570  0.03859  0.33556  1.49385 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  1.00282    0.27832   3.603 0.000416 ***
## GROUP1      -0.02190    0.26455  -0.083 0.934122    
## WAVE        -0.31731    0.11754  -2.700 0.007663 ** 
## BASELINE     0.85518    0.05427  15.758  < 2e-16 ***
## GROUP1:WAVE  0.27064    0.16721   1.619 0.107450    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.545 on 165 degrees of freedom
## Multiple R-squared:  0.6163, Adjusted R-squared:  0.607 
## F-statistic: 66.25 on 4 and 165 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 3.79 0.93   3.92    3.86 1.11 1.79 5.04  3.25 -0.51
##          kurtosis   se
## BASELINE    -0.90 0.08
## meanLET     -0.85 0.10
## -------------------------------------------------------- 
## group: 1
##          vars  n mean   sd median trimmed  mad  min  max range  skew
## BASELINE    1 84 3.67 0.84   3.67    3.69 0.76 2.00 5.00  3.00 -0.17
## meanLET     2 84 4.05 0.78   4.15    4.11 0.77 2.12 5.56  3.44 -0.57
##          kurtosis   se
## BASELINE    -0.97 0.09
## meanLET     -0.28 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=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
##   283.3391 305.2897 -134.6696
## 
## Random effects:
##  Formula: ~1 | ID
##         (Intercept)  Residual
## StdDev:   0.1998569 0.4983305
## 
## Fixed effects: meanLET ~ GROUP * WAVE + BASELINE 
##                  Value  Std.Error DF   t-value p-value
## (Intercept)  1.0027758 0.28205360 83  3.555267  0.0006
## GROUP1      -0.0218990 0.24950769 83 -0.087769  0.9303
## WAVE        -0.3173126 0.10908893 82 -2.908752  0.0047
## BASELINE     0.8551930 0.05788937 82 14.772884  0.0000
## GROUP1:WAVE  0.2706393 0.15519447 82  1.743872  0.0849
##  Correlation: 
##             (Intr) GROUP1 WAVE   BASELI
## GROUP1      -0.469                     
## WAVE        -0.580  0.656              
## BASELINE    -0.784  0.041  0.000       
## GROUP1:WAVE  0.413 -0.933 -0.703 -0.007
## 
## Standardized Within-Group Residuals:
##          Min           Q1          Med           Q3          Max 
## -3.087989741 -0.526219330 -0.003624056  0.620868507  2.759472504 
## 
## Number of Observations: 170
## Number of Groups: 85

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 80 3.64 0.85   3.67    3.65 0.86 2.00   5  3.00 -0.11
## meanLET     2 51 4.20 0.73   4.33    4.30 0.74 2.17   5  2.83 -1.00
##          kurtosis   se
## BASELINE    -1.00 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 60 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 64 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
##   161.6138 180.2579 -73.80689
## 
## Random effects:
##  Formula: ~1 | ID
##         (Intercept)  Residual
## StdDev:   0.3672036 0.3595585
## 
## Fixed effects: meanLET ~ GROUP * WAVE + BASELINE 
##                  Value Std.Error DF   t-value p-value
## (Intercept)  1.0734496 0.3422886 63  3.136095  0.0026
## GROUP1      -0.0332918 0.2466452 63 -0.134979  0.8931
## WAVE        -0.0979069 0.1053771 38 -0.929110  0.3587
## BASELINE     0.7886371 0.0788471 63 10.002111  0.0000
## GROUP1:WAVE  0.2396382 0.1552839 38  1.543226  0.1311
##  Correlation: 
##             (Intr) GROUP1 WAVE   BASELI
## GROUP1      -0.358                     
## WAVE        -0.379  0.585              
## BASELINE    -0.876  0.041 -0.050       
## GROUP1:WAVE  0.260 -0.872 -0.678  0.031
## 
## Standardized Within-Group Residuals:
##         Min          Q1         Med          Q3         Max 
## -3.37395454 -0.36961077  0.04445371  0.41109427  1.91612566 
## 
## Number of Observations: 106
## Number of Groups: 66
Table with P-value

|             |       Value|  Std.Error|  DF|     t-value|    p-value|
|:------------|-----------:|----------:|---:|-----------:|----------:|
|(Intercept)  |   1.0734496|  0.3422886|  63|   3.1360947|  0.0026018|
|GROUP1       |  -0.0332918|  0.2466452|  63|  -0.1349786|  0.8930591|
|WAVE         |  -0.0979069|  0.1053771|  38|  -0.9291096|  0.3586970|
|BASELINE     |   0.7886371|  0.0788471|  63|  10.0021105|  0.0000000|
|GROUP1:WAVE  |   0.2396382|  0.1552839|  38|   1.5432261|  0.1310637|

Table with confidence intervals

est. lower upper
(Intercept) 1.0734496 0.4057675 1.7411317
GROUP1 -0.0332918 -0.5144080 0.4478244
WAVE -0.0979069 -0.3061397 0.1103259
BASELINE 0.7886371 0.6348348 0.9424394
GROUP1:WAVE 0.2396382 -0.0672141 0.5464905