Loading the dataset

data.test4 <- read.csv("/Volumes/TOSHIBA EXT/Dropbox/ADULT STUDY/adult_study011615.csv")
# Load the psych package
library(psych)
data.test4$LS <- apply(data.test4[, c("LS1","LS2", "LS3", "LS4",  "LS5")], 1, mean, na.rm = TRUE)
library(reshape2); library(car); library(Amelia);library(mitools);library(nlme);library(predictmeans)
## 
## Attaching package: 'car'
## 
## The following object is masked from 'package:psych':
## 
##     logit
## 
## Loading required package: Rcpp
## ## 
## ## Amelia II: Multiple Imputation
## ## (Version 1.7.3, built: 2014-11-14)
## ## Copyright (C) 2005-2015 James Honaker, Gary King and Matthew Blackwell
## ## Refer to http://gking.harvard.edu/amelia/ for more information
## ## 
## Loading required package: lme4
## Loading required package: Matrix
## 
## Attaching package: 'lme4'
## 
## The following object is masked from 'package:nlme':
## 
##     lmList
#Remove the LS and ID Group and wave from data.test4 and create a new #dataset with only those variables.
data <- data.test4[,c("ID", "GROUP", "wave", "LS")]
#Use dcast to cnage from long-format data to wide format data
data <- dcast(data, ID + GROUP ~ wave, mean, value.var = "LS")
# 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 LS and wave 2 and 3 of LS. 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", "LS", "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
## 
## -- Imputation 2 --
## 
##   1  2  3  4  5  6
## 
## -- Imputation 3 --
## 
##   1  2  3  4  5  6  7
## 
## -- Imputation 4 --
## 
##   1  2  3  4  5  6
## 
## -- Imputation 5 --
## 
##   1  2  3  4  5  6  7
## 
## -- Imputation 6 --
## 
##   1  2  3  4  5  6  7  8  9
## 
## -- Imputation 7 --
## 
##   1  2  3  4  5  6  7
## 
## -- Imputation 8 --
## 
##   1  2  3  4  5  6  7  8
## 
## -- Imputation 9 --
## 
##   1  2  3  4  5  6  7  8  9
## 
## -- Imputation 10 --
## 
##   1  2  3  4  5  6  7  8
## 
## -- 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
## 
## -- Imputation 14 --
## 
##   1  2  3  4  5  6  7
## 
## -- Imputation 15 --
## 
##   1  2  3  4  5  6  7
## 
## -- Imputation 16 --
## 
##   1  2  3  4  5  6  7  8  9 10
## 
## -- Imputation 17 --
## 
##   1  2  3  4  5  6  7  8  9 10
## 
## -- Imputation 18 --
## 
##   1  2  3  4  5  6  7
## 
## -- Imputation 19 --
## 
##   1  2  3  4  5  6
## 
## -- Imputation 20 --
## 
##   1  2  3  4  5  6  7
## 
## -- Imputation 21 --
## 
##   1  2  3  4  5  6  7  8
## 
## -- Imputation 22 --
## 
##   1  2  3  4  5  6
## 
## -- Imputation 23 --
## 
##   1  2  3  4  5  6  7  8
## 
## -- Imputation 24 --
## 
##   1  2  3  4  5  6  7  8  9
## 
## -- Imputation 25 --
## 
##   1  2  3  4  5  6  7
## 
## -- Imputation 26 --
## 
##   1  2  3  4  5
## 
## -- Imputation 27 --
## 
##   1  2  3  4  5  6  7
## 
## -- Imputation 28 --
## 
##   1  2  3  4  5  6  7  8
## 
## -- Imputation 29 --
## 
##   1  2  3  4  5  6  7  8  9 10
## 
## -- Imputation 30 --
## 
##   1  2  3  4  5  6  7  8
## 
## -- Imputation 31 --
## 
##   1  2  3  4  5  6  7
## 
## -- Imputation 32 --
## 
##   1  2  3  4  5  6  7  8  9
## 
## -- Imputation 33 --
## 
##   1  2  3  4  5  6  7  8  9 10 11 12
## 
## -- Imputation 34 --
## 
##   1  2  3  4  5  6  7  8  9
## 
## -- Imputation 35 --
## 
##   1  2  3  4  5  6
## 
## -- Imputation 36 --
## 
##   1  2  3  4  5  6  7  8
## 
## -- Imputation 37 --
## 
##   1  2  3  4  5  6  7  8
## 
## -- Imputation 38 --
## 
##   1  2  3  4  5  6  7
## 
## -- Imputation 39 --
## 
##   1  2  3  4  5  6  7  8  9 10 11
## 
## -- Imputation 40 --
## 
##   1  2  3  4  5  6  7  8
## 
## -- Imputation 41 --
## 
##   1  2  3  4  5  6  7  8
## 
## -- Imputation 42 --
## 
##   1  2  3  4  5  6
## 
## -- 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  7  8  9 10
## 
## -- Imputation 46 --
## 
##   1  2  3  4  5  6  7  8  9
## 
## -- Imputation 47 --
## 
##   1  2  3  4  5  6  7
## 
## -- Imputation 48 --
## 
##   1  2  3  4  5  6  7  8  9 10
## 
## -- Imputation 49 --
## 
##   1  2  3  4  5  6  7  8  9 10
## 
## -- 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(LS ~ 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.94828014 0.46676303  1.0305106 2.8660497     36 %
## GROUP1       0.60102583 0.45378707 -0.2900567 1.4921083     28 %
## WAVE         0.02225812 0.21840447 -0.4074096 0.4519258     39 %
## BASELINE     0.58257085 0.07776137  0.4292915 0.7358503     48 %
## GROUP1:WAVE -0.06362630 0.29400908 -0.6415940 0.5143414     35 %

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(LS ~ 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.94104048  0.4763352  4.07494641 5.509874e-05
## GROUP1       0.60078111  0.4929801  1.21867209 2.232898e-01
## WAVE         0.02225812  0.2387719  0.09321916 9.257697e-01
## BASELINE     0.58419562  0.0742848  7.86426841 3.064407e-13
## GROUP1:WAVE -0.06365690  0.3232194 -0.19694640 8.439369e-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(LS ~ 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 
## -3.2231 -0.4806  0.1140  0.5835  2.2379 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  2.17489    0.38485   5.651 6.44e-08 ***
## GROUP1       0.81887    0.43209   1.895   0.0597 .  
## WAVE         0.11712    0.19646   0.596   0.5518    
## BASELINE     0.51906    0.05099  10.180  < 2e-16 ***
## GROUP1:WAVE -0.25202    0.27330  -0.922   0.3577    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.9109 on 173 degrees of freedom
## Multiple R-squared:  0.4099, Adjusted R-squared:  0.3962 
## F-statistic: 30.04 on 4 and 173 DF,  p-value: < 2.2e-16

Describe the LS 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.46 1.29   4.60    4.51 1.19 1.6   7   5.4 -0.30
## LS          2 86 4.66 1.29   4.92    4.76 1.06 1.0   7   6.0 -0.73
##          kurtosis   se
## BASELINE    -0.46 0.14
## LS           0.20 0.14
## -------------------------------------------------------- 
## group: 1
##          vars  n mean   sd median trimmed  mad min  max range  skew
## BASELINE    1 92 4.66 1.40    4.9    4.71 1.33 1.6 8.14  6.54 -0.27
## LS          2 92 5.21 0.98    5.4    5.29 0.89 2.2 7.06  4.86 -0.79
##          kurtosis   se
## BASELINE    -0.50 0.15
## LS           0.37 0.10

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

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

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

Comparing Basline to Wave 2 and 3 by Group.

fullModeldata1 <- lme(LS ~ 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
##   469.4423 491.7148 -227.7212
## 
## Random effects:
##  Formula: ~1 | ID
##         (Intercept)  Residual
## StdDev:   0.5295191 0.7254564
## 
## Fixed effects: LS ~ GROUP * WAVE + BASELINE 
##                  Value Std.Error DF   t-value p-value
## (Intercept)  2.2201630 0.3695873 87  6.007140  0.0000
## GROUP1       0.8195845 0.3671797 87  2.232107  0.0282
## WAVE         0.1171242 0.1587008 86  0.738019  0.4625
## BASELINE     0.5089034 0.0580569 86  8.765604  0.0000
## GROUP1:WAVE -0.2511311 0.2208057 86 -1.137340  0.2586
##  Correlation: 
##             (Intr) GROUP1 WAVE   BASELI
## GROUP1      -0.506                     
## WAVE        -0.644  0.648              
## BASELINE    -0.700 -0.011  0.000       
## GROUP1:WAVE  0.479 -0.901 -0.719 -0.023
## 
## Standardized Within-Group Residuals:
##         Min          Q1         Med          Q3         Max 
## -3.16654934 -0.39550579  0.08800447  0.54275163  2.31747974 
## 
## Number of Observations: 178
## Number of Groups: 89

Another random imputation

data10=MI$imputations[[10]]
library("Zelig")
zelig.fitdata10 <- zelig(LS ~ 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 
## -3.2889 -0.4706 -0.0324  0.5131  2.5327 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  1.75706    0.37818   4.646 6.68e-06 ***
## GROUP1       0.70911    0.42205   1.680   0.0947 .  
## WAVE        -0.03806    0.19187  -0.198   0.8430    
## BASELINE     0.63176    0.05067  12.467  < 2e-16 ***
## GROUP1:WAVE -0.09314    0.26689  -0.349   0.7275    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.8897 on 173 degrees of freedom
## Multiple R-squared:  0.5124, Adjusted R-squared:  0.5011 
## F-statistic: 45.45 on 4 and 173 DF,  p-value: < 2.2e-16

Describe the LS 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.46 1.29   4.60    4.51 1.19 1.60   7  5.40 -0.30
## LS          2 86 4.51 1.33   4.65    4.61 1.12 0.61   7  6.39 -0.73
##          kurtosis   se
## BASELINE    -0.46 0.14
## LS           0.39 0.14
## -------------------------------------------------------- 
## group: 1
##          vars  n mean   sd median trimmed  mad min  max range  skew
## BASELINE    1 92 4.63 1.36    4.8    4.69 1.33 1.6 6.80  5.20 -0.40
## LS          2 92 5.19 1.10    5.4    5.25 0.90 2.2 8.76  6.56 -0.32
##          kurtosis   se
## BASELINE    -0.66 0.14
## LS           0.89 0.11

Create a plot that visualizes LS variable by the GROUP variable

library(ggplot2)
library(influence.ME)

Take a look at the residuals

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

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

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

Comparing Basline to Wave 2 and 3 by Group.

fullModeldata10 <- lme(LS ~ 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
##   465.8144 488.0869 -225.9072
## 
## Random effects:
##  Formula: ~1 | ID
##         (Intercept) Residual
## StdDev:   0.4542161 0.750325
## 
## Fixed effects: LS ~ GROUP * WAVE + BASELINE 
##                  Value Std.Error DF   t-value p-value
## (Intercept)  1.7605248 0.3695967 87  4.763367  0.0000
## GROUP1       0.7092195 0.3740799 87  1.895904  0.0613
## WAVE        -0.0380638 0.1641410 86 -0.231897  0.8172
## BASELINE     0.6309832 0.0569128 86 11.086834  0.0000
## GROUP1:WAVE -0.0931220 0.2283175 86 -0.407862  0.6844
##  Correlation: 
##             (Intr) GROUP1 WAVE   BASELI
## GROUP1      -0.508                     
## WAVE        -0.666  0.658              
## BASELINE    -0.686 -0.022  0.000       
## GROUP1:WAVE  0.482 -0.915 -0.719 -0.005
## 
## Standardized Within-Group Residuals:
##         Min          Q1         Med          Q3         Max 
## -3.33765148 -0.47397296 -0.01404168  0.58964947  2.85613991 
## 
## Number of Observations: 178
## Number of Groups: 89

Another random imputation

data15=MI$imputations[[15]]
library("Zelig")
zelig.fitdata15 <- zelig(LS ~ 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 
## -3.1660 -0.5311  0.0510  0.6116  2.6180 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  2.10077    0.42265   4.971  1.6e-06 ***
## GROUP1       0.79612    0.47266   1.684   0.0939 .  
## WAVE        -0.07060    0.21490  -0.329   0.7429    
## BASELINE     0.55388    0.05641   9.819  < 2e-16 ***
## GROUP1:WAVE -0.13492    0.29893  -0.451   0.6523    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.9965 on 173 degrees of freedom
## Multiple R-squared:  0.4059, Adjusted R-squared:  0.3921 
## F-statistic: 29.55 on 4 and 173 DF,  p-value: < 2.2e-16

Describe the LS 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.46 1.29    4.6    4.51 1.19 1.6 7.00  5.40 -0.30
## LS          2 86 4.46 1.42    4.6    4.51 1.41 1.0 7.01  6.01 -0.36
##          kurtosis   se
## BASELINE    -0.46 0.14
## LS          -0.54 0.15
## -------------------------------------------------------- 
## group: 1
##          vars  n mean   sd median trimmed  mad min  max range  skew
## BASELINE    1 92 4.62 1.37   4.90    4.68 1.33 1.6 6.80  5.20 -0.43
## LS          2 92 5.15 1.03   5.39    5.23 0.88 2.2 8.17  5.97 -0.53
##          kurtosis   se
## BASELINE    -0.68 0.14
## LS           0.51 0.11

Create a plot that visualizes LS variable by the GROUP variable

library(ggplot2)
library(influence.ME)

Take a look at the residuals

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

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

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

Comparing Basline to Wave 2 and 3 by Group.

fullModeldata15 <- lme(LS ~ 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
##   509.373 531.6455 -247.6865
## 
## Random effects:
##  Formula: ~1 | ID
##         (Intercept)  Residual
## StdDev:   0.4334834 0.8815661
## 
## Fixed effects: LS ~ GROUP * WAVE + BASELINE 
##                  Value Std.Error DF   t-value p-value
## (Intercept)  2.1195658 0.4151700 87  5.105296  0.0000
## GROUP1       0.7964825 0.4343060 87  1.833920  0.0701
## WAVE        -0.0705970 0.1928513 86 -0.366069  0.7152
## BASELINE     0.5496627 0.0614168 86  8.949712  0.0000
## GROUP1:WAVE -0.1347087 0.2682676 86 -0.502143  0.6169
##  Correlation: 
##             (Intr) GROUP1 WAVE   BASELI
## GROUP1      -0.533                     
## WAVE        -0.697  0.666              
## BASELINE    -0.659 -0.012  0.000       
## GROUP1:WAVE  0.509 -0.926 -0.719 -0.012
## 
## Standardized Within-Group Residuals:
##         Min          Q1         Med          Q3         Max 
## -2.93561677 -0.55407865  0.03719312  0.59524414  2.57195368 
## 
## Number of Observations: 178
## Number of Groups: 89

Another random imputation

data25=MI$imputations[[25]]

library("Zelig")
zelig.fitdata25 <- zelig(LS ~ 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 
## -3.15534 -0.55346  0.04248  0.59014  2.71149 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  1.57101    0.39250   4.003 9.28e-05 ***
## GROUP1       0.77476    0.43822   1.768   0.0788 .  
## WAVE         0.02079    0.19924   0.104   0.9170    
## BASELINE     0.64588    0.05254  12.293  < 2e-16 ***
## GROUP1:WAVE -0.19530    0.27716  -0.705   0.4820    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.9238 on 173 degrees of freedom
## Multiple R-squared:  0.4967, Adjusted R-squared:  0.4851 
## F-statistic: 42.69 on 4 and 173 DF,  p-value: < 2.2e-16

Describe the LS 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.46 1.29   4.60    4.51 1.19 1.6 7.00  5.40 -0.30
## LS          2 86 4.48 1.35   4.65    4.54 1.12 1.0 7.21  6.21 -0.46
##          kurtosis   se
## BASELINE    -0.46 0.14
## LS          -0.13 0.15
## -------------------------------------------------------- 
## group: 1
##          vars  n mean   sd median trimmed  mad  min  max range  skew
## BASELINE    1 92 4.65 1.36   4.90    4.72 1.33 1.60 6.80  5.20 -0.42
## LS          2 92 5.09 1.16   5.25    5.17 0.90 1.81 7.96  6.16 -0.62
##          kurtosis   se
## BASELINE    -0.64 0.14
## LS           0.59 0.12

Create a plot that visualizes LS variable by the GROUP variable

library(ggplot2)
library(influence.ME)

Take a look at the residuals

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

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

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

Comparing Basline to Wave 2 and 3 by Group.

fullModeldata25 <- lme(LS ~ 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
##   484.7606 507.0331 -235.3803
## 
## Random effects:
##  Formula: ~1 | ID
##         (Intercept)  Residual
## StdDev:   0.3038293 0.8586006
## 
## Fixed effects: LS ~ GROUP * WAVE + BASELINE 
##                  Value Std.Error DF   t-value p-value
## (Intercept)  1.5721244 0.3887931 87  4.043602  0.0001
## GROUP1       0.7747803 0.4182601 87  1.852389  0.0674
## WAVE         0.0207851 0.1878274 86  0.110660  0.9121
## BASELINE     0.6456333 0.0553170 86 11.671525  0.0000
## GROUP1:WAVE -0.1952793 0.2612896 86 -0.747367  0.4569
##  Correlation: 
##             (Intr) GROUP1 WAVE   BASELI
## GROUP1      -0.548                     
## WAVE        -0.725  0.674              
## BASELINE    -0.634 -0.012  0.000       
## GROUP1:WAVE  0.530 -0.937 -0.719 -0.015
## 
## Standardized Within-Group Residuals:
##         Min          Q1         Med          Q3         Max 
## -3.26837133 -0.54300737  0.05739964  0.57012011  3.03433081 
## 
## Number of Observations: 178
## Number of Groups: 89

Check assumptions on model without any imputations

Describe the LS 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 kurtosis
## BASELINE    1 86 4.46 1.29    4.6    4.51 1.19 1.6   7   5.4 -0.3    -0.46
## LS          2 59 4.68 1.39    5.0    4.78 0.89 1.0   7   6.0 -0.8     0.19
##            se
## BASELINE 0.14
## LS       0.18
## -------------------------------------------------------- 
## group: 1
##          vars  n mean   sd median trimmed  mad min max range  skew
## BASELINE    1 88 4.63 1.37    4.9    4.69 1.33 1.6 6.8   5.2 -0.41
## LS          2 54 5.33 0.91    5.6    5.45 0.59 2.2 6.6   4.4 -1.48
##          kurtosis   se
## BASELINE    -0.69 0.15
## LS           2.34 0.12

Create a plot that visualizes LS variable by the GROUP variable

library(ggplot2)
library(influence.ME)

Take a look at the residuals

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

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

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

Comparing Basline to Wave 2 and 3 by Group.

fullModel <- lme(LS ~ 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
##   285.2613 304.1007 -135.6306
## 
## Random effects:
##  Formula: ~1 | ID
##         (Intercept) Residual
## StdDev:   0.7555963   0.5452
## 
## Fixed effects: LS ~ GROUP * WAVE + BASELINE 
##                  Value Std.Error DF   t-value p-value
## (Intercept)  1.8976696 0.4445103 66  4.269124  0.0001
## GROUP1       0.8492713 0.3954000 66  2.147879  0.0354
## WAVE         0.0861473 0.1633497 38  0.527379  0.6010
## BASELINE     0.5622375 0.0812697 66  6.918170  0.0000
## GROUP1:WAVE -0.1410534 0.2381139 38 -0.592378  0.5571
##  Correlation: 
##             (Intr) GROUP1 WAVE   BASELI
## GROUP1      -0.337                     
## WAVE        -0.458  0.569              
## BASELINE    -0.799 -0.087 -0.055       
## GROUP1:WAVE  0.320 -0.832 -0.686  0.031
## 
## Standardized Within-Group Residuals:
##          Min           Q1          Med           Q3          Max 
## -3.416279437 -0.327526605  0.002786168  0.285883223  1.985580204 
## 
## Number of Observations: 109
## Number of Groups: 69
Table with P-value

|             |       Value|  Std.Error|  DF|     t-value|    p-value|
|:------------|-----------:|----------:|---:|-----------:|----------:|
|(Intercept)  |   1.8976696|  0.4445103|  66|   4.2691243|  0.0000642|
|GROUP1       |   0.8492713|  0.3954000|  66|   2.1478790|  0.0353990|
|WAVE         |   0.0861473|  0.1633497|  38|   0.5273794|  0.6009960|
|BASELINE     |   0.5622375|  0.0812697|  66|   6.9181697|  0.0000000|
|GROUP1:WAVE  |  -0.1410534|  0.2381139|  38|  -0.5923779|  0.5571063|

Table with confidence intervals

est. lower upper
(Intercept) 1.8976696 1.0307705 2.7645687
GROUP1 0.8492713 0.0781488 1.6203938
WAVE 0.0861473 -0.2368634 0.4091579
BASELINE 0.5622375 0.4037426 0.7207324
GROUP1:WAVE -0.1410534 -0.6119042 0.3297973