Loading the dataset

data.test4 <- read.csv("~/Final_Adult_Study_R_Docs/adult_study011615.csv")
# Load the psych package
library(psych)
data.test4$meanRES <- apply(data.test4[, c("RES1",  "RES2", "RES3")], 1, mean, na.rm = TRUE)
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 meanRES and ID Group and wave from data.test4 and create a new #dataset with only those variables.
data <- data.test4[,c("ID", "GROUP", "wave", "meanRES")]
#Use dcast to cnage from long-format data to wide format data
data <- dcast(data, ID + GROUP ~ wave, mean, value.var = "meanRES")
# 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 meanRES and wave 2 and 3 of meanRES. 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", "meanRES", "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  8
## 
## -- Imputation 2 --
## 
##   1  2  3  4  5  6  7  8  9 10 11 12 13
## 
## -- Imputation 3 --
## 
##   1  2  3  4  5  6
## 
## -- Imputation 4 --
## 
##   1  2  3  4  5  6  7  8  9
## 
## -- Imputation 5 --
## 
##   1  2  3  4  5  6  7  8  9 10
## 
## -- Imputation 6 --
## 
##   1  2  3  4  5  6  7  8
## 
## -- Imputation 7 --
## 
##   1  2  3  4  5  6  7
## 
## -- Imputation 8 --
## 
##   1  2  3  4  5  6  7  8
## 
## -- Imputation 9 --
## 
##   1  2  3  4  5  6  7  8  9 10 11
## 
## -- Imputation 10 --
## 
##   1  2  3  4  5  6  7  8  9 10
## 
## -- Imputation 11 --
## 
##   1  2  3  4  5  6  7  8  9 10
## 
## -- Imputation 12 --
## 
##   1  2  3  4  5  6  7  8  9
## 
## -- Imputation 13 --
## 
##   1  2  3  4  5  6  7  8  9
## 
## -- Imputation 14 --
## 
##   1  2  3  4  5  6  7  8  9 10
## 
## -- Imputation 15 --
## 
##   1  2  3  4  5  6  7  8  9 10 11 12 13
## 
## -- Imputation 16 --
## 
##   1  2  3  4  5  6  7  8  9
## 
## -- Imputation 17 --
## 
##   1  2  3  4  5  6  7  8  9
## 
## -- Imputation 18 --
## 
##   1  2  3  4  5  6  7  8  9 10
## 
## -- Imputation 19 --
## 
##   1  2  3  4  5  6  7  8  9
## 
## -- Imputation 20 --
## 
##   1  2  3  4  5  6
## 
## -- Imputation 21 --
## 
##   1  2  3  4  5  6  7  8  9
## 
## -- Imputation 22 --
## 
##   1  2  3  4  5  6  7  8  9
## 
## -- Imputation 23 --
## 
##   1  2  3  4  5  6  7  8
## 
## -- Imputation 24 --
## 
##   1  2  3  4  5  6  7  8  9 10
## 
## -- Imputation 25 --
## 
##   1  2  3  4  5  6  7  8  9 10
## 
## -- Imputation 26 --
## 
##   1  2  3  4  5  6  7  8  9 10 11
## 
## -- Imputation 27 --
## 
##   1  2  3  4  5  6  7
## 
## -- Imputation 28 --
## 
##   1  2  3  4  5  6  7  8  9
## 
## -- Imputation 29 --
## 
##   1  2  3  4  5  6  7  8
## 
## -- Imputation 30 --
## 
##   1  2  3  4  5  6  7  8  9 10
## 
## -- Imputation 31 --
## 
##   1  2  3  4  5  6  7  8  9 10
## 
## -- Imputation 32 --
## 
##   1  2  3  4  5  6  7  8  9
## 
## -- Imputation 33 --
## 
##   1  2  3  4  5  6  7
## 
## -- Imputation 34 --
## 
##   1  2  3  4  5  6  7  8  9
## 
## -- Imputation 35 --
## 
##   1  2  3  4  5  6  7  8
## 
## -- Imputation 36 --
## 
##   1  2  3  4  5  6  7  8
## 
## -- Imputation 37 --
## 
##   1  2  3  4  5  6  7  8  9 10
## 
## -- Imputation 38 --
## 
##   1  2  3  4  5  6  7  8
## 
## -- Imputation 39 --
## 
##   1  2  3  4  5  6  7
## 
## -- Imputation 40 --
## 
##   1  2  3  4  5  6  7  8  9 10
## 
## -- Imputation 41 --
## 
##   1  2  3  4  5  6  7  8  9 10 11
## 
## -- Imputation 42 --
## 
##   1  2  3  4  5  6  7  8
## 
## -- Imputation 43 --
## 
##   1  2  3  4  5  6  7  8
## 
## -- Imputation 44 --
## 
##   1  2  3  4  5  6  7  8  9 10 11
## 
## -- Imputation 45 --
## 
##   1  2  3  4  5  6  7  8  9
## 
## -- Imputation 46 --
## 
##   1  2  3  4  5  6  7  8  9 10 11
## 
## -- Imputation 47 --
## 
##   1  2  3  4  5  6  7  8  9
## 
## -- Imputation 48 --
## 
##   1  2  3  4  5  6  7  8  9 10 11 12
## 
## -- Imputation 49 --
## 
##   1  2  3  4  5  6  7  8  9
## 
## -- Imputation 50 --
## 
##   1  2  3  4  5  6  7  8  9

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(meanRES ~ 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.0271050 0.66072045 -0.2738149 2.3280249     43 %
## GROUP1      -0.2135718 0.59871930 -1.3884113 0.9612676     22 %
## WAVE         0.2632860 0.28685338 -0.3002616 0.8268336     31 %
## BASELINE     0.8248579 0.08035961  0.6658420 0.9838738     63 %
## GROUP1:WAVE  0.1837799 0.39331851 -0.5886972 0.9562571     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(meanRES ~ 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.0241802 0.66779766  1.5336684 1.262473e-01
## GROUP1      -0.2136624 0.62291450 -0.3430044 7.316545e-01
## WAVE         0.2632860 0.29912069  0.8801998 3.790986e-01
## BASELINE     0.8253065 0.07946188 10.3861943 1.653024e-18
## GROUP1:WAVE  0.1837680 0.41062360  0.4475339 6.546282e-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(meanRES ~ 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 
## -4.3325 -0.7106  0.0322  0.8306  3.0941 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  1.34557    0.54118   2.486   0.0139 *  
## GROUP1      -0.08889    0.59425  -0.150   0.8813    
## WAVE         0.53166    0.27014   1.968   0.0507 .  
## BASELINE     0.74068    0.05097  14.531   <2e-16 ***
## GROUP1:WAVE -0.08614    0.37576  -0.229   0.8190    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.253 on 173 degrees of freedom
## Multiple R-squared:  0.5579, Adjusted R-squared:  0.5476 
## F-statistic: 54.57 on 4 and 173 DF,  p-value: < 2.2e-16

Describe the meanRES 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 6.52 2.01   7.00    6.59 1.98 2.67 10.00  7.33 -0.38
## meanRES     2 86 6.97 1.99   7.33    7.06 0.99 0.67 12.13 11.46 -0.56
##          kurtosis   se
## BASELINE    -0.94 0.22
## meanRES      0.73 0.21
## -------------------------------------------------------- 
## group: 1
##          vars  n mean   sd median trimmed  mad  min   max range  skew
## BASELINE    1 92 6.81 1.70   7.00    6.90 1.48 2.67 10.00  7.33 -0.48
## meanRES     2 92 6.97 1.74   7.05    7.09 1.59 2.33 10.36  8.03 -0.49
##          kurtosis   se
## BASELINE    -0.41 0.18
## meanRES     -0.40 0.18

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

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

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

Comparing Basline to Wave 2 and 3 by Group.

fullModeldata1 <- lme(meanRES ~ 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
##   590.9248 613.1973 -288.4624
## 
## Random effects:
##  Formula: ~1 | ID
##         (Intercept) Residual
## StdDev:   0.5404128 1.110359
## 
## Fixed effects: meanRES ~ GROUP * WAVE + BASELINE 
##                  Value Std.Error DF   t-value p-value
## (Intercept)  1.3555209 0.5346814 87  2.535194  0.0130
## GROUP1      -0.0885296 0.5468859 87 -0.161879  0.8718
## WAVE         0.5316604 0.2429020 86  2.188786  0.0313
## BASELINE     0.7391517 0.0556006 86 13.293955  0.0000
## GROUP1:WAVE -0.0860835 0.3378739 86 -0.254780  0.7995
##  Correlation: 
##             (Intr) GROUP1 WAVE   BASELI
## GROUP1      -0.512                     
## WAVE        -0.681  0.666              
## BASELINE    -0.678 -0.024  0.000       
## GROUP1:WAVE  0.494 -0.927 -0.719 -0.006
## 
## Standardized Within-Group Residuals:
##         Min          Q1         Med          Q3         Max 
## -3.14617133 -0.52981052  0.01388978  0.56475712  2.63406133 
## 
## Number of Observations: 178
## Number of Groups: 89

Another random imputation

data10=MI$imputations[[10]]
library("Zelig")
zelig.fitdata10 <- zelig(meanRES ~ 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 
## -4.6235 -0.6328  0.1470  0.6860  2.8041 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  1.22719    0.49072   2.501   0.0133 *  
## GROUP1      -0.11339    0.53661  -0.211   0.8329    
## WAVE         0.52643    0.24392   2.158   0.0323 *  
## BASELINE     0.76094    0.04654  16.350   <2e-16 ***
## GROUP1:WAVE  0.06512    0.33929   0.192   0.8480    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.131 on 173 degrees of freedom
## Multiple R-squared:  0.6173, Adjusted R-squared:  0.6084 
## F-statistic: 69.76 on 4 and 173 DF,  p-value: < 2.2e-16

Describe the meanRES 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 6.52 2.01   7.00    6.59 1.98 2.67  10  7.33 -0.38
## meanRES     2 86 6.98 1.85   7.36    7.11 1.22 0.67  10  9.33 -0.83
##          kurtosis   se
## BASELINE    -0.94 0.22
## meanRES      0.63 0.20
## -------------------------------------------------------- 
## group: 1
##          vars  n mean   sd median trimmed  mad  min   max range  skew
## BASELINE    1 92 6.77 1.65   7.00    6.86 1.48 2.67 10.00  7.33 -0.53
## meanRES     2 92 7.15 1.77   7.63    7.27 1.71 2.33 11.22  8.89 -0.54
##          kurtosis   se
## BASELINE    -0.33 0.17
## meanRES     -0.24 0.18

Create a plot that visualizes meanRES variable by the GROUP variable

library(ggplot2)
library(influence.ME)

Take a look at the residuals

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

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

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

Comparing Basline to Wave 2 and 3 by Group.

fullModeldata10 <- lme(meanRES ~ 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
##   556.7672 579.0397 -271.3836
## 
## Random effects:
##  Formula: ~1 | ID
##         (Intercept) Residual
## StdDev:   0.3739858 1.050431
## 
## Fixed effects: meanRES ~ GROUP * WAVE + BASELINE 
##                  Value Std.Error DF   t-value p-value
## (Intercept)  1.2274227 0.4875767 87  2.517394  0.0137
## GROUP1      -0.1133764 0.5119261 87 -0.221470  0.8252
## WAVE         0.5264335 0.2297921 86  2.290912  0.0244
## BASELINE     0.7609037 0.0490780 86 15.503958  0.0000
## GROUP1:WAVE  0.0651190 0.3196338 86  0.203730  0.8390
##  Correlation: 
##             (Intr) GROUP1 WAVE   BASELI
## GROUP1      -0.525                     
## WAVE        -0.707  0.673              
## BASELINE    -0.656 -0.026  0.000       
## GROUP1:WAVE  0.506 -0.937 -0.719  0.003
## 
## Standardized Within-Group Residuals:
##        Min         Q1        Med         Q3        Max 
## -3.8570766 -0.4974058  0.1049904  0.5687370  2.6181529 
## 
## Number of Observations: 178
## Number of Groups: 89

Another random imputation

data15=MI$imputations[[15]]
library("Zelig")
zelig.fitdata15 <- zelig(meanRES ~ 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 
## -4.6359 -0.7149  0.1456  0.6815  3.8616 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   1.8262     0.5687   3.211  0.00158 ** 
## GROUP1       -0.4035     0.6222  -0.648  0.51753    
## WAVE          0.1522     0.2828   0.538  0.59122    
## BASELINE      0.7344     0.0539  13.626  < 2e-16 ***
## GROUP1:WAVE   0.3839     0.3934   0.976  0.33053    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.311 on 173 degrees of freedom
## Multiple R-squared:  0.5276, Adjusted R-squared:  0.5167 
## F-statistic: 48.31 on 4 and 173 DF,  p-value: < 2.2e-16

Describe the meanRES 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 6.52 2.01   7.00    6.59 1.98 2.67 10.00  7.33 -0.38
## meanRES     2 86 6.84 1.91   7.04    6.93 1.43 0.67 11.41 10.75 -0.56
##          kurtosis   se
## BASELINE    -0.94 0.22
## meanRES      0.45 0.21
## -------------------------------------------------------- 
## group: 1
##          vars  n mean   sd median trimmed  mad  min  max range  skew
## BASELINE    1 92 6.77 1.66   7.00    6.86 1.48 2.67 10.0  7.33 -0.52
## meanRES     2 92 7.20 1.86   7.46    7.31 1.96 1.44 10.4  8.97 -0.61
##          kurtosis   se
## BASELINE    -0.36 0.17
## meanRES      0.07 0.19

Create a plot that visualizes meanRES variable by the GROUP variable

library(ggplot2)
library(influence.ME)

Take a look at the residuals

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

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

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

Comparing Basline to Wave 2 and 3 by Group.

fullModeldata15 <- lme(meanRES ~ 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
##   610.5805 632.853 -298.2902
## 
## Random effects:
##  Formula: ~1 | ID
##         (Intercept) Residual
## StdDev:   0.1675616 1.282038
## 
## Fixed effects: meanRES ~ GROUP * WAVE + BASELINE 
##                  Value Std.Error DF   t-value p-value
## (Intercept)  1.8260974 0.5682031 87  3.213811  0.0018
## GROUP1      -0.4034717 0.6179857 87 -0.652882  0.5156
## WAVE         0.1521878 0.2804585 86  0.542639  0.5888
## BASELINE     0.7344467 0.0543483 86 13.513710  0.0000
## GROUP1:WAVE  0.3839045 0.3901089 86  0.984096  0.3278
##  Correlation: 
##             (Intr) GROUP1 WAVE   BASELI
## GROUP1      -0.550                     
## WAVE        -0.740  0.681              
## BASELINE    -0.624 -0.020  0.000       
## GROUP1:WAVE  0.534 -0.947 -0.719 -0.002
## 
## Standardized Within-Group Residuals:
##        Min         Q1        Med         Q3        Max 
## -3.5421682 -0.5524806  0.1054307  0.5113592  2.9833331 
## 
## Number of Observations: 178
## Number of Groups: 89

Another random imputation

data25=MI$imputations[[25]]

library("Zelig")
zelig.fitdata25 <- zelig(meanRES ~ 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 
## -4.2226 -0.7073  0.0577  0.7612  3.2153 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  1.72056    0.50929   3.378 0.000901 ***
## GROUP1      -0.37840    0.55798  -0.678 0.498579    
## WAVE         0.05181    0.25368   0.204 0.838399    
## BASELINE     0.77125    0.04814  16.021  < 2e-16 ***
## GROUP1:WAVE  0.18907    0.35286   0.536 0.592773    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.176 on 173 degrees of freedom
## Multiple R-squared:  0.5988, Adjusted R-squared:  0.5895 
## F-statistic: 64.55 on 4 and 173 DF,  p-value: < 2.2e-16

Describe the meanRES 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 6.52 2.01   7.00    6.59 1.98 2.67 10.00  7.33 -0.38
## meanRES     2 86 6.83 1.82   7.01    6.93 1.47 0.67 10.27  9.60 -0.69
##          kurtosis   se
## BASELINE    -0.94 0.22
## meanRES      0.59 0.20
## -------------------------------------------------------- 
## group: 1
##          vars  n mean   sd median trimmed  mad  min max range  skew
## BASELINE    1 92 6.76 1.67      7    6.84 1.48 2.67  10  7.33 -0.45
## meanRES     2 92 6.92 1.86      7    7.02 1.98 2.33  10  7.67 -0.41
##          kurtosis   se
## BASELINE    -0.37 0.17
## meanRES     -0.77 0.19

Create a plot that visualizes meanRES variable by the GROUP variable

library(ggplot2)
library(influence.ME)

Take a look at the residuals

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

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

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

Comparing Basline to Wave 2 and 3 by Group.

fullModeldata25 <- lme(meanRES ~ 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
##   571.0713 593.3438 -278.5356
## 
## Random effects:
##  Formula: ~1 | ID
##         (Intercept) Residual
## StdDev:   0.3555108 1.103773
## 
## Fixed effects: meanRES ~ GROUP * WAVE + BASELINE 
##                  Value Std.Error DF   t-value p-value
## (Intercept)  1.7260087 0.5062630 87  3.409312  0.0010
## GROUP1      -0.3782460 0.5366099 87 -0.704881  0.4828
## WAVE         0.0518131 0.2414614 86  0.214581  0.8306
## BASELINE     0.7704099 0.0502963 86 15.317413  0.0000
## GROUP1:WAVE  0.1890996 0.3358694 86  0.563015  0.5749
##  Correlation: 
##             (Intr) GROUP1 WAVE   BASELI
## GROUP1      -0.536                     
## WAVE        -0.715  0.675              
## BASELINE    -0.648 -0.017  0.000       
## GROUP1:WAVE  0.518 -0.939 -0.719 -0.006
## 
## Standardized Within-Group Residuals:
##         Min          Q1         Med          Q3         Max 
## -3.42041126 -0.57931853  0.04822548  0.59758321  2.83428234 
## 
## Number of Observations: 178
## Number of Groups: 89

Check assumptions on model without any imputations

Describe the meanRES 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 6.52 2.01   7.00    6.59 1.98 2.67  10  7.33 -0.38
## meanRES     2 59 6.95 1.90   7.33    7.11 0.99 0.67  10  9.33 -0.97
##          kurtosis   se
## BASELINE    -0.94 0.22
## meanRES      0.95 0.25
## -------------------------------------------------------- 
## group: 1
##          vars  n mean   sd median trimmed  mad  min max range  skew
## BASELINE    1 88 6.73 1.68      7    6.81 1.48 2.67  10  7.33 -0.48
## meanRES     2 54 7.31 1.73      8    7.48 1.48 2.33  10  7.67 -0.86
##          kurtosis   se
## BASELINE    -0.42 0.18
## meanRES      0.16 0.24

Create a plot that visualizes meanRES variable by the GROUP variable

library(ggplot2)
library(influence.ME)

Take a look at the residuals

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

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

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

Comparing Basline to Wave 2 and 3 by Group.

fullModel <- lme(meanRES ~ 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
##   355.5365 374.3759 -170.7682
## 
## Random effects:
##  Formula: ~1 | ID
##         (Intercept) Residual
## StdDev:   0.5619731 1.027192
## 
## Fixed effects: meanRES ~ GROUP * WAVE + BASELINE 
##                  Value Std.Error DF   t-value p-value
## (Intercept)  1.0216798 0.6513498 66  1.568558  0.1215
## GROUP1      -0.2956335 0.6435147 66 -0.459404  0.6475
## WAVE         0.1958078 0.2882305 38  0.679345  0.5010
## BASELINE     0.8389021 0.0733463 66 11.437555  0.0000
## GROUP1:WAVE  0.2415480 0.4232303 38  0.570725  0.5715
##  Correlation: 
##             (Intr) GROUP1 WAVE   BASELI
## GROUP1      -0.406                     
## WAVE        -0.602  0.625              
## BASELINE    -0.744 -0.061 -0.018       
## GROUP1:WAVE  0.402 -0.922 -0.681  0.023
## 
## Standardized Within-Group Residuals:
##        Min         Q1        Med         Q3        Max 
## -3.3505420 -0.3916260  0.0791085  0.4685979  2.7998940 
## 
## Number of Observations: 109
## Number of Groups: 69
Table with P-value

|             |       Value|  Std.Error|  DF|     t-value|    p-value|
|:------------|-----------:|----------:|---:|-----------:|----------:|
|(Intercept)  |   1.0216798|  0.6513498|  66|   1.5685578|  0.1215333|
|GROUP1       |  -0.2956335|  0.6435147|  66|  -0.4594044|  0.6474534|
|WAVE         |   0.1958078|  0.2882305|  38|   0.6793447|  0.5010368|
|BASELINE     |   0.8389021|  0.0733463|  66|  11.4375551|  0.0000000|
|GROUP1:WAVE  |   0.2415480|  0.4232303|  38|   0.5707248|  0.5715460|

Table with confidence intervals

est. lower upper
(Intercept) 1.0216798 -0.2486048 2.2919644
GROUP1 -0.2956335 -1.5506378 0.9593708
WAVE 0.1958078 -0.3741443 0.7657600
BASELINE 0.8389021 0.6958596 0.9819445
GROUP1:WAVE 0.2415480 -0.5953552 1.0784513