Purpose in Life (MLQ-P) with Imputed Data Intention to treat model (ITT)
Loading the dataset
data.test4 <- read.csv("/Volumes/TOSHIBA EXT/Dropbox/ADULT STUDY/adult_study011615.csv")
# Load the psych package
library(psych)
items <- c("MLQ1" ,"MLQ4", "MLQ5", "MLQ6", "MLQ9")
scaleKey <- c(1, 1, 1,1,-1)
data.test4$meanmlq <- scoreItems(scaleKey, items=data.test4[,items], delete=FALSE)$score
library(reshape2); library(car); library(Amelia);library(mitools);library(nlme);library(predictmeans)##
## Attaching package: 'car'
##
## The following object is masked from 'package:psych':
##
## logit
##
## Loading required package: Rcpp
## ##
## ## Amelia II: Multiple Imputation
## ## (Version 1.7.3, built: 2014-11-14)
## ## Copyright (C) 2005-2015 James Honaker, Gary King and Matthew Blackwell
## ## Refer to http://gking.harvard.edu/amelia/ for more information
## ##
## Loading required package: lme4
## Loading required package: Matrix
##
## Attaching package: 'lme4'
##
## The following object is masked from 'package:nlme':
##
## lmListdata <- data.test4[,c("ID", "GROUP", "wave", "meanmlq")]
data <- dcast(data, ID + GROUP ~ wave, mean, value.var = "meanmlq")
data[,3:5] <- apply(data[,3:5],2,function(x) recode(x, "NaN = NA") )Create new data set with ID Group baseline meanmlq and wave so that we have Baseline, time 1 and 2 to compare to
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", "meanmlq", "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"] <- 1Make 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 8 9 10
##
## -- Imputation 4 --
##
## 1 2 3 4 5 6 7
##
## -- Imputation 5 --
##
## 1 2 3 4 5 6 7 8 9 10
##
## -- Imputation 6 --
##
## 1 2 3 4 5 6 7
##
## -- Imputation 7 --
##
## 1 2 3 4 5 6 7 8
##
## -- Imputation 8 --
##
## 1 2 3 4 5 6 7 8 9 10 11 12 13
##
## -- 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 9
##
## -- Imputation 12 --
##
## 1 2 3 4 5 6 7
##
## -- Imputation 13 --
##
## 1 2 3 4 5 6 7 8 9 10
##
## -- 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
##
## -- Imputation 18 --
##
## 1 2 3 4 5 6 7
##
## -- Imputation 19 --
##
## 1 2 3 4 5 6 7 8 9
##
## -- Imputation 20 --
##
## 1 2 3 4 5 6 7
##
## -- Imputation 21 --
##
## 1 2 3 4 5 6 7 8 9 10 11 12 13
##
## -- Imputation 22 --
##
## 1 2 3 4 5 6 7 8 9 10 11
##
## -- Imputation 23 --
##
## 1 2 3 4 5 6 7 8
##
## -- Imputation 24 --
##
## 1 2 3 4
##
## -- Imputation 25 --
##
## 1 2 3 4 5 6 7 8 9
##
## -- Imputation 26 --
##
## 1 2 3 4 5 6 7 8
##
## -- Imputation 27 --
##
## 1 2 3 4 5 6 7
##
## -- Imputation 28 --
##
## 1 2 3 4 5 6 7
##
## -- 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
##
## -- Imputation 32 --
##
## 1 2 3 4 5 6 7
##
## -- Imputation 33 --
##
## 1 2 3 4 5 6 7 8
##
## -- Imputation 34 --
##
## 1 2 3 4 5 6 7 8 9
##
## -- Imputation 35 --
##
## 1 2 3 4 5 6 7 8 9 10
##
## -- Imputation 36 --
##
## 1 2 3 4 5 6 7
##
## -- Imputation 37 --
##
## 1 2 3 4 5 6 7 8 9 10
##
## -- Imputation 38 --
##
## 1 2 3 4 5 6 7
##
## -- Imputation 39 --
##
## 1 2 3 4 5 6 7
##
## -- Imputation 40 --
##
## 1 2 3 4 5 6 7
##
## -- Imputation 41 --
##
## 1 2 3 4 5
##
## -- Imputation 42 --
##
## 1 2 3 4 5 6
##
## -- 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
##
## -- Imputation 45 --
##
## 1 2 3 4 5 6 7
##
## -- Imputation 46 --
##
## 1 2 3 4 5 6 7 8 9 10 11
##
## -- Imputation 47 --
##
## 1 2 3 4 5 6 7 8 9
##
## -- Imputation 48 --
##
## 1 2 3 4 5 6 7 8
##
## -- Imputation 49 --
##
## 1 2 3 4 5
##
## -- Imputation 50 --
##
## 1 2 3 4 5 6 7 8 9 10 11Creating new dataset with missing data imputed
data(MI$imputations)## Warning in data(MI$imputations): data set 'MI$imputations' not foundallimplogreg<-lapply(MI$imputations,function(X) {lme(meanmlq ~ 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.72370499 0.36738733 1.00175663 2.4456534 33 %
## GROUP1 0.78712436 0.37345276 0.05344626 1.5208024 31 %
## WAVE 0.10143086 0.16828571 -0.22928928 0.4321510 33 %
## BASELINE 0.65552817 0.05873157 0.53978496 0.7712714 47 %
## GROUP1:WAVE -0.02969385 0.23520186 -0.49196776 0.4325801 34 %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':
##
## citezelig.fit <- zelig(meanmlq ~ 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.71638966 0.37376699 4.5921382 5.531921e-06
## GROUP1 0.78744990 0.40103417 1.9635481 4.998745e-02
## WAVE 0.10143086 0.18347232 0.5528401 5.805665e-01
## BASELINE 0.65706711 0.05629104 11.6726773 3.497378e-24
## GROUP1:WAVE -0.02969923 0.25624800 -0.1159003 9.077695e-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).Describe the meanmlq 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.75 1.27 4.8 4.82 1.19 1.4 7 5.6 -0.43
## meanmlq 2 59 5.17 1.15 5.4 5.22 1.19 2.4 7 4.6 -0.45
## kurtosis se
## BASELINE -0.20 0.14
## meanmlq -0.38 0.15
## --------------------------------------------------------
## group: 1
## vars n mean sd median trimmed mad min max range skew
## BASELINE 1 88 4.50 1.47 4.4 4.50 1.63 1.8 7 5.2 0.05
## meanmlq 2 54 5.76 1.03 6.0 5.84 0.89 3.8 7 3.2 -0.61
## kurtosis se
## BASELINE -1.05 0.16
## meanmlq -0.82 0.14Create a plot that visualizes meanmlq 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':
##
## influenceTake a look at the residuals
residual <- lm(meanmlq ~ BASELINE, data=data2)$residualPlot 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$meanmlq))
sel2 <- which(!is.na(data2$BASELINE))
data2$residual[intersect(sel1,sel2)] <- residual
qplot(GROUP, meanmlq, 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 meanmlq and the Residuals
# Load the nlme package
library(nlme)
with(data2, boxplot(meanmlq ~ WAVE + GROUP))
with(data2, boxplot(residual ~ WAVE + GROUP))
Comparing Basline to Wave 2 and 3 by Group.
fullModel <- lme(meanmlq ~ 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")
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
## 240.2404 259.0798 -113.1202
##
## Random effects:
## Formula: ~1 | ID
## (Intercept) Residual
## StdDev: 0.5321759 0.4951252
##
## Fixed effects: meanmlq ~ GROUP * WAVE + BASELINE
## Value Std.Error DF t-value p-value
## (Intercept) 1.6901896 0.3918256 66 4.313628 0.0001
## GROUP1 0.8795864 0.3380075 66 2.602269 0.0114
## WAVE 0.0773449 0.1454611 38 0.531722 0.5980
## BASELINE 0.6589066 0.0654276 66 10.070773 0.0000
## GROUP1:WAVE -0.0842401 0.2125839 38 -0.396268 0.6941
## Correlation:
## (Intr) GROUP1 WAVE BASELI
## GROUP1 -0.426
## WAVE -0.484 0.590
## BASELINE -0.811 0.037 -0.032
## GROUP1:WAVE 0.333 -0.869 -0.684 0.021
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -2.4086341 -0.4440518 0.1101479 0.3570213 2.0486129
##
## Number of Observations: 109
## Number of Groups: 69Table with P-value
| | Value| Std.Error| DF| t-value| p-value|
|:------------|-----------:|----------:|---:|-----------:|----------:|
|(Intercept) | 1.6901896| 0.3918256| 66| 4.3136276| 0.0000549|
|GROUP1 | 0.8795864| 0.3380075| 66| 2.6022691| 0.0114235|
|WAVE | 0.0773449| 0.1454611| 38| 0.5317219| 0.5980136|
|BASELINE | 0.6589066| 0.0654276| 66| 10.0707734| 0.0000000|
|GROUP1:WAVE | -0.0842401| 0.2125839| 38| -0.3962676| 0.6941241|
Table with confidence intervals
| est. | lower | upper | |
|---|---|---|---|
| (Intercept) | 1.6901896 | 0.9260380 | 2.4543411 |
| GROUP1 | 0.8795864 | 0.2203927 | 1.5387801 |
| WAVE | 0.0773449 | -0.2102926 | 0.3649823 |
| BASELINE | 0.6589066 | 0.5313075 | 0.7865058 |
| GROUP1:WAVE | -0.0842401 | -0.5046074 | 0.3361272 |