Loading the dataset

data.test4 <- read.csv("~/Dropbox/ADULT STUDY/adult_study011615.csv")
# Load the psych package
library(psych)
items <- c("PWB1", "PWB2", "PWB3", "PWB4", "PWB5", "PWB6", "PWB7","PWB8", "PWB9")
scaleKey <- c(-1, -1, -1, -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
## 
## -- Imputation 2 --
## 
##   1  2  3  4  5  6  7  8  9 10 11
## 
## -- Imputation 3 --
## 
##   1  2  3  4  5  6  7  8  9 10
## 
## -- Imputation 4 --
## 
##   1  2  3  4  5  6  7  8
## 
## -- Imputation 5 --
## 
##   1  2  3  4  5  6  7  8
## 
## -- Imputation 6 --
## 
##   1  2  3  4  5  6  7  8  9 10 11 12 13 14
## 
## -- Imputation 7 --
## 
##   1  2  3  4  5  6  7
## 
## -- 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
## 
## -- Imputation 10 --
## 
##   1  2  3  4  5  6  7  8
## 
## -- Imputation 11 --
## 
##   1  2  3  4  5  6  7  8  9
## 
## -- Imputation 12 --
## 
##   1  2  3  4  5  6  7  8  9 10 11
## 
## -- Imputation 13 --
## 
##   1  2  3  4  5  6  7  8  9 10 11 12
## 
## -- Imputation 14 --
## 
##   1  2  3  4  5  6  7  8  9
## 
## -- Imputation 15 --
## 
##   1  2  3  4  5  6
## 
## -- Imputation 16 --
## 
##   1  2  3  4  5  6  7  8  9
## 
## -- Imputation 17 --
## 
##   1  2  3  4  5  6  7  8  9 10 11 12 13
## 
## -- Imputation 18 --
## 
##   1  2  3  4  5  6  7  8  9 10
## 
## -- Imputation 19 --
## 
##   1  2  3  4  5  6  7  8  9 10
## 
## -- Imputation 20 --
## 
##   1  2  3  4  5  6  7  8  9 10 11
## 
## -- 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 10
## 
## -- 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
## 
## -- 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
## 
## -- Imputation 29 --
## 
##   1  2  3  4  5  6  7  8  9 10 11
## 
## -- Imputation 30 --
## 
##   1  2  3  4  5  6  7  8  9
## 
## -- Imputation 31 --
## 
##   1  2  3  4  5  6  7  8  9 10 11 12 13 14
## 
## -- Imputation 32 --
## 
##   1  2  3  4  5  6  7  8  9 10
## 
## -- Imputation 33 --
## 
##   1  2  3  4  5  6  7  8  9 10
## 
## -- Imputation 34 --
## 
##   1  2  3  4  5  6
## 
## -- Imputation 35 --
## 
##   1  2  3  4  5  6  7
## 
## -- Imputation 36 --
## 
##   1  2  3  4  5  6  7
## 
## -- Imputation 37 --
## 
##   1  2  3  4  5  6
## 
## -- Imputation 38 --
## 
##   1  2  3  4  5  6  7
## 
## -- Imputation 39 --
## 
##   1  2  3  4  5  6  7  8  9 10
## 
## -- Imputation 40 --
## 
##   1  2  3  4  5  6  7  8  9
## 
## -- Imputation 41 --
## 
##   1  2  3  4  5  6
## 
## -- Imputation 42 --
## 
##   1  2  3  4  5  6  7  8
## 
## -- Imputation 43 --
## 
##   1  2  3  4  5  6  7  8
## 
## -- Imputation 44 --
## 
##   1  2  3  4  5  6  7  8  9 10 11 12 13
## 
## -- Imputation 45 --
## 
##   1  2  3  4  5  6  7  8  9 10
## 
## -- Imputation 46 --
## 
##   1  2  3  4  5  6  7  8  9
## 
## -- Imputation 47 --
## 
##   1  2  3  4  5  6  7  8  9
## 
## -- Imputation 48 --
## 
##   1  2  3  4  5  6  7  8  9
## 
## -- Imputation 49 --
## 
##   1  2  3  4  5  6  7  8
## 
## -- 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.7533391 0.36466672  1.0347943 2.4718840     47 %
## GROUP1      0.1687542 0.27452289 -0.3703624 0.7078707     28 %
## WAVE        0.1179548 0.13588371 -0.1495689 0.3854786     43 %
## BASELINE    0.5994604 0.06787587  0.4655365 0.7333843     52 %
## GROUP1:WAVE 0.1464798 0.18075706 -0.2090116 0.5019713     38 %

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.7471866 0.36517870 4.7844701 3.078826e-06
## GROUP1      0.1690321 0.29656533 0.5699660 5.688538e-01
## WAVE        0.1179548 0.14704561 0.8021650 4.229710e-01
## BASELINE    0.6008504 0.06557747 9.1624519 2.688601e-16
## GROUP1:WAVE 0.1464970 0.19694138 0.7438609 4.573106e-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.65354 -0.35022 -0.00523  0.34284  1.31518 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  1.65998    0.28380   5.849 2.41e-08 ***
## GROUP1       0.07009    0.27256   0.257    0.797    
## WAVE         0.09973    0.12385   0.805    0.422    
## BASELINE     0.62959    0.04640  13.567  < 2e-16 ***
## GROUP1:WAVE  0.18168    0.17228   1.055    0.293    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.5743 on 173 degrees of freedom
## Multiple R-squared:   0.53,  Adjusted R-squared:  0.5191 
## F-statistic: 48.76 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 0.97   4.67    4.47 0.82 2.33 6.00  3.67 -0.43
## meanPWB     2 86 4.60 0.88   4.61    4.63 0.90 2.44 6.47  4.03 -0.29
##          kurtosis   se
## BASELINE    -0.75 0.10
## meanPWB     -0.41 0.09
## -------------------------------------------------------- 
## group: 1
##          vars  n mean   sd median trimmed  mad  min  max range  skew
## BASELINE    1 92 4.19 0.90   4.06    4.24 0.91 1.00 5.56  4.56 -0.76
## meanPWB     2 92 4.79 0.77   4.89    4.82 0.66 2.89 6.71  3.82 -0.36
##          kurtosis   se
## BASELINE     1.53 0.09
## meanPWB     -0.03 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
##   314.9247 337.1972 -150.4624
## 
## Random effects:
##  Formula: ~1 | ID
##         (Intercept)  Residual
## StdDev:    0.210012 0.5257757
## 
## Fixed effects: meanPWB ~ GROUP * WAVE + BASELINE 
##                 Value  Std.Error DF   t-value p-value
## (Intercept) 1.6599404 0.28654535 87  5.792941  0.0000
## GROUP1      0.0700933 0.25717092 87  0.272555  0.7858
## WAVE        0.0997344 0.11501865 86  0.867115  0.3883
## BASELINE    0.6295970 0.04948599 86 12.722733  0.0000
## GROUP1:WAVE 0.1816786 0.15998890 86  1.135570  0.2593
##  Correlation: 
##             (Intr) GROUP1 WAVE   BASELI
## GROUP1      -0.494                     
## WAVE        -0.602  0.671              
## BASELINE    -0.764  0.040  0.000       
## GROUP1:WAVE  0.429 -0.933 -0.719  0.005
## 
## Standardized Within-Group Residuals:
##         Min          Q1         Med          Q3         Max 
## -2.51521653 -0.60939208 -0.01638222  0.60507931  2.34816698 
## 
## 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.66954 -0.26087  0.03462  0.30551  1.19623 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  1.95468    0.26162   7.472 3.76e-12 ***
## GROUP1       0.03325    0.25290   0.131    0.896    
## WAVE         0.03553    0.11495   0.309    0.758    
## BASELINE     0.58731    0.04251  13.815  < 2e-16 ***
## GROUP1:WAVE  0.19472    0.15989   1.218    0.225    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.533 on 173 degrees of freedom
## Multiple R-squared:  0.5374, Adjusted R-squared:  0.5267 
## F-statistic: 50.24 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 0.97   4.67    4.47 0.82 2.33 6.00  3.67 -0.43
## meanPWB     2 86 4.61 0.81   4.69    4.64 0.89 2.44 6.23  3.79 -0.43
##          kurtosis   se
## BASELINE    -0.75 0.10
## meanPWB     -0.31 0.09
## -------------------------------------------------------- 
## group: 1
##          vars  n mean   sd median trimmed  mad  min  max range  skew
## BASELINE    1 92 4.22 0.92   4.06    4.25 0.99 1.00 6.12  5.12 -0.66
## meanPWB     2 92 4.81 0.73   4.89    4.88 0.66 2.83 6.00  3.17 -0.83
##          kurtosis   se
## BASELINE     1.34 0.10
## meanPWB      0.36 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=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
##   289.5508 311.8232 -137.7754
## 
## Random effects:
##  Formula: ~1 | ID
##         (Intercept)  Residual
## StdDev:   0.1443651 0.5052341
## 
## Fixed effects: meanPWB ~ GROUP * WAVE + BASELINE 
##                 Value  Std.Error DF   t-value p-value
## (Intercept) 1.9558686 0.26285735 87  7.440799  0.0000
## GROUP1      0.0332093 0.24516091 87  0.135459  0.8926
## WAVE        0.0355349 0.11052497 86  0.321510  0.7486
## BASELINE    0.5870465 0.04407217 86 13.320116  0.0000
## GROUP1:WAVE 0.1947136 0.15374254 86  1.266492  0.2088
##  Correlation: 
##             (Intr) GROUP1 WAVE   BASELI
## GROUP1      -0.503                     
## WAVE        -0.631  0.676              
## BASELINE    -0.742  0.029  0.000       
## GROUP1:WAVE  0.447 -0.940 -0.719  0.009
## 
## Standardized Within-Group Residuals:
##          Min           Q1          Med           Q3          Max 
## -2.876829997 -0.480030628  0.007902466  0.557933455  2.283222799 
## 
## 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.63138 -0.32913 -0.01152  0.34285  1.54517 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  1.66868    0.27265   6.120 6.07e-09 ***
## GROUP1       0.17888    0.26224   0.682    0.496    
## WAVE         0.08467    0.11919   0.710    0.478    
## BASELINE     0.62938    0.04451  14.139  < 2e-16 ***
## GROUP1:WAVE  0.18306    0.16579   1.104    0.271    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.5526 on 173 degrees of freedom
## Multiple R-squared:  0.5583, Adjusted R-squared:  0.5481 
## F-statistic: 54.66 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 0.97   4.67    4.47 0.82 2.33   6  3.67 -0.43
## meanPWB     2 86 4.58 0.85   4.66    4.63 0.99 2.44   6  3.56 -0.46
##          kurtosis   se
## BASELINE    -0.75 0.10
## meanPWB     -0.59 0.09
## -------------------------------------------------------- 
## group: 1
##          vars  n mean   sd median trimmed  mad  min  max range  skew
## BASELINE    1 92 4.19 0.90   4.06    4.23 0.91 1.00 5.64  4.64 -0.71
## meanPWB     2 92 4.89 0.77   4.89    4.91 0.66 2.44 6.78  4.34 -0.41
##          kurtosis   se
## BASELINE     1.47 0.09
## 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=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
##   294.0226 316.2951 -140.0113
## 
## Random effects:
##  Formula: ~1 | ID
##         (Intercept)  Residual
## StdDev:   0.3028331 0.4529162
## 
## Fixed effects: meanPWB ~ GROUP * WAVE + BASELINE 
##                 Value  Std.Error DF   t-value p-value
## (Intercept) 1.6778529 0.27800569 87  6.035319  0.0000
## GROUP1      0.1785217 0.22761471 87  0.784315  0.4350
## WAVE        0.0846736 0.09907990 86  0.854599  0.3951
## BASELINE    0.6273074 0.05079457 86 12.349893  0.0000
## GROUP1:WAVE 0.1829861 0.13783008 86  1.327621  0.1878
##  Correlation: 
##             (Intr) GROUP1 WAVE   BASELI
## GROUP1      -0.454                     
## WAVE        -0.535  0.653              
## BASELINE    -0.809  0.039  0.000       
## GROUP1:WAVE  0.373 -0.908 -0.719  0.014
## 
## Standardized Within-Group Residuals:
##         Min          Q1         Med          Q3         Max 
## -2.36235446 -0.55729423 -0.04853592  0.56033793  2.73090870 
## 
## 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.74383 -0.27627 -0.00778  0.32377  1.18086 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   1.7350     0.2484   6.985 5.89e-11 ***
## GROUP1        0.1834     0.2405   0.763   0.4467    
## WAVE          0.1822     0.1093   1.667   0.0974 .  
## BASELINE      0.5875     0.0403  14.578  < 2e-16 ***
## GROUP1:WAVE   0.1324     0.1521   0.871   0.3851    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.5069 on 173 degrees of freedom
## Multiple R-squared:  0.5758, Adjusted R-squared:  0.566 
## F-statistic: 58.71 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 0.97   4.67    4.47 0.82 2.33   6  3.67 -0.43
## meanPWB     2 86 4.61 0.78   4.78    4.65 0.83 2.44   6  3.56 -0.51
##          kurtosis   se
## BASELINE    -0.75 0.10
## meanPWB     -0.30 0.08
## -------------------------------------------------------- 
## group: 1
##          vars  n mean   sd median trimmed  mad min  max range  skew
## BASELINE    1 92 4.22 0.93   4.06    4.26 0.95 1.0 6.17  5.17 -0.64
## meanPWB     2 92 4.87 0.74   5.00    4.91 0.51 2.7 6.55  3.85 -0.59
##          kurtosis   se
## BASELINE     1.29 0.10
## meanPWB      0.52 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
##   262.1543 284.4268 -124.0771
## 
## Random effects:
##  Formula: ~1 | ID
##         (Intercept)  Residual
## StdDev:   0.2854322 0.4101631
## 
## Fixed effects: meanPWB ~ GROUP * WAVE + BASELINE 
##                 Value  Std.Error DF   t-value p-value
## (Intercept) 1.7449815 0.25302526 87  6.896471  0.0000
## GROUP1      0.1830895 0.20679706 87  0.885358  0.3784
## WAVE        0.1821764 0.08972723 86  2.030336  0.0454
## BASELINE    0.5852553 0.04626933 86 12.648882  0.0000
## GROUP1:WAVE 0.1323161 0.12481738 86  1.060078  0.2921
##  Correlation: 
##             (Intr) GROUP1 WAVE   BASELI
## GROUP1      -0.450                     
## WAVE        -0.532  0.651              
## BASELINE    -0.809  0.035  0.000       
## GROUP1:WAVE  0.372 -0.905 -0.719  0.013
## 
## Standardized Within-Group Residuals:
##         Min          Q1         Med          Q3         Max 
## -2.70369568 -0.50337479  0.02480347  0.55678331  2.24801498 
## 
## 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 0.97   4.67    4.47 0.82 2.33   6  3.67 -0.43
## meanPWB     2 59 4.69 0.82   4.78    4.75 0.82 2.44   6  3.56 -0.70
##          kurtosis   se
## BASELINE    -0.75 0.10
## meanPWB      0.03 0.11
## -------------------------------------------------------- 
## group: 1
##          vars  n mean   sd median trimmed  mad  min  max range  skew
## BASELINE    1 88 4.17 0.91   4.00    4.21 0.99 1.00 5.56  4.56 -0.69
## meanPWB     2 54 4.88 0.64   4.94    4.93 0.49 2.89 6.00  3.11 -0.83
##          kurtosis   se
## BASELINE     1.41 0.10
## meanPWB      0.69 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=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
##   176.3976 195.237 -81.19878
## 
## Random effects:
##  Formula: ~1 | ID
##         (Intercept)  Residual
## StdDev:   0.4091197 0.3619596
## 
## Fixed effects: meanPWB ~ GROUP * WAVE + BASELINE 
##                 Value Std.Error DF  t-value p-value
## (Intercept) 1.8931604 0.3810992 66 4.967632  0.0000
## GROUP1      0.0601033 0.2504277 66 0.240003  0.8111
## WAVE        0.0681424 0.1066892 38 0.638700  0.5268
## BASELINE    0.5790232 0.0749272 66 7.727815  0.0000
## GROUP1:WAVE 0.2198848 0.1559119 38 1.410314  0.1666
##  Correlation: 
##             (Intr) GROUP1 WAVE   BASELI
## GROUP1      -0.380                     
## WAVE        -0.383  0.585              
## BASELINE    -0.896  0.090 -0.002       
## GROUP1:WAVE  0.258 -0.860 -0.684  0.007
## 
## Standardized Within-Group Residuals:
##         Min          Q1         Med          Q3         Max 
## -2.29290462 -0.47354223  0.08279482  0.52335417  1.88597819 
## 
## Number of Observations: 109
## Number of Groups: 69
Table with P-value

|             |      Value|  Std.Error|  DF|    t-value|    p-value|
|:------------|----------:|----------:|---:|----------:|----------:|
|(Intercept)  |  1.8931604|  0.3810992|  66|  4.9676319|  0.0000051|
|GROUP1       |  0.0601033|  0.2504277|  66|  0.2400027|  0.8110720|
|WAVE         |  0.0681424|  0.1066892|  38|  0.6387000|  0.5268493|
|BASELINE     |  0.5790232|  0.0749272|  66|  7.7278147|  0.0000000|
|GROUP1:WAVE  |  0.2198848|  0.1559119|  38|  1.4103142|  0.1665805|

Table with confidence intervals

est. lower upper
(Intercept) 1.8931604 1.1499278 2.6363930
GROUP1 0.0601033 -0.4282893 0.5484960
WAVE 0.0681424 -0.1428268 0.2791116
BASELINE 0.5790232 0.4328977 0.7251487
GROUP1:WAVE 0.2198848 -0.0884182 0.5281878