Loading the dataset

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

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(meanAPSI ~ 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.49703786 0.26605247  0.97422356 2.0198522     33 %
## GROUP1       0.45718563 0.20043036  0.06377612 0.8505951     24 %
## WAVE         0.08705878 0.09719714 -0.10406545 0.2781830     37 %
## BASELINE     0.60187205 0.05814329  0.48747338 0.7162707     40 %
## GROUP1:WAVE -0.04733202 0.12961342 -0.30196108 0.2072970     31 %

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(meanAPSI ~ 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.49245765 0.26484205  5.6352745 3.090154e-08
## GROUP1       0.45742881 0.21442010  2.1333299 3.312012e-02
## WAVE         0.08705878 0.10441262  0.8337956 4.047977e-01
## BASELINE     0.60302803 0.05607629 10.7537074 1.143795e-22
## GROUP1:WAVE -0.04733942 0.14004871 -0.3380212 7.354471e-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(meanAPSI ~ 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.26130 -0.25303 -0.02228  0.24864  1.09548 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  1.22522    0.21658   5.657 6.26e-08 ***
## GROUP1       0.41811    0.18946   2.207   0.0286 *  
## WAVE         0.10678    0.08608   1.241   0.2165    
## BASELINE     0.66387    0.04252  15.613  < 2e-16 ***
## GROUP1:WAVE -0.06021    0.11973  -0.503   0.6157    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.3991 on 173 degrees of freedom
## Multiple R-squared:  0.5963, Adjusted R-squared:  0.587 
## F-statistic:  63.9 on 4 and 173 DF,  p-value: < 2.2e-16

Describe the meanAPSI 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.96 0.63   4.00    4.00 0.56 2.38   5  2.62 -0.48
## meanAPSI    2 86 4.02 0.63   4.12    4.05 0.60 2.25   5  2.75 -0.49
##          kurtosis   se
## BASELINE    -0.17 0.07
## meanAPSI    -0.31 0.07
## -------------------------------------------------------- 
## group: 1
##          vars  n mean   sd median trimmed  mad  min  max range  skew
## BASELINE    1 92 3.75 0.77   3.75    3.79 0.74 1.62 5.00  3.38 -0.36
## meanAPSI    2 92 4.21 0.61   4.25    4.24 0.64 2.36 5.22  2.86 -0.58
##          kurtosis   se
## BASELINE    -0.11 0.08
## meanAPSI     0.09 0.06

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

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

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

Comparing Basline to Wave 2 and 3 by Group.

fullModeldata1 <- lme(meanAPSI ~ 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
##   179.0357 201.3081 -82.51783
## 
## Random effects:
##  Formula: ~1 | ID
##         (Intercept)  Residual
## StdDev:   0.2134814 0.3305451
## 
## Fixed effects: meanAPSI ~ GROUP * WAVE + BASELINE 
##                  Value  Std.Error DF   t-value p-value
## (Intercept)  1.2406389 0.22524234 87  5.508018  0.0000
## GROUP1       0.4174162 0.16575777 87  2.518230  0.0136
## WAVE         0.1067824 0.07231002 86  1.476731  0.1434
## BASELINE     0.6599836 0.04826536 86 13.674064  0.0000
## GROUP1:WAVE -0.0602796 0.10058481 86 -0.599291  0.5506
##  Correlation: 
##             (Intr) GROUP1 WAVE   BASELI
## GROUP1      -0.424                     
## WAVE        -0.482  0.654              
## BASELINE    -0.849  0.052  0.000       
## GROUP1:WAVE  0.338 -0.910 -0.719  0.009
## 
## Standardized Within-Group Residuals:
##         Min          Q1         Med          Q3         Max 
## -2.85505130 -0.58461252 -0.06396883  0.56136518  2.68866193 
## 
## Number of Observations: 178
## Number of Groups: 89

Another random imputation

data10=MI$imputations[[10]]
library("Zelig")
zelig.fitdata10 <- zelig(meanAPSI ~ 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.31354 -0.26400 -0.03461  0.31217  1.08470 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  1.58732    0.22089   7.186 1.91e-11 ***
## GROUP1       0.55702    0.19430   2.867  0.00466 ** 
## WAVE         0.14824    0.08822   1.680  0.09470 .  
## BASELINE     0.56251    0.04322  13.014  < 2e-16 ***
## GROUP1:WAVE -0.10711    0.12272  -0.873  0.38398    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.4091 on 173 degrees of freedom
## Multiple R-squared:  0.5287, Adjusted R-squared:  0.5178 
## F-statistic: 48.52 on 4 and 173 DF,  p-value: < 2.2e-16

Describe the meanAPSI 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.96 0.63   4.00    4.00 0.56 2.38   5  2.62 -0.48
## meanAPSI    2 86 4.04 0.60   4.07    4.07 0.63 2.25   5  2.75 -0.45
##          kurtosis   se
## BASELINE    -0.17 0.07
## meanAPSI    -0.08 0.06
## -------------------------------------------------------- 
## group: 1
##          vars  n mean   sd median trimmed  mad  min max range  skew
## BASELINE    1 92 3.76 0.78   3.75    3.80 0.93 1.62 5.0  3.38 -0.35
## meanAPSI    2 92 4.32 0.55   4.38    4.36 0.56 2.51 5.5  2.98 -0.59
##          kurtosis   se
## BASELINE    -0.20 0.08
## meanAPSI     0.46 0.06

Create a plot that visualizes meanAPSI variable by the GROUP variable

library(ggplot2)
library(influence.ME)

Take a look at the residuals

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

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

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

Comparing Basline to Wave 2 and 3 by Group.

fullModeldata10 <- lme(meanAPSI ~ 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
##   193.0425 215.315 -89.52126
## 
## Random effects:
##  Formula: ~1 | ID
##         (Intercept)  Residual
## StdDev:   0.1694766 0.3659558
## 
## Fixed effects: meanAPSI ~ GROUP * WAVE + BASELINE 
##                  Value  Std.Error DF   t-value p-value
## (Intercept)  1.5926083 0.22607075 87  7.044734  0.0000
## GROUP1       0.5567035 0.18014172 87  3.090364  0.0027
## WAVE         0.1482438 0.08005647 86  1.851740  0.0675
## BASELINE     0.5611779 0.04680907 86 11.988659  0.0000
## GROUP1:WAVE -0.1070792 0.11136145 86 -0.961547  0.3390
##  Correlation: 
##             (Intr) GROUP1 WAVE   BASELI
## GROUP1      -0.460                     
## WAVE        -0.531  0.667              
## BASELINE    -0.820  0.061  0.000       
## GROUP1:WAVE  0.390 -0.928 -0.719 -0.010
## 
## Standardized Within-Group Residuals:
##         Min          Q1         Med          Q3         Max 
## -2.82811979 -0.63439172 -0.07315838  0.72087545  2.53840665 
## 
## Number of Observations: 178
## Number of Groups: 89

Another random imputation

data15=MI$imputations[[15]]
library("Zelig")
zelig.fitdata15 <- zelig(meanAPSI ~ 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.24566 -0.28418 -0.03953  0.29258  1.27079 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  1.54518    0.22226   6.952 7.06e-11 ***
## GROUP1       0.42640    0.19414   2.196   0.0294 *  
## WAVE         0.04631    0.08820   0.525   0.6002    
## BASELINE     0.60848    0.04368  13.931  < 2e-16 ***
## GROUP1:WAVE -0.01681    0.12269  -0.137   0.8912    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.409 on 173 degrees of freedom
## Multiple R-squared:  0.5526, Adjusted R-squared:  0.5423 
## F-statistic: 53.43 on 4 and 173 DF,  p-value: < 2.2e-16

Describe the meanAPSI 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.96 0.63      4    4.00 0.56 2.38   5  2.62 -0.48
## meanAPSI    2 86 4.03 0.60      4    4.05 0.57 2.25   5  2.75 -0.40
##          kurtosis   se
## BASELINE    -0.17 0.07
## meanAPSI    -0.09 0.06
## -------------------------------------------------------- 
## group: 1
##          vars  n mean   sd median trimmed  mad  min  max range  skew
## BASELINE    1 92 3.75 0.77   3.75    3.78 0.74 1.62 5.00  3.38 -0.34
## meanAPSI    2 92 4.30 0.59   4.37    4.31 0.57 2.83 5.55  2.73 -0.24
##          kurtosis   se
## BASELINE    -0.08 0.08
## meanAPSI    -0.33 0.06

Create a plot that visualizes meanAPSI variable by the GROUP variable

library(ggplot2)
library(influence.ME)

Take a look at the residuals

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

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

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

Comparing Basline to Wave 2 and 3 by Group.

fullModeldata15 <- lme(meanAPSI ~ 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
##   194.1094 216.3819 -90.05472
## 
## Random effects:
##  Formula: ~1 | ID
##         (Intercept)  Residual
## StdDev:   0.1484388 0.3748579
## 
## Fixed effects: meanAPSI ~ GROUP * WAVE + BASELINE 
##                  Value  Std.Error DF   t-value p-value
## (Intercept)  1.5489679 0.22653638 87  6.837612  0.0000
## GROUP1       0.4262227 0.18335647 87  2.324558  0.0224
## WAVE         0.0463059 0.08200389 86  0.564679  0.5738
## BASELINE     0.6075265 0.04652364 86 13.058448  0.0000
## GROUP1:WAVE -0.0168332 0.11406906 86 -0.147571  0.8830
##  Correlation: 
##             (Intr) GROUP1 WAVE   BASELI
## GROUP1      -0.455                     
## WAVE        -0.543  0.671              
## BASELINE    -0.814  0.047  0.000       
## GROUP1:WAVE  0.383 -0.933 -0.719  0.009
## 
## Standardized Within-Group Residuals:
##        Min         Q1        Med         Q3        Max 
## -2.8540388 -0.5924343 -0.1248252  0.6500930  2.9698728 
## 
## Number of Observations: 178
## Number of Groups: 89

Another random imputation

data25=MI$imputations[[25]]

library("Zelig")
zelig.fitdata25 <- zelig(meanAPSI ~ 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.26903 -0.21839 -0.02437  0.23882  1.05949 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  1.60897    0.19949   8.065 1.17e-13 ***
## GROUP1       0.56878    0.17510   3.248  0.00139 ** 
## WAVE         0.07852    0.07952   0.987  0.32483    
## BASELINE     0.58143    0.03909  14.874  < 2e-16 ***
## GROUP1:WAVE -0.14414    0.11061  -1.303  0.19425    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.3687 on 173 degrees of freedom
## Multiple R-squared:  0.5824, Adjusted R-squared:  0.5728 
## F-statistic: 60.32 on 4 and 173 DF,  p-value: < 2.2e-16

Describe the meanAPSI 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.96 0.63   4.00    4.00 0.56 2.38   5  2.62 -0.48
## meanAPSI    2 86 4.03 0.59   4.08    4.05 0.63 2.25   5  2.75 -0.41
##          kurtosis   se
## BASELINE    -0.17 0.07
## meanAPSI     0.03 0.06
## -------------------------------------------------------- 
## group: 1
##          vars  n mean   sd median trimmed  mad  min max range  skew
## BASELINE    1 92 3.77 0.78   3.75    3.80 0.83 1.62 5.0  3.38 -0.37
## meanAPSI    2 92 4.27 0.52   4.26    4.29 0.58 3.00 5.4  2.40 -0.28
##          kurtosis   se
## BASELINE    -0.14 0.08
## meanAPSI    -0.53 0.05

Create a plot that visualizes meanAPSI variable by the GROUP variable

library(ggplot2)
library(influence.ME)

Take a look at the residuals

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

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

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

Comparing Basline to Wave 2 and 3 by Group.

fullModeldata25 <- lme(meanAPSI ~ 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
##   155.1489 177.4214 -70.57447
## 
## Random effects:
##  Formula: ~1 | ID
##         (Intercept)  Residual
## StdDev:    0.163582 0.3246095
## 
## Fixed effects: meanAPSI ~ GROUP * WAVE + BASELINE 
##                  Value  Std.Error DF   t-value p-value
## (Intercept)  1.6136952 0.20494630 87  7.873746  0.0000
## GROUP1       0.5685211 0.16037256 87  3.545002  0.0006
## WAVE         0.0785161 0.07101154 86  1.105680  0.2719
## BASELINE     0.5802339 0.04279840 86 13.557373  0.0000
## GROUP1:WAVE -0.1441231 0.09877697 86 -1.459076  0.1482
##  Correlation: 
##             (Intr) GROUP1 WAVE   BASELI
## GROUP1      -0.452                     
## WAVE        -0.520  0.664              
## BASELINE    -0.827  0.059  0.000       
## GROUP1:WAVE  0.379 -0.924 -0.719 -0.007
## 
## Standardized Within-Group Residuals:
##         Min          Q1         Med          Q3         Max 
## -3.38292777 -0.60530852 -0.06756387  0.63721907  2.73486235 
## 
## Number of Observations: 178
## Number of Groups: 89

Check assumptions on model without any imputations

Describe the meanAPSI 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.96 0.63   4.00    4.00 0.56 2.38   5  2.62 -0.48
## meanAPSI    2 59 4.09 0.65   4.12    4.14 0.56 2.25   5  2.75 -0.68
##          kurtosis   se
## BASELINE    -0.17 0.07
## meanAPSI     0.03 0.08
## -------------------------------------------------------- 
## group: 1
##          vars  n mean   sd median trimmed  mad  min max range  skew
## BASELINE    1 88 3.74 0.78   3.69    3.77 0.83 1.62   5  3.38 -0.32
## meanAPSI    2 54 4.31 0.48   4.38    4.33 0.56 3.00   5  2.00 -0.42
##          kurtosis   se
## BASELINE    -0.18 0.08
## meanAPSI    -0.56 0.07

Create a plot that visualizes meanAPSI variable by the GROUP variable

library(ggplot2)
library(influence.ME)

Take a look at the residuals

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

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

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

Comparing Basline to Wave 2 and 3 by Group.

fullModel <- lme(meanAPSI ~ 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
##   116.0995 134.9389 -51.04973
## 
## Random effects:
##  Formula: ~1 | ID
##         (Intercept)  Residual
## StdDev:   0.2831168 0.2911441
## 
## Fixed effects: meanAPSI ~ GROUP * WAVE + BASELINE 
##                  Value  Std.Error DF   t-value p-value
## (Intercept)  1.5245895 0.29719767 66  5.129884  0.0000
## GROUP1       0.4979467 0.19532844 66  2.549279  0.0131
## WAVE         0.0928050 0.08500811 38  1.091719  0.2818
## BASELINE     0.5917203 0.06740176 66  8.779004  0.0000
## GROUP1:WAVE -0.0865771 0.12429379 38 -0.696552  0.4903
##  Correlation: 
##             (Intr) GROUP1 WAVE   BASELI
## GROUP1      -0.350                     
## WAVE        -0.338  0.594              
## BASELINE    -0.895  0.055 -0.062       
## GROUP1:WAVE  0.222 -0.877 -0.685  0.052
## 
## Standardized Within-Group Residuals:
##         Min          Q1         Med          Q3         Max 
## -2.21220762 -0.54771380 -0.03644676  0.59967572  2.44696288 
## 
## Number of Observations: 109
## Number of Groups: 69
Table with P-value

|             |       Value|  Std.Error|  DF|     t-value|    p-value|
|:------------|-----------:|----------:|---:|-----------:|----------:|
|(Intercept)  |   1.5245895|  0.2971977|  66|   5.1298836|  0.0000027|
|GROUP1       |   0.4979467|  0.1953284|  66|   2.5492793|  0.0131262|
|WAVE         |   0.0928050|  0.0850081|  38|   1.0917193|  0.2818302|
|BASELINE     |   0.5917203|  0.0674018|  66|   8.7790042|  0.0000000|
|GROUP1:WAVE  |  -0.0865771|  0.1242938|  38|  -0.6965522|  0.4903223|

Table with confidence intervals

est. lower upper
(Intercept) 1.5245895 0.9449844 2.1041945
GROUP1 0.4979467 0.1170106 0.8788829
WAVE 0.0928050 -0.0752916 0.2609016
BASELINE 0.5917203 0.4602711 0.7231696
GROUP1:WAVE -0.0865771 -0.3323579 0.1592037