Loading the dataset
setwd("~/Dropbox/ADULT STUDY")
data.test4 <- read.csv("adult_study011615.csv")
# Load the psych package
library(psych)
items <- c("GRIT8", "GRIT9", "GRIT10", "GRIT11")
scaleKey <- c(1, 1, 1, 1)
data.test4$meanGRIT <- 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':
##
## lmList
#Remove the meanGRIT and ID Group and wave from dtat.test4 and create a new #dataset with only those variables.
data <- data.test4[,c("ID", "GROUP", "wave", "meanGRIT")]
#Use dcast to cnage from long-format data to wide format data
data <- dcast(data, ID + GROUP ~ wave, mean, value.var = "meanGRIT")
# 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 meanGRIT and wave 2 and 3 of meanGRIT. 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 compaGRIT Time 2 to Baseline) “wave 1” and then the second 89 we call “wave 2” which compaGRIT 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", "meanGRIT", "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). This line of data makes Group 2 become Group 1 so that Group 2 which were the people who dropped out become Group 1 i.e. part of the treatment group.
data2[which(data2$GROUP ==2), "GROUP"] <- 1
data2$GROUP <-as.factor(data2$GROUP)
data2$ID <-as.factor(data2$ID)
Imputing missing data. 50 datasets are created.
MI <- amelia(data2, 50, idvars = c("ID"), ords = "GROUP")
## -- Imputation 1 --
##
## 1 2 3 4 5 6 7 8 9
##
## -- 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 7 8
##
## -- Imputation 5 --
##
## 1 2 3 4 5 6 7
##
## -- Imputation 6 --
##
## 1 2 3 4 5 6 7 8 9 10
##
## -- Imputation 7 --
##
## 1 2 3 4 5 6 7 8
##
## -- 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 9
##
## -- Imputation 11 --
##
## 1 2 3 4 5 6 7 8 9
##
## -- Imputation 12 --
##
## 1 2 3 4 5 6 7 8 9 10 11 12
##
## -- Imputation 13 --
##
## 1 2 3 4 5 6 7 8 9
##
## -- Imputation 14 --
##
## 1 2 3 4 5 6 7 8
##
## -- Imputation 15 --
##
## 1 2 3 4 5 6 7 8
##
## -- Imputation 16 --
##
## 1 2 3 4 5 6 7
##
## -- Imputation 17 --
##
## 1 2 3 4 5 6 7 8
##
## -- Imputation 18 --
##
## 1 2 3 4 5 6 7 8
##
## -- Imputation 19 --
##
## 1 2 3 4 5 6
##
## -- Imputation 20 --
##
## 1 2 3 4 5 6
##
## -- Imputation 21 --
##
## 1 2 3 4 5 6 7 8 9 10 11
##
## -- Imputation 22 --
##
## 1 2 3 4 5 6 7 8
##
## -- Imputation 23 --
##
## 1 2 3 4 5 6 7 8 9 10
##
## -- Imputation 24 --
##
## 1 2 3 4 5 6
##
## -- Imputation 25 --
##
## 1 2 3 4 5 6 7
##
## -- Imputation 26 --
##
## 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15
##
## -- Imputation 27 --
##
## 1 2 3 4 5 6 7 8 9
##
## -- Imputation 28 --
##
## 1 2 3 4 5 6 7 8 9 10
##
## -- Imputation 29 --
##
## 1 2 3 4 5 6 7
##
## -- Imputation 30 --
##
## 1 2 3 4 5 6 7
##
## -- Imputation 31 --
##
## 1 2 3 4 5 6 7 8
##
## -- Imputation 32 --
##
## 1 2 3 4 5 6 7
##
## -- Imputation 33 --
##
## 1 2 3 4 5 6 7 8 9
##
## -- Imputation 34 --
##
## 1 2 3 4 5 6 7 8 9
##
## -- Imputation 35 --
##
## 1 2 3 4 5 6 7 8 9 10
##
## -- Imputation 36 --
##
## 1 2 3 4 5 6 7 8
##
## -- Imputation 37 --
##
## 1 2 3 4 5 6 7 8 9
##
## -- Imputation 38 --
##
## 1 2 3 4 5 6 7 8 9
##
## -- Imputation 39 --
##
## 1 2 3 4 5 6
##
## -- Imputation 40 --
##
## 1 2 3 4 5 6
##
## -- Imputation 41 --
##
## 1 2 3 4 5 6 7 8 9 10
##
## -- Imputation 42 --
##
## 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16
##
## -- 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
##
## -- Imputation 46 --
##
## 1 2 3 4 5 6 7
##
## -- Imputation 47 --
##
## 1 2 3 4 5 6
##
## -- Imputation 48 --
##
## 1 2 3 4 5 6 7 8 9
##
## -- Imputation 49 --
##
## 1 2 3 4 5 6 7 8 9
##
## -- Imputation 50 --
##
## 1 2 3 4 5 6 7 8
Creating new dataset with missing data imputed. On the second line of code a repeated measure analysis is condicted on the data set which has the data imputed.
data(MI$imputations)
## Warning in data(MI$imputations): data set 'MI$imputations' not found
allimplogreg<-lapply(MI$imputations,function(X) {lme(meanGRIT ~ 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<-summary(MIcombine(betas,vars))
## Multiple imputation results:
## MIcombine.default(betas, vars)
## results se (lower upper) missInfo
## (Intercept) 1.76409494 0.29812789 1.1768711 2.3513188 45 %
## GROUP1 0.02044326 0.22281974 -0.4171392 0.4580257 29 %
## WAVE -0.05103264 0.10877538 -0.2650811 0.1630158 41 %
## BASELINE 0.58089510 0.06459203 0.4534286 0.7083616 53 %
## GROUP1:WAVE 0.07062662 0.14706391 -0.2185730 0.3598262 37 %
summary
## results se (lower upper) missInfo
## (Intercept) 1.76409494 0.29812789 1.1768711 2.3513188 45 %
## GROUP1 0.02044326 0.22281974 -0.4171392 0.4580257 29 %
## WAVE -0.05103264 0.10877538 -0.2650811 0.1630158 41 %
## BASELINE 0.58089510 0.06459203 0.4534286 0.7083616 53 %
## GROUP1:WAVE 0.07062662 0.14706391 -0.2185730 0.3598262 37 %
library(pander)
Table
##
## ----------------------------------------------------------------
## results se (lower upper) missInfo
## ----------------- --------- ------- -------- -------- ----------
## **(Intercept)** 1.764 0.2981 1.177 2.351 45 %
##
## **GROUP1** 0.02044 0.2228 -0.4171 0.458 29 %
##
## **WAVE** -0.05103 0.1088 -0.2651 0.163 41 %
##
## **BASELINE** 0.5809 0.06459 0.4534 0.7084 53 %
##
## **GROUP1:WAVE** 0.07063 0.1471 -0.2186 0.3598 37 %
## ----------------------------------------------------------------
Check GRITults 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(meanGRIT ~ GROUP * WAVE + BASELINE, random = ~1 | ID, data = MI$imputations, model = "ls", digits = 4, cite = F)
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.75896241 0.29817229 5.89914782 1.214795e-08
## GROUP1 0.02082638 0.23598218 0.08825406 9.296977e-01
## WAVE -0.05103264 0.11546262 -0.44198406 6.587484e-01
## BASELINE 0.58219716 0.06291863 9.25317599 1.292924e-16
## GROUP1:WAVE 0.07059381 0.15662775 0.45071075 6.524064e-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).
summary1<-summary(zelig.fit)
Table with p-values
##
## ----------------------------------------------------------
## Value Std. Error t-stat p-value
## ----------------- -------- ------------ -------- ---------
## **(Intercept)** 1.759 0.2982 5.899 1.215e-08
##
## **GROUP1** 0.02083 0.236 0.08825 0.9297
##
## **WAVE** -0.05103 0.1155 -0.442 0.6587
##
## **BASELINE** 0.5822 0.06292 9.253 1.293e-16
##
## **GROUP1:WAVE** 0.07059 0.1566 0.4507 0.6524
## ----------------------------------------------------------
Check assumptions with Random Computations. Zailig fit with just one of the imputed data sets.
data1=MI$imputations[[1]]
library("Zelig")
zelig.fitdata1 <- zelig(meanGRIT ~ 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.1321 -0.2304 -0.0114 0.2777 1.4154
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.58608 0.20510 7.733 8.26e-13 ***
## GROUP1 -0.15156 0.18857 -0.804 0.423
## WAVE -0.03061 0.08564 -0.357 0.721
## BASELINE 0.61909 0.03908 15.842 < 2e-16 ***
## GROUP1:WAVE 0.19247 0.11913 1.616 0.108
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.3971 on 173 degrees of freedom
## Multiple R-squared: 0.5953, Adjusted R-squared: 0.586
## F-statistic: 63.63 on 4 and 173 DF, p-value: < 2.2e-16
Describe the meanGRIT 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 3.94 0.66 4 3.96 0.74 2.50 5.00 2.50 -0.28
## meanGRIT 2 86 3.98 0.62 4 3.99 0.69 2.25 5.47 3.22 -0.29
## kurtosis se
## BASELINE -0.71 0.07
## meanGRIT 0.06 0.07
## --------------------------------------------------------
## group: 1
## vars n mean sd median trimmed mad min max range skew
## BASELINE 1 92 3.67 0.85 3.75 3.76 0.74 1.0 5.00 4.00 -0.86
## meanGRIT 2 92 3.95 0.61 4.00 4.00 0.61 1.5 5.03 3.53 -0.99
## kurtosis se
## BASELINE 0.67 0.09
## meanGRIT 1.92 0.06
Create a plot that visualizes meanGRIT variable by the GROUP variable. Load the proper packedes.
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. Of a random selected dataset with imputed data.
residual <- lm(meanGRIT ~ BASELINE, data=data1)$residual
Plot the residuals to see that they are random
# A density plot
plot(density(residual))
# A quantile normal plot to checking normality
qqnorm(residual)
qqline(residual)
Checking the different between intervention and control groups residuals within the selected imputed dataset. This allows us to control for individual unsystematic differences.
data2$residual <- NA
sel1 <- which(!is.na(data1$meanGRIT))
sel2 <- which(!is.na(data1$BASELINE))
data1$residual[intersect(sel1,sel2)] <- residual
qplot(GROUP, meanGRIT, data=data1, geom="boxplot")
Plot of the difference between intervention and control groups within the selected imputed dataset.
qplot(GROUP, residual, data=data1, geom="boxplot")
Two way repeated measuGRIT on dataset Randomly Selected Imputed Data ======================================================== Graphing the Two-Way Interaction. Both meanGRIT and the residuals
# nlme package
with(data1, boxplot(meanGRIT ~ WAVE + GROUP))
with(data1, boxplot(residual ~ WAVE + GROUP))
Comparing Basline to Wave 2 and 3 by Group.
fullModeldata1 <- lme(meanGRIT ~ GROUP * WAVE + BASELINE, random = ~1 | ID, data = data1, method = "ML", na.action = "na.omit")
CookD(fullModeldata1)
plot(fullModeldata1, 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(fullModeldata1)
## Linear mixed-effects model fit by maximum likelihood
## Data: data1
## AIC BIC logLik
## 180.9777 203.2502 -83.48887
##
## Random effects:
## Formula: ~1 | ID
## (Intercept) Residual
## StdDev: 0.1824228 0.3463772
##
## Fixed effects: meanGRIT ~ GROUP * WAVE + BASELINE
## Value Std.Error DF t-value p-value
## (Intercept) 1.5897450 0.20973565 87 7.579756 0.0000
## GROUP1 -0.1517799 0.17150648 87 -0.884981 0.3786
## WAVE -0.0306054 0.07577343 86 -0.403907 0.6873
## BASELINE 0.6181644 0.04308111 86 14.348851 0.0000
## GROUP1:WAVE 0.1924443 0.10540287 86 1.825797 0.0714
## Correlation:
## (Intr) GROUP1 WAVE BASELI
## GROUP1 -0.469
## WAVE -0.542 0.663
## BASELINE -0.810 0.059 0.000
## GROUP1:WAVE 0.382 -0.921 -0.719 0.009
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -2.21925514 -0.59553981 -0.01969639 0.63980059 2.83755615
##
## Number of Observations: 178
## Number of Groups: 89
Another random selected imputation
data10=MI$imputations[[10]]
library("Zelig")
zelig.fitdata10 <- zelig(meanGRIT ~ 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.0995 -0.3062 -0.0277 0.2784 1.3040
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.69284 0.22969 7.370 6.73e-12 ***
## GROUP1 -0.02067 0.21132 -0.098 0.922
## WAVE -0.05563 0.09600 -0.580 0.563
## BASELINE 0.58930 0.04373 13.475 < 2e-16 ***
## GROUP1:WAVE 0.15086 0.13355 1.130 0.260
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.4451 on 173 degrees of freedom
## Multiple R-squared: 0.5141, Adjusted R-squared: 0.5029
## F-statistic: 45.76 on 4 and 173 DF, p-value: < 2.2e-16
Describe the meanGRIT 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 3.94 0.66 4 3.96 0.74 2.50 5.00 2.50 -0.28
## meanGRIT 2 86 3.93 0.63 4 3.95 0.74 2.25 5.02 2.77 -0.29
## kurtosis se
## BASELINE -0.71 0.07
## meanGRIT -0.33 0.07
## --------------------------------------------------------
## group: 1
## vars n mean sd median trimmed mad min max range skew
## BASELINE 1 92 3.68 0.85 3.75 3.76 0.74 1.0 5.00 4.00 -0.86
## meanGRIT 2 92 3.98 0.64 4.00 4.01 0.66 1.5 5.78 4.28 -0.56
## kurtosis se
## BASELINE 0.65 0.09
## meanGRIT 1.94 0.07
Create a plot that visualizes meanGRIT variable by the GROUP variable
library(ggplot2)
library(influence.ME)
Take a look at the residuals
residual <- lm(meanGRIT ~ 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.
data10$residual <- NA
sel1 <- which(!is.na(data10$meanGRIT))
sel2 <- which(!is.na(data10$BASELINE))
data10$residual[intersect(sel1,sel2)] <- residual
qplot(GROUP, meanGRIT, data=data10, geom="boxplot")
Plot of the difference between intervention and control groups.
qplot(GROUP, residual, data=data10, geom="boxplot")
Two way repeated measuGRIT ======================================================== Graphing the Two-Way Interaction. Both meanGRIT and the residuals
# Load the nlme package
library(nlme)
with(data10, boxplot(meanGRIT ~ WAVE + GROUP))
with(data10, boxplot(residual ~ WAVE + GROUP))
Comparing Basline to Wave 2 and 3 by Group.
fullModeldata10 <- lme(meanGRIT ~ GROUP * WAVE + BASELINE, random = ~1 | ID, data = data10, method = "ML", na.action = "na.omit")
CookD(fullModeldata10)
plot(fullModeldata10, 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(fullModeldata10)
## Linear mixed-effects model fit by maximum likelihood
## Data: data10
## AIC BIC logLik
## 224.0421 246.3146 -105.0211
##
## Random effects:
## Formula: ~1 | ID
## (Intercept) Residual
## StdDev: 0.1671221 0.4057728
##
## Fixed effects: meanGRIT ~ GROUP * WAVE + BASELINE
## Value Std.Error DF t-value p-value
## (Intercept) 1.6999015 0.23301534 87 7.295234 0.0000
## GROUP1 -0.0210392 0.19873705 87 -0.105865 0.9159
## WAVE -0.0556325 0.08876683 86 -0.626726 0.5325
## BASELINE 0.5875116 0.04672874 86 12.572809 0.0000
## GROUP1:WAVE 0.1507823 0.12348626 86 1.221045 0.2254
## Correlation:
## (Intr) GROUP1 WAVE BASELI
## GROUP1 -0.478
## WAVE -0.571 0.670
## BASELINE -0.790 0.048 0.000
## GROUP1:WAVE 0.399 -0.931 -0.719 0.015
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -2.24224335 -0.62544958 -0.05411535 0.60854423 2.68605503
##
## Number of Observations: 178
## Number of Groups: 89
Another random selected imputation
data15=MI$imputations[[15]]
library("Zelig")
zelig.fitdata15 <- zelig(meanGRIT ~ 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.22808 -0.24351 0.00976 0.28835 1.16042
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.10044 0.22465 9.350 <2e-16 ***
## GROUP1 -0.15311 0.20752 -0.738 0.462
## WAVE -0.08642 0.09414 -0.918 0.360
## BASELINE 0.51682 0.04268 12.108 <2e-16 ***
## GROUP1:WAVE 0.19502 0.13095 1.489 0.138
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.4365 on 173 degrees of freedom
## Multiple R-squared: 0.4632, Adjusted R-squared: 0.4508
## F-statistic: 37.32 on 4 and 173 DF, p-value: < 2.2e-16
Describe the meanGRIT 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 3.94 0.66 4 3.96 0.74 2.50 5.00 2.50 -0.28
## meanGRIT 2 86 4.01 0.60 4 4.02 0.65 2.25 5.23 2.98 -0.39
## kurtosis se
## BASELINE -0.71 0.07
## meanGRIT 0.15 0.06
## --------------------------------------------------------
## group: 1
## vars n mean sd median trimmed mad min max range skew
## BASELINE 1 92 3.65 0.86 3.75 3.73 0.74 1.0 5.00 4.00 -0.79
## meanGRIT 2 92 4.00 0.59 4.00 4.05 0.44 1.5 5.24 3.74 -1.15
## kurtosis se
## BASELINE 0.47 0.09
## meanGRIT 3.06 0.06
Create a plot that visualizes meanGRIT variable by the GROUP variable
library(ggplot2)
library(influence.ME)
Take a look at the residuals
residual <- lm(meanGRIT ~ 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$meanGRIT))
sel2 <- which(!is.na(data15$BASELINE))
data15$residual[intersect(sel1,sel2)] <- residual
qplot(GROUP, meanGRIT, data=data15, geom="boxplot")
Plot of the difference between intervention and control groups.
qplot(GROUP, residual, data=data15, geom="boxplot")
Two way repeated measuGRIT ======================================================== Graphing the Two-Way Interaction. Both meanGRIT and the residuals
# Load the nlme package
library(nlme)
with(data15, boxplot(meanGRIT ~ WAVE + GROUP))
with(data15, boxplot(residual ~ WAVE + GROUP))
Comparing Basline to Wave 2 and 3 by Group.
fullModeldata15 <- lme(meanGRIT ~ GROUP * WAVE + BASELINE, random = ~1 | ID, data = data15, method = "ML", na.action = "na.omit")
CookD(fullModeldata15)
plot(fullModeldata15, 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(fullModeldata15)
## Linear mixed-effects model fit by maximum likelihood
## Data: data15
## AIC BIC logLik
## 218.5184 240.7908 -102.2592
##
## Random effects:
## Formula: ~1 | ID
## (Intercept) Residual
## StdDev: 0.1145134 0.4148188
##
## Fixed effects: meanGRIT ~ GROUP * WAVE + BASELINE
## Value Std.Error DF t-value p-value
## (Intercept) 2.1028731 0.22624173 87 9.294807 0.0000
## GROUP1 -0.1533159 0.20162285 87 -0.760409 0.4491
## WAVE -0.0864175 0.09074571 86 -0.952304 0.3436
## BASELINE 0.5162085 0.04414785 86 11.692720 0.0000
## GROUP1:WAVE 0.1950332 0.12622990 86 1.545063 0.1260
## Correlation:
## (Intr) GROUP1 WAVE BASELI
## GROUP1 -0.514
## WAVE -0.602 0.675
## BASELINE -0.769 0.072 0.000
## GROUP1:WAVE 0.440 -0.940 -0.719 -0.010
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -2.73305358 -0.57448212 0.04902422 0.64147796 2.65659916
##
## Number of Observations: 178
## Number of Groups: 89
Another randomly selected imputation
data25=MI$imputations[[25]]
library("Zelig")
zelig.fitdata25 <- zelig(meanGRIT ~ 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.22104 -0.29279 0.00718 0.31784 1.19863
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.96973 0.21905 8.992 4.16e-16 ***
## GROUP1 -0.06262 0.20147 -0.311 0.756
## WAVE -0.05836 0.09147 -0.638 0.524
## BASELINE 0.53934 0.04173 12.923 < 2e-16 ***
## GROUP1:WAVE 0.12393 0.12724 0.974 0.331
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.4241 on 173 degrees of freedom
## Multiple R-squared: 0.4927, Adjusted R-squared: 0.481
## F-statistic: 42.01 on 4 and 173 DF, p-value: < 2.2e-16
Describe the meanGRIT 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 3.94 0.66 4 3.96 0.74 2.50 5 2.50 -0.28
## meanGRIT 2 86 4.01 0.57 4 4.03 0.54 2.25 5 2.75 -0.47
## kurtosis se
## BASELINE -0.71 0.07
## meanGRIT 0.35 0.06
## --------------------------------------------------------
## group: 1
## vars n mean sd median trimmed mad min max range skew
## BASELINE 1 92 3.69 0.85 3.75 3.77 0.74 1.0 5.00 4.00 -0.90
## meanGRIT 2 92 3.99 0.60 4.00 4.03 0.60 1.5 5.11 3.61 -0.82
## kurtosis se
## BASELINE 0.72 0.09
## meanGRIT 1.79 0.06
Create a plot that visualizes meanGRIT variable by the GROUP variable
library(ggplot2)
library(influence.ME)
Take a look at the residuals
residual <- lm(meanGRIT ~ 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.
data25$residual <- NA
sel1 <- which(!is.na(data25$meanGRIT))
sel2 <- which(!is.na(data25$BASELINE))
data25$residual[intersect(sel1,sel2)] <- residual
qplot(GROUP, meanGRIT, data=data25, geom="boxplot")
Plot of the difference between intervention and control groups.
qplot(GROUP, residual, data=data25, geom="boxplot")
Two way repeated measuGRIT ======================================================== Graphing the Two-Way Interaction. Both meanGRIT and the residuals
# Load the nlme package
library(nlme)
with(data25, boxplot(meanGRIT ~ WAVE + GROUP))
with(data25, boxplot(residual ~ WAVE + GROUP))
Comparing Basline to Wave 2 and 3 by Group.
fullModeldata25 <- lme(meanGRIT ~ GROUP * WAVE + BASELINE, random = ~1 | ID, data = data25, method = "ML", na.action = "na.omit")
CookD(fullModeldata25)
plot(fullModeldata25, 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(fullModeldata25)
## Linear mixed-effects model fit by maximum likelihood
## Data: data25
## AIC BIC logLik
## 208.1828 230.4553 -97.09139
##
## Random effects:
## Formula: ~1 | ID
## (Intercept) Residual
## StdDev: 0.1173026 0.4013539
##
## Fixed effects: meanGRIT ~ GROUP * WAVE + BASELINE
## Value Std.Error DF t-value p-value
## (Intercept) 1.9720326 0.22084070 87 8.929661 0.0000
## GROUP1 -0.0627678 0.19506607 87 -0.321777 0.7484
## WAVE -0.0583623 0.08780014 86 -0.664717 0.5080
## BASELINE 0.5387598 0.04332724 86 12.434667 0.0000
## GROUP1:WAVE 0.1239315 0.12212692 86 1.014776 0.3131
## Correlation:
## (Intr) GROUP1 WAVE BASELI
## GROUP1 -0.500
## WAVE -0.596 0.675
## BASELINE -0.773 0.058 0.000
## GROUP1:WAVE 0.429 -0.939 -0.719 -0.001
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -2.86149102 -0.68802592 0.01113375 0.72529090 2.60597303
##
## Number of Observations: 178
## Number of Groups: 89
Check assumptions on model without any imputations
Describe the meanGRIT 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 3.94 0.66 4 3.96 0.74 2.50 5 2.50 -0.28
## meanGRIT 2 59 4.04 0.62 4 4.08 0.74 2.25 5 2.75 -0.60
## kurtosis se
## BASELINE -0.71 0.07
## meanGRIT 0.32 0.08
## --------------------------------------------------------
## group: 1
## vars n mean sd median trimmed mad min max range skew
## BASELINE 1 88 3.68 0.87 3.75 3.76 0.74 1.0 5.00 4.00 -0.87
## meanGRIT 2 54 4.00 0.60 4.00 4.05 0.56 1.5 4.75 3.25 -1.34
## kurtosis se
## BASELINE 0.59 0.09
## meanGRIT 3.66 0.08
Create a plot that visualizes meanGRIT variable by the GROUP variable
library(ggplot2)
library(influence.ME)
Take a look at the residuals
residual <- lm(meanGRIT ~ 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.
data$residual <- NA
sel1 <- which(!is.na(data2$meanGRIT))
sel2 <- which(!is.na(data2$BASELINE))
data2$residual[intersect(sel1,sel2)] <- residual
qplot(GROUP, meanGRIT, 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 measuGRIT ======================================================== Graphing the Two-Way Interaction. Both meanGRIT and the residuals
# Load the nlme package
library(nlme)
with(data2, boxplot(meanGRIT ~ WAVE + GROUP))
with(data2, boxplot(residual ~ WAVE + GROUP))
Comparing Basline to Wave 2 and 3 by Group.
fullModel <- lme(meanGRIT ~ GROUP * WAVE + BASELINE, random = ~1 | ID, data = data2, method = "ML", na.action = "na.omit")
CookD(fullModel)
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
## 136.4122 155.2516 -61.20609
##
## Random effects:
## Formula: ~1 | ID
## (Intercept) Residual
## StdDev: 0.2483191 0.3550153
##
## Fixed effects: meanGRIT ~ GROUP * WAVE + BASELINE
## Value Std.Error DF t-value p-value
## (Intercept) 1.6179408 0.28660152 66 5.645262 0.0000
## GROUP1 0.0346653 0.22791092 66 0.152100 0.8796
## WAVE -0.0354711 0.10128172 38 -0.350222 0.7281
## BASELINE 0.6098979 0.06251865 66 9.755455 0.0000
## GROUP1:WAVE 0.0682428 0.14816740 38 0.460579 0.6477
## Correlation:
## (Intr) GROUP1 WAVE BASELI
## GROUP1 -0.432
## WAVE -0.431 0.608
## BASELINE -0.844 0.082 -0.067
## GROUP1:WAVE 0.331 -0.906 -0.681 0.002
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -2.1510859 -0.5247612 0.1533065 0.5991881 2.1739981
##
## Number of Observations: 109
## Number of Groups: 69
Table with confidence intervals
Table with confidence intervals
##
## ------------------------------------------
## lower est. upper
## ----------------- ------- -------- -------
## **(Intercept)** 1.059 1.618 2.177
##
## **GROUP1** -0.4098 0.03467 0.4791
##
## **WAVE** -0.2357 -0.03547 0.1648
##
## **BASELINE** 0.488 0.6099 0.7318
##
## **GROUP1:WAVE** -0.2247 0.06824 0.3612
## ------------------------------------------
Table with p-values
##
## ---------------------------------------------------------------
## Value Std.Error DF t-value p-value
## ----------------- -------- ----------- ---- --------- ---------
## **(Intercept)** 1.618 0.2866 66 5.645 3.762e-07
##
## **GROUP1** 0.03467 0.2279 66 0.1521 0.8796
##
## **WAVE** -0.03547 0.1013 38 -0.3502 0.7281
##
## **BASELINE** 0.6099 0.06252 66 9.755 2.02e-14
##
## **GROUP1:WAVE** 0.06824 0.1482 38 0.4606 0.6477
## ---------------------------------------------------------------