Loading the dataset

data.test4 <- read.csv("/Volumes/TOSHIBA EXT/Dropbox/ADULT STUDY/adult_study011615.csv")
# Load the psych package
library(psych)
items <- c("PWB1", "PWB2","PWB3","PWB4","PWB5","PWB9")
scaleKey <- c(-1,-1,-1,-1,-1,-1)
data.test4[,items] <- apply(data.test4[,items], 2, as.numeric)
data.test4$meanPWB <- scoreItems(scaleKey, items = data.test4[, items])$score
library(reshape2); library(car); library(Amelia);library(mitools);library(nlme);library(predictmeans)
## 
## Attaching package: 'car'
## 
## The following object is masked from 'package:psych':
## 
##     logit
## 
## Loading required package: Rcpp
## ## 
## ## Amelia II: Multiple Imputation
## ## (Version 1.7.3, built: 2014-11-14)
## ## Copyright (C) 2005-2015 James Honaker, Gary King and Matthew Blackwell
## ## Refer to http://gking.harvard.edu/amelia/ for more information
## ## 
## Loading required package: lme4
## Loading required package: Matrix
## 
## Attaching package: 'lme4'
## 
## The following object is masked from 'package:nlme':
## 
##     lmList
#Remove the meanPWB and ID Group and wave from data.test4 and create a new #dataset with only those variables.
data <- data.test4[,c("ID", "GROUP", "wave", "meanPWB")]
#Use dcast to cnage from long-format data to wide format data
data <- dcast(data, ID + GROUP ~ wave, mean, value.var = "meanPWB")
# 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 meanPWB and wave 2 and 3 of meanPWB. 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", "meanPWB", "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
## 
## -- Imputation 3 --
## 
##   1  2  3  4  5  6  7  8  9 10 11
## 
## -- Imputation 4 --
## 
##   1  2  3  4  5  6  7  8  9 10 11 12
## 
## -- Imputation 5 --
## 
##   1  2  3  4  5  6  7  8  9
## 
## -- Imputation 6 --
## 
##   1  2  3  4  5  6  7  8  9
## 
## -- Imputation 7 --
## 
##   1  2  3  4  5  6  7  8
## 
## -- Imputation 8 --
## 
##   1  2  3  4  5  6
## 
## -- Imputation 9 --
## 
##   1  2  3  4  5  6  7  8
## 
## -- Imputation 10 --
## 
##   1  2  3  4  5  6  7
## 
## -- Imputation 11 --
## 
##   1  2  3  4  5  6  7  8  9 10
## 
## -- Imputation 12 --
## 
##   1  2  3  4  5  6  7
## 
## -- Imputation 13 --
## 
##   1  2  3  4  5  6  7  8  9
## 
## -- Imputation 14 --
## 
##   1  2  3  4  5  6  7  8  9
## 
## -- Imputation 15 --
## 
##   1  2  3  4  5  6  7  8
## 
## -- Imputation 16 --
## 
##   1  2  3  4  5  6
## 
## -- Imputation 17 --
## 
##   1  2  3  4  5  6  7  8  9
## 
## -- Imputation 18 --
## 
##   1  2  3  4  5  6
## 
## -- Imputation 19 --
## 
##   1  2  3  4  5  6  7  8  9
## 
## -- Imputation 20 --
## 
##   1  2  3  4  5  6  7  8  9
## 
## -- Imputation 21 --
## 
##   1  2  3  4  5  6  7
## 
## -- 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  8  9 10 11
## 
## -- Imputation 25 --
## 
##   1  2  3  4  5  6  7  8  9 10 11
## 
## -- Imputation 26 --
## 
##   1  2  3  4  5  6  7  8  9 10 11
## 
## -- Imputation 27 --
## 
##   1  2  3  4  5  6  7  8  9
## 
## -- Imputation 28 --
## 
##   1  2  3  4  5  6  7  8  9 10 11
## 
## -- Imputation 29 --
## 
##   1  2  3  4  5  6  7  8
## 
## -- Imputation 30 --
## 
##   1  2  3  4  5  6  7  8
## 
## -- Imputation 31 --
## 
##   1  2  3  4  5  6  7  8  9
## 
## -- Imputation 32 --
## 
##   1  2  3  4  5  6  7  8  9 10 11
## 
## -- Imputation 33 --
## 
##   1  2  3  4  5  6  7  8  9
## 
## -- Imputation 34 --
## 
##   1  2  3  4  5  6  7  8  9
## 
## -- Imputation 35 --
## 
##   1  2  3  4  5  6  7  8  9 10
## 
## -- Imputation 36 --
## 
##   1  2  3  4  5
## 
## -- 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
## 
## -- Imputation 41 --
## 
##   1  2  3  4  5  6  7  8  9 10
## 
## -- Imputation 42 --
## 
##   1  2  3  4  5  6  7  8  9
## 
## -- Imputation 43 --
## 
##   1  2  3  4  5  6  7  8
## 
## -- Imputation 44 --
## 
##   1  2  3  4  5  6  7
## 
## -- Imputation 45 --
## 
##   1  2  3  4  5  6
## 
## -- Imputation 46 --
## 
##   1  2  3  4  5  6  7  8  9 10
## 
## -- Imputation 47 --
## 
##   1  2  3  4  5  6  7  8
## 
## -- Imputation 48 --
## 
##   1  2  3  4  5  6  7
## 
## -- Imputation 49 --
## 
##   1  2  3  4  5  6  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(meanPWB ~ 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.9716470 0.37047461  1.2428440 2.7004500     39 %
## GROUP1      0.1291271 0.32487138 -0.5089366 0.7671909     29 %
## WAVE        0.1672063 0.14403617 -0.1156290 0.4500416     28 %
## BASELINE    0.5490000 0.06573357  0.4194794 0.6785207     47 %
## GROUP1:WAVE 0.1302287 0.20235431 -0.2672015 0.5276589     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(meanPWB ~ 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.9685921  0.3726396 5.2828317 2.288249e-07
## GROUP1      0.1293501  0.3401852 0.3802343 7.038870e-01
## WAVE        0.1672063  0.1525058 1.0963936 2.732315e-01
## BASELINE    0.5496890  0.0642843 8.5509060 2.522338e-15
## GROUP1:WAVE 0.1302306  0.2140289 0.6084721 5.430634e-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(meanPWB ~ 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.81537 -0.30592  0.06771  0.37360  1.31440 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  2.02893    0.28583   7.098 3.12e-11 ***
## GROUP1       0.45954    0.28371   1.620   0.1071    
## WAVE         0.27730    0.12887   2.152   0.0328 *  
## BASELINE     0.50684    0.04521  11.212  < 2e-16 ***
## GROUP1:WAVE -0.05996    0.17926  -0.334   0.7384    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.5975 on 173 degrees of freedom
## Multiple R-squared:  0.4441, Adjusted R-squared:  0.4312 
## F-statistic: 34.55 on 4 and 173 DF,  p-value: < 2.2e-16

Describe the meanPWB 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 4.43 1.00   4.50    4.47 0.99 2.33 6.00  3.67 -0.28
## meanPWB     2 86 4.69 0.83   4.71    4.74 0.80 2.33 6.03  3.69 -0.45
##          kurtosis   se
## BASELINE    -0.80 0.11
## meanPWB     -0.14 0.09
## -------------------------------------------------------- 
## group: 1
##          vars  n mean   sd median trimmed  mad  min  max range  skew
## BASELINE    1 92 4.11 0.99   4.00    4.16 0.99 1.00 5.67  4.67 -0.57
## meanPWB     2 92 4.90 0.74   4.95    4.95 0.67 2.67 6.43  3.76 -0.67
##          kurtosis   se
## BASELINE     0.59 0.10
## meanPWB      0.41 0.08

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

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

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

Comparing Basline to Wave 2 and 3 by Group.

fullModeldata1 <- lme(meanPWB ~ 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
##   321.0479 343.3204 -153.524
## 
## Random effects:
##  Formula: ~1 | ID
##         (Intercept) Residual
## StdDev:   0.3340004 0.485254
## 
## Fixed effects: meanPWB ~ GROUP * WAVE + BASELINE 
##                  Value  Std.Error DF   t-value p-value
## (Intercept)  2.0321604 0.28950428 87  7.019449  0.0000
## GROUP1       0.4593385 0.24471358 87  1.877046  0.0639
## WAVE         0.2773020 0.10615411 86  2.612258  0.0106
## BASELINE     0.5061116 0.05190576 86  9.750586  0.0000
## GROUP1:WAVE -0.0599798 0.14766429 86 -0.406190  0.6856
##  Correlation: 
##             (Intr) GROUP1 WAVE   BASELI
## GROUP1      -0.483                     
## WAVE        -0.550  0.651              
## BASELINE    -0.795  0.060  0.000       
## GROUP1:WAVE  0.387 -0.904 -0.719  0.010
## 
## Standardized Within-Group Residuals:
##         Min          Q1         Med          Q3         Max 
## -2.95683850 -0.47154910  0.09122617  0.52258548  1.93272260 
## 
## Number of Observations: 178
## Number of Groups: 89

Another random imputation

data10=MI$imputations[[10]]
library("Zelig")
zelig.fitdata10 <- zelig(meanPWB ~ 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.64045 -0.42092  0.02494  0.42795  1.34005 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  1.73704    0.29935   5.803 3.05e-08 ***
## GROUP1       0.23822    0.30011   0.794    0.428    
## WAVE         0.18562    0.13636   1.361    0.175    
## BASELINE     0.58677    0.04683  12.528  < 2e-16 ***
## GROUP1:WAVE  0.06364    0.18969   0.335    0.738    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.6323 on 173 degrees of freedom
## Multiple R-squared:  0.4867, Adjusted R-squared:  0.4748 
## F-statistic:    41 on 4 and 173 DF,  p-value: < 2.2e-16

Describe the meanPWB 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 4.43 1.00   4.50    4.47 0.99 2.33 6.00  3.67 -0.28
## meanPWB     2 86 4.62 0.89   4.66    4.65 0.75 2.33 6.66  4.32 -0.29
##          kurtosis   se
## BASELINE    -0.80 0.11
## meanPWB     -0.27 0.10
## -------------------------------------------------------- 
## group: 1
##          vars  n mean   sd median trimmed  mad  min  max range  skew
## BASELINE    1 92 4.12 1.03   4.00    4.17 0.99 1.00 6.09  5.09 -0.43
## meanPWB     2 92 4.77 0.85   4.88    4.84 0.76 2.39 6.34  3.95 -0.79
##          kurtosis   se
## BASELINE     0.32 0.11
## meanPWB      0.22 0.09

Create a plot that visualizes meanPWB variable by the GROUP variable

library(ggplot2)
library(influence.ME)

Take a look at the residuals

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

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

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

Comparing Basline to Wave 2 and 3 by Group.

fullModeldata10 <- lme(meanPWB ~ 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
##   350.6625 372.935 -168.3313
## 
## Random effects:
##  Formula: ~1 | ID
##         (Intercept)  Residual
## StdDev:   0.1365625 0.6081857
## 
## Fixed effects: meanPWB ~ GROUP * WAVE + BASELINE 
##                 Value  Std.Error DF   t-value p-value
## (Intercept) 1.7420386 0.29968094 87  5.812977  0.0000
## GROUP1      0.2379492 0.29430683 87  0.808507  0.4210
## WAVE        0.1856181 0.13304664 86  1.395135  0.1666
## BASELINE    0.5856377 0.04789882 86 12.226557  0.0000
## GROUP1:WAVE 0.0635785 0.18507927 86  0.343520  0.7320
##  Correlation: 
##             (Intr) GROUP1 WAVE   BASELI
## GROUP1      -0.534                     
## WAVE        -0.666  0.678              
## BASELINE    -0.709  0.039  0.000       
## GROUP1:WAVE  0.469 -0.943 -0.719  0.013
## 
## Standardized Within-Group Residuals:
##         Min          Q1         Med          Q3         Max 
## -2.56804392 -0.64470824  0.03749645  0.62035339  2.17184902 
## 
## Number of Observations: 178
## Number of Groups: 89

Another random imputation

data15=MI$imputations[[15]]
library("Zelig")
zelig.fitdata15 <- zelig(meanPWB ~ 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.63368 -0.30825  0.04471  0.38354  1.46160 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  2.14417    0.28897   7.420 5.05e-12 ***
## GROUP1      -0.23187    0.28854  -0.804   0.4227    
## WAVE         0.07800    0.13107   0.595   0.5526    
## BASELINE     0.53673    0.04542  11.818  < 2e-16 ***
## GROUP1:WAVE  0.41924    0.18232   2.299   0.0227 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.6078 on 173 degrees of freedom
## Multiple R-squared:  0.4822, Adjusted R-squared:  0.4703 
## F-statistic: 40.28 on 4 and 173 DF,  p-value: < 2.2e-16

Describe the meanPWB 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 4.43 1.00   4.50    4.47 0.99 2.33   6  3.67 -0.28
## meanPWB     2 86 4.64 0.82   4.67    4.69 0.79 2.33   6  3.67 -0.55
##          kurtosis   se
## BASELINE    -0.80 0.11
## meanPWB     -0.07 0.09
## -------------------------------------------------------- 
## group: 1
##          vars  n mean   sd median trimmed  mad  min  max range  skew
## BASELINE    1 92 4.13 1.02      4    4.17 0.99 1.00 6.44  5.44 -0.47
## meanPWB     2 92 4.87 0.84      5    4.95 0.63 1.87 7.32  5.46 -0.90
##          kurtosis   se
## BASELINE     0.58 0.11
## meanPWB      1.88 0.09

Create a plot that visualizes meanPWB variable by the GROUP variable

library(ggplot2)
library(influence.ME)

Take a look at the residuals

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

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

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

Comparing Basline to Wave 2 and 3 by Group.

fullModeldata15 <- lme(meanPWB ~ 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
##   336.7923 359.0647 -161.3961
## 
## Random effects:
##  Formula: ~1 | ID
##          (Intercept)  Residual
## StdDev: 5.609596e-05 0.5991647
## 
## Fixed effects: meanPWB ~ GROUP * WAVE + BASELINE 
##                  Value  Std.Error DF   t-value p-value
## (Intercept)  2.1441733 0.28896825 87  7.420100  0.0000
## GROUP1      -0.2318735 0.28854135 87 -0.803606  0.4238
## WAVE         0.0780014 0.13107322 86  0.595098  0.5533
## BASELINE     0.5367327 0.04541515 86 11.818360  0.0000
## GROUP1:WAVE  0.4192383 0.18232099 86  2.299452  0.0239
##  Correlation: 
##             (Intr) GROUP1 WAVE   BASELI
## GROUP1      -0.545                     
## WAVE        -0.680  0.681              
## BASELINE    -0.697  0.043  0.000       
## GROUP1:WAVE  0.485 -0.948 -0.719  0.006
## 
## Standardized Within-Group Residuals:
##         Min          Q1         Med          Q3         Max 
## -2.72659786 -0.51447421  0.07462827  0.64011836  2.43940301 
## 
## Number of Observations: 178
## Number of Groups: 89

Another random imputation

data25=MI$imputations[[25]]

library("Zelig")
zelig.fitdata25 <- zelig(meanPWB ~ 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.71155 -0.32531  0.07105  0.35041  1.32603 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  2.09792    0.27761   7.557  2.3e-12 ***
## GROUP1       0.11795    0.27573   0.428    0.669    
## WAVE         0.08069    0.12524   0.644    0.520    
## BASELINE     0.55199    0.04388  12.579  < 2e-16 ***
## GROUP1:WAVE  0.13237    0.17420   0.760    0.448    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.5807 on 173 degrees of freedom
## Multiple R-squared:  0.4867, Adjusted R-squared:  0.4748 
## F-statistic: 41.01 on 4 and 173 DF,  p-value: < 2.2e-16

Describe the meanPWB 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 4.43 1.00    4.5    4.47 0.99 2.33 6.00  3.67 -0.28
## meanPWB     2 86 4.67 0.84    4.8    4.71 0.78 2.33 6.41  4.07 -0.48
##          kurtosis   se
## BASELINE    -0.80 0.11
## meanPWB      0.07 0.09
## -------------------------------------------------------- 
## group: 1
##          vars  n mean   sd median trimmed  mad  min  max range  skew
## BASELINE    1 92 4.11 1.00   4.00    4.17 0.99 1.00 5.67  4.67 -0.58
## meanPWB     2 92 4.81 0.76   4.99    4.88 0.63 2.67 6.08  3.42 -0.87
##          kurtosis   se
## BASELINE     0.57 0.10
## meanPWB      0.53 0.08

Create a plot that visualizes meanPWB variable by the GROUP variable

library(ggplot2)
library(influence.ME)

Take a look at the residuals

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

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

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

Comparing Basline to Wave 2 and 3 by Group.

fullModeldata25 <- lme(meanPWB ~ 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
##   320.0995 342.372 -153.0498
## 
## Random effects:
##  Formula: ~1 | ID
##         (Intercept)  Residual
## StdDev:   0.1544576 0.5512481
## 
## Fixed effects: meanPWB ~ GROUP * WAVE + BASELINE 
##                 Value  Std.Error DF   t-value p-value
## (Intercept) 2.0983168 0.27844795 87  7.535759  0.0000
## GROUP1      0.1179220 0.26762551 87  0.440623  0.6606
## WAVE        0.0806946 0.12059099 86  0.669159  0.5052
## BASELINE    0.5518971 0.04544581 86 12.144070  0.0000
## GROUP1:WAVE 0.1323686 0.16773964 86  0.789131  0.4322
##  Correlation: 
##             (Intr) GROUP1 WAVE   BASELI
## GROUP1      -0.532                     
## WAVE        -0.650  0.676              
## BASELINE    -0.724  0.050  0.000       
## GROUP1:WAVE  0.464 -0.940 -0.719  0.005
## 
## Standardized Within-Group Residuals:
##        Min         Q1        Med         Q3        Max 
## -2.9308728 -0.5520893  0.0893892  0.5390023  2.2301829 
## 
## Number of Observations: 178
## Number of Groups: 89

Check assumptions on model without any imputations

Describe the meanPWB 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 4.43 1.0   4.50    4.47 0.99 2.33   6  3.67 -0.28
## meanPWB     2 59 4.74 0.8   4.83    4.80 0.74 2.33   6  3.67 -0.70
##          kurtosis   se
## BASELINE    -0.80 0.11
## meanPWB      0.47 0.10
## -------------------------------------------------------- 
## group: 1
##          vars  n mean   sd median trimmed  mad  min  max range  skew
## BASELINE    1 88 4.08 1.00   3.92    4.13 0.86 1.00 5.67  4.67 -0.52
## meanPWB     2 54 4.86 0.75   5.00    4.94 0.49 2.67 6.00  3.33 -1.00
##          kurtosis   se
## BASELINE     0.55 0.11
## meanPWB      0.89 0.10

Create a plot that visualizes meanPWB variable by the GROUP variable

library(ggplot2)
library(influence.ME)

Take a look at the residuals

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

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

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

Comparing Basline to Wave 2 and 3 by Group.

fullModel <- lme(meanPWB ~ 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
##   204.2625 223.102 -95.13127
## 
## Random effects:
##  Formula: ~1 | ID
##         (Intercept)  Residual
## StdDev:   0.4211136 0.4381331
## 
## Fixed effects: meanPWB ~ GROUP * WAVE + BASELINE 
##                  Value Std.Error DF   t-value p-value
## (Intercept)  2.1628728 0.4051872 66  5.337960  0.0000
## GROUP1      -0.0035302 0.2945290 66 -0.011986  0.9905
## WAVE         0.0815437 0.1276001 38  0.639056  0.5266
## BASELINE     0.5306399 0.0764377 66  6.942122  0.0000
## GROUP1:WAVE  0.2161482 0.1866448 38  1.158072  0.2541
##  Correlation: 
##             (Intr) GROUP1 WAVE   BASELI
## GROUP1      -0.422                     
## WAVE        -0.453  0.599              
## BASELINE    -0.872  0.104  0.022       
## GROUP1:WAVE  0.297 -0.879 -0.683  0.000
## 
## Standardized Within-Group Residuals:
##         Min          Q1         Med          Q3         Max 
## -1.96683329 -0.47239535  0.09610513  0.48030483  1.49367351 
## 
## Number of Observations: 109
## Number of Groups: 69
Table with P-value

|             |       Value|  Std.Error|  DF|     t-value|    p-value|
|:------------|-----------:|----------:|---:|-----------:|----------:|
|(Intercept)  |   2.1628728|  0.4051872|  66|   5.3379597|  0.0000012|
|GROUP1       |  -0.0035302|  0.2945290|  66|  -0.0119859|  0.9904730|
|WAVE         |   0.0815437|  0.1276001|  38|   0.6390565|  0.5266198|
|BASELINE     |   0.5306399|  0.0764377|  66|   6.9421218|  0.0000000|
|GROUP1:WAVE  |   0.2161482|  0.1866448|  38|   1.1580723|  0.2540615|

Table with confidence intervals

est. lower upper
(Intercept) 2.1628728 1.3726630 2.9530826
GROUP1 -0.0035302 -0.5779308 0.5708704
WAVE 0.0815437 -0.1707751 0.3338624
BASELINE 0.5306399 0.3815685 0.6797114
GROUP1:WAVE 0.2161482 -0.1529267 0.5852231