Loading the dataset

data.test4 <- read.csv("/Volumes/TOSHIBA EXT/Dropbox/ADULT STUDY/adult_study011615.csv")
# Load the psych package
library(psych)
items <- c("PWB6", "PWB7","PWB8")
scaleKey <- c(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 11 12
## 
## -- Imputation 3 --
## 
##   1  2  3  4  5  6  7  8  9
## 
## -- Imputation 4 --
## 
##   1  2  3  4  5
## 
## -- Imputation 5 --
## 
##   1  2  3  4  5  6  7  8
## 
## -- Imputation 6 --
## 
##   1  2  3  4  5  6  7  8  9 10 11
## 
## -- Imputation 7 --
## 
##   1  2  3  4  5  6
## 
## -- Imputation 8 --
## 
##   1  2  3  4  5  6  7  8  9 10 11 12 13 14 15
## 
## -- Imputation 9 --
## 
##   1  2  3  4  5  6  7  8  9 10
## 
## -- Imputation 10 --
## 
##   1  2  3  4  5  6  7  8  9 10 11
## 
## -- Imputation 11 --
## 
##   1  2  3  4  5  6  7  8
## 
## -- Imputation 12 --
## 
##   1  2  3  4  5  6  7
## 
## -- Imputation 13 --
## 
##   1  2  3  4  5  6  7  8
## 
## -- Imputation 14 --
## 
##   1  2  3  4  5  6  7  8  9 10 11 12
## 
## -- Imputation 15 --
## 
##   1  2  3  4  5  6  7  8
## 
## -- Imputation 16 --
## 
##   1  2  3  4  5  6  7  8  9 10
## 
## -- Imputation 17 --
## 
##   1  2  3  4  5  6  7  8  9 10 11
## 
## -- Imputation 18 --
## 
##   1  2  3  4  5  6  7  8  9 10 11 12 13 14
## 
## -- Imputation 19 --
## 
##   1  2  3  4  5  6  7  8  9 10 11 12
## 
## -- Imputation 20 --
## 
##   1  2  3  4  5  6  7  8
## 
## -- Imputation 21 --
## 
##   1  2  3  4  5  6  7  8  9
## 
## -- Imputation 22 --
## 
##   1  2  3  4  5  6  7  8  9 10
## 
## -- Imputation 23 --
## 
##   1  2  3  4  5  6  7
## 
## -- Imputation 24 --
## 
##   1  2  3  4  5  6  7  8
## 
## -- Imputation 25 --
## 
##   1  2  3  4  5  6  7  8  9 10 11
## 
## -- Imputation 26 --
## 
##   1  2  3  4  5  6  7
## 
## -- Imputation 27 --
## 
##   1  2  3  4  5  6  7  8  9
## 
## -- Imputation 28 --
## 
##   1  2  3  4  5  6  7  8
## 
## -- Imputation 29 --
## 
##   1  2  3  4  5  6  7  8  9
## 
## -- Imputation 30 --
## 
##   1  2  3  4  5  6
## 
## -- Imputation 31 --
## 
##   1  2  3  4  5  6
## 
## -- Imputation 32 --
## 
##   1  2  3  4  5  6  7  8  9 10
## 
## -- Imputation 33 --
## 
##   1  2  3  4  5  6  7  8
## 
## -- Imputation 34 --
## 
##   1  2  3  4  5  6  7
## 
## -- Imputation 35 --
## 
##   1  2  3  4  5  6  7  8  9
## 
## -- Imputation 36 --
## 
##   1  2  3  4
## 
## -- Imputation 37 --
## 
##   1  2  3  4  5  6  7  8  9
## 
## -- Imputation 38 --
## 
##   1  2  3  4  5  6  7  8
## 
## -- Imputation 39 --
## 
##   1  2  3  4  5  6  7  8
## 
## -- Imputation 40 --
## 
##   1  2  3  4  5  6  7
## 
## -- Imputation 41 --
## 
##   1  2  3  4  5  6  7  8  9 10 11 12 13
## 
## -- Imputation 42 --
## 
##   1  2  3  4  5  6  7  8  9 10 11 12
## 
## -- Imputation 43 --
## 
##   1  2  3  4  5  6  7  8  9 10 11 12
## 
## -- Imputation 44 --
## 
##   1  2  3  4  5  6  7  8  9 10 11 12 13 14 15
## 
## -- Imputation 45 --
## 
##   1  2  3  4  5  6  7
## 
## -- Imputation 46 --
## 
##   1  2  3  4  5  6  7  8  9 10 11 12 13
## 
## -- Imputation 47 --
## 
##   1  2  3  4  5  6  7  8  9
## 
## -- Imputation 48 --
## 
##   1  2  3  4  5  6  7
## 
## -- Imputation 49 --
## 
##   1  2  3  4  5  6  7  8
## 
## -- Imputation 50 --
## 
##   1  2  3  4  5  6  7

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.60793415 0.41486694  0.7919748 2.4238935     38 %
## GROUP1      0.34114782 0.39302151 -0.4311993 1.1134950     33 %
## WAVE        0.10411135 0.17898332 -0.2477955 0.4560182     36 %
## BASELINE    0.61195820 0.07935049  0.4552177 0.7686987     57 %
## GROUP1:WAVE 0.04686015 0.24994526 -0.4446116 0.5383319     37 %

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.60130631  0.4182633 3.8284647 1.511604e-04
## GROUP1      0.34121986  0.4184118 0.8155120 4.151093e-01
## WAVE        0.10411135  0.1928079 0.5399746 5.894463e-01
## BASELINE    0.61346082  0.0768415 7.9834570 4.478782e-13
## GROUP1:WAVE 0.04684212  0.2691091 0.1740637 8.618854e-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 
## -2.50633 -0.40505  0.02167  0.47750  2.02231 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  1.70235    0.37577   4.530 1.09e-05 ***
## GROUP1       0.47365    0.39679   1.194    0.234    
## WAVE         0.08302    0.18041   0.460    0.646    
## BASELINE     0.58814    0.05546  10.606  < 2e-16 ***
## GROUP1:WAVE -0.01492    0.25096  -0.059    0.953    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.8365 on 173 degrees of freedom
## Multiple R-squared:  0.4187, Adjusted R-squared:  0.4052 
## F-statistic: 31.15 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.41 1.15   4.67    4.49 0.99 1.67 6.00  4.33 -0.58
## meanPWB     2 86 4.42 1.10   4.64    4.48 1.02 1.67 6.77  5.11 -0.49
##          kurtosis   se
## BASELINE    -0.64 0.12
## meanPWB     -0.29 0.12
## -------------------------------------------------------- 
## group: 1
##          vars  n mean   sd median trimmed  mad  min max range  skew
## BASELINE    1 92 4.38 1.12   4.67    4.41 1.24 1.00 6.0  5.00 -0.55
## meanPWB     2 92 4.85 1.04   5.00    4.92 0.99 1.28 7.4  6.12 -0.75
##          kurtosis   se
## BASELINE    -0.10 0.12
## meanPWB      1.32 0.11

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
##   444.3389 466.6114 -215.1694
## 
## Random effects:
##  Formula: ~1 | ID
##         (Intercept) Residual
## StdDev:   0.4198965 0.709806
## 
## Fixed effects: meanPWB ~ GROUP * WAVE + BASELINE 
##                  Value Std.Error DF   t-value p-value
## (Intercept)  1.7099584 0.3734907 87  4.578316  0.0000
## GROUP1       0.4737210 0.3532608 87  1.340995  0.1834
## WAVE         0.0830250 0.1552771 86  0.534689  0.5942
## BASELINE     0.5864158 0.0620877 86  9.444966  0.0000
## GROUP1:WAVE -0.0150085 0.2160072 86 -0.069481  0.9448
##  Correlation: 
##             (Intr) GROUP1 WAVE   BASELI
## GROUP1      -0.483                     
## WAVE        -0.624  0.659              
## BASELINE    -0.733 -0.007  0.000       
## GROUP1:WAVE  0.438 -0.917 -0.719  0.014
## 
## Standardized Within-Group Residuals:
##         Min          Q1         Med          Q3         Max 
## -3.01369482 -0.55687626  0.08657019  0.63990176  2.01437269 
## 
## 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 
## -2.24249 -0.37987  0.07047  0.43611  1.76890 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  1.89153    0.33363   5.670 5.89e-08 ***
## GROUP1       0.50848    0.35234   1.443   0.1508    
## WAVE         0.34334    0.16019   2.143   0.0335 *  
## BASELINE     0.49910    0.04923  10.138  < 2e-16 ***
## GROUP1:WAVE -0.21119    0.22283  -0.948   0.3446    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.7428 on 173 degrees of freedom
## Multiple R-squared:   0.39,  Adjusted R-squared:  0.3759 
## F-statistic: 27.65 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.41 1.15   4.67    4.49 0.99 1.67 6.00  4.33 -0.58
## meanPWB     2 86 4.61 1.03   4.67    4.70 0.99 1.67 6.14  4.48 -0.79
##          kurtosis   se
## BASELINE    -0.64 0.12
## meanPWB      0.16 0.11
## -------------------------------------------------------- 
## group: 1
##          vars  n mean   sd median trimmed  mad  min  max range  skew
## BASELINE    1 92 4.36 1.12   4.50    4.39 1.24 1.00 6.08  5.08 -0.50
## meanPWB     2 92 4.78 0.84   4.87    4.81 0.98 2.67 6.08  3.42 -0.34
##          kurtosis   se
## BASELINE    -0.11 0.12
## meanPWB     -0.70 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
##   406.2107 428.4832 -196.1053
## 
## Random effects:
##  Formula: ~1 | ID
##         (Intercept)  Residual
## StdDev:   0.2826776 0.6755069
## 
## Fixed effects: meanPWB ~ GROUP * WAVE + BASELINE 
##                  Value Std.Error DF   t-value p-value
## (Intercept)  1.9035673 0.3322516 87  5.729295  0.0000
## GROUP1       0.5081979 0.3306881 87  1.536789  0.1280
## WAVE         0.3433449 0.1477738 86  2.323449  0.0225
## BASELINE     0.4963737 0.0526280 86  9.431743  0.0000
## GROUP1:WAVE -0.2110830 0.2055581 86 -1.026878  0.3074
##  Correlation: 
##             (Intr) GROUP1 WAVE   BASELI
## GROUP1      -0.526                     
## WAVE        -0.667  0.670              
## BASELINE    -0.699  0.017  0.000       
## GROUP1:WAVE  0.486 -0.932 -0.719 -0.010
## 
## Standardized Within-Group Residuals:
##         Min          Q1         Med          Q3         Max 
## -2.98618429 -0.46326207  0.06400816  0.56598691  2.56552143 
## 
## 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 
## -2.18551 -0.33840  0.00679  0.47391  1.82379 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  1.55084    0.32152   4.823 3.08e-06 ***
## GROUP1      -0.02195    0.33965  -0.065    0.949    
## WAVE        -0.02734    0.15443  -0.177    0.860    
## BASELINE     0.67113    0.04742  14.152  < 2e-16 ***
## GROUP1:WAVE  0.34227    0.21482   1.593    0.113    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.7161 on 173 degrees of freedom
## Multiple R-squared:  0.5632, Adjusted R-squared:  0.5531 
## F-statistic: 55.77 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.41 1.15   4.67    4.49 0.99 1.67 6.00  4.33 -0.58
## meanPWB     2 86 4.47 1.16   4.67    4.54 0.99 1.67 6.86  5.19 -0.47
##          kurtosis   se
## BASELINE    -0.64 0.12
## meanPWB     -0.41 0.12
## -------------------------------------------------------- 
## group: 1
##          vars  n mean   sd median trimmed  mad  min  max range  skew
## BASELINE    1 92 4.38 1.13    4.5    4.41 1.23 1.00 6.37  5.37 -0.52
## meanPWB     2 92 4.94 0.94    5.0    5.00 0.99 2.67 6.47  3.81 -0.47
##          kurtosis   se
## BASELINE    -0.04 0.12
## meanPWB     -0.40 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=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
##   389.5325 411.805 -187.7662
## 
## Random effects:
##  Formula: ~1 | ID
##         (Intercept)  Residual
## StdDev:   0.3514034 0.6122592
## 
## Fixed effects: meanPWB ~ GROUP * WAVE + BASELINE 
##                  Value Std.Error DF   t-value p-value
## (Intercept)  1.5645213 0.3194337 87  4.897797  0.0000
## GROUP1      -0.0217530 0.3041370 87 -0.071524  0.9431
## WAVE        -0.0273355 0.1339378 86 -0.204091  0.8388
## BASELINE     0.6680227 0.0527979 86 12.652450  0.0000
## GROUP1:WAVE  0.3420845 0.1863293 86  1.835913  0.0698
##  Correlation: 
##             (Intr) GROUP1 WAVE   BASELI
## GROUP1      -0.484                     
## WAVE        -0.629  0.661              
## BASELINE    -0.729 -0.011  0.000       
## GROUP1:WAVE  0.440 -0.919 -0.719  0.017
## 
## Standardized Within-Group Residuals:
##         Min          Q1         Med          Q3         Max 
## -2.57327354 -0.49390244 -0.01965743  0.58478052  2.41406374 
## 
## 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 
## -2.04944 -0.32486  0.02334  0.41755  1.70181 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  1.74585    0.31024   5.627 7.25e-08 ***
## GROUP1       0.14660    0.32730   0.448    0.655    
## WAVE         0.04578    0.14882   0.308    0.759    
## BASELINE     0.59804    0.04584  13.046  < 2e-16 ***
## GROUP1:WAVE  0.23860    0.20701   1.153    0.251    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.69 on 173 degrees of freedom
## Multiple R-squared:  0.5307, Adjusted R-squared:  0.5198 
## F-statistic:  48.9 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.41 1.15   4.67    4.49 0.99 1.67   6  4.33 -0.58
## meanPWB     2 86 4.45 1.06   4.67    4.52 0.99 1.67   6  4.33 -0.59
##          kurtosis   se
## BASELINE    -0.64 0.12
## meanPWB     -0.36 0.11
## -------------------------------------------------------- 
## group: 1
##          vars  n mean   sd median trimmed  mad  min  max range  skew
## BASELINE    1 92 4.38 1.12   4.67    4.42 1.11 1.00 6.00  5.00 -0.56
## meanPWB     2 92 4.94 0.88   5.00    4.97 0.99 2.67 6.91  4.25 -0.26
##          kurtosis   se
## BASELINE    -0.07 0.12
## meanPWB     -0.55 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=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
##   379.8882 402.1607 -182.9441
## 
## Random effects:
##  Formula: ~1 | ID
##         (Intercept)  Residual
## StdDev:   0.2661288 0.6260737
## 
## Fixed effects: meanPWB ~ GROUP * WAVE + BASELINE 
##                 Value Std.Error DF   t-value p-value
## (Intercept) 1.7472448 0.3092580 87  5.649797  0.0000
## GROUP1      0.1466137 0.3066149 87  0.478169  0.6337
## WAVE        0.0457796 0.1369598 86  0.334256  0.7390
## BASELINE    0.5977225 0.0491769 86 12.154538  0.0000
## GROUP1:WAVE 0.2385862 0.1905170 86  1.252309  0.2139
##  Correlation: 
##             (Intr) GROUP1 WAVE   BASELI
## GROUP1      -0.509                     
## WAVE        -0.664  0.670              
## BASELINE    -0.701 -0.005  0.000       
## GROUP1:WAVE  0.470 -0.932 -0.719  0.011
## 
## Standardized Within-Group Residuals:
##         Min          Q1         Med          Q3         Max 
## -2.79041074 -0.50594274  0.01358961  0.67732521  2.30821606 
## 
## 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.41 1.15   4.67    4.49 0.99 1.67   6  4.33 -0.58
## meanPWB     2 59 4.59 1.08   4.67    4.69 0.99 1.67   6  4.33 -0.83
##          kurtosis   se
## BASELINE    -0.64 0.12
## meanPWB      0.08 0.14
## -------------------------------------------------------- 
## group: 1
##          vars  n mean   sd median trimmed  mad  min max range  skew
## BASELINE    1 88 4.36 1.13    4.5    4.39 1.24 1.00   6  5.00 -0.53
## meanPWB     2 54 4.93 0.80    5.0    4.98 0.99 2.67   6  3.33 -0.61
##          kurtosis   se
## BASELINE    -0.11 0.12
## meanPWB     -0.23 0.11

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
##   249.0491 267.8885 -117.5246
## 
## Random effects:
##  Formula: ~1 | ID
##         (Intercept)  Residual
## StdDev:   0.5145779 0.5396038
## 
## Fixed effects: meanPWB ~ GROUP * WAVE + BASELINE 
##                 Value Std.Error DF  t-value p-value
## (Intercept) 1.7101552 0.4436641 66 3.854617  0.0003
## GROUP1      0.1435482 0.3603730 66 0.398332  0.6917
## WAVE        0.0531790 0.1572023 38 0.338284  0.7370
## BASELINE    0.5962042 0.0832944 66 7.157798  0.0000
## GROUP1:WAVE 0.2218678 0.2297903 38 0.965523  0.3404
##  Correlation: 
##             (Intr) GROUP1 WAVE   BASELI
## GROUP1      -0.392                     
## WAVE        -0.447  0.599              
## BASELINE    -0.835  0.022 -0.048       
## GROUP1:WAVE  0.318 -0.884 -0.683  0.019
## 
## Standardized Within-Group Residuals:
##         Min          Q1         Med          Q3         Max 
## -2.06474848 -0.37527738  0.05915006  0.53436891  1.69783161 
## 
## Number of Observations: 109
## Number of Groups: 69
Table with P-value

|             |      Value|  Std.Error|  DF|    t-value|    p-value|
|:------------|----------:|----------:|---:|----------:|----------:|
|(Intercept)  |  1.7101552|  0.4436641|  66|  3.8546169|  0.0002652|
|GROUP1       |  0.1435482|  0.3603730|  66|  0.3983323|  0.6916709|
|WAVE         |  0.0531790|  0.1572023|  38|  0.3382835|  0.7370115|
|BASELINE     |  0.5962042|  0.0832944|  66|  7.1577982|  0.0000000|
|GROUP1:WAVE  |  0.2218678|  0.2297903|  38|  0.9655229|  0.3403878|

Table with confidence intervals

est. lower upper
(Intercept) 1.7101552 0.8449063 2.5754041
GROUP1 0.1435482 -0.5592635 0.8463600
WAVE 0.0531790 -0.2576758 0.3640337
BASELINE 0.5962042 0.4337607 0.7586477
GROUP1:WAVE 0.2218678 -0.2325237 0.6762592