title: “Happiness (Subjective Happiness (Lyubomirsky & Lepper, 1999, only first three questions) with Imputed Data Intention to treat model (ITT)”

Loading the dataset

setwd("~/Dropbox/ADULT STUDY")
data.test4 <- read.csv("adult_study011615.csv")
# Load the psych package
library(psych)
items <- c("HAPPI1" ,"HAPPI2", "HAPPI3")
scaleKey <- c(1, 1, 1)
data.test4$meanHAPPI  <- 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 meanHAPPI and ID Group and wave from dtat.test4 and create a new #dataset with only those variables.
data <- data.test4[,c("ID", "GROUP", "wave", "meanHAPPI")]
#Use dcast to cnage from long-format data to wide format data
data <- dcast(data, ID + GROUP ~ wave, mean, value.var = "meanHAPPI")
# 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 meanHAPPI and wave 2 and 3 of meanHAPPI. 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", "meanHAPPI", "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

Make GROUP and ID a factor

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
## 
## -- Imputation 2 --
## 
##   1  2  3  4  5  6  7  8
## 
## -- Imputation 3 --
## 
##   1  2  3  4  5  6  7  8
## 
## -- Imputation 4 --
## 
##   1  2  3  4  5  6  7
## 
## -- Imputation 5 --
## 
##   1  2  3  4  5  6  7  8
## 
## -- Imputation 6 --
## 
##   1  2  3  4  5  6  7  8  9 10
## 
## -- Imputation 7 --
## 
##   1  2  3  4  5  6  7  8  9 10
## 
## -- Imputation 8 --
## 
##   1  2  3  4  5  6
## 
## -- 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
## 
## -- Imputation 12 --
## 
##   1  2  3  4  5  6  7  8
## 
## -- Imputation 13 --
## 
##   1  2  3  4  5  6
## 
## -- Imputation 14 --
## 
##   1  2  3  4  5  6  7  8
## 
## -- Imputation 15 --
## 
##   1  2  3  4  5  6  7  8  9 10
## 
## -- Imputation 16 --
## 
##   1  2  3  4  5  6  7  8  9 10 11 12
## 
## -- Imputation 17 --
## 
##   1  2  3  4  5  6  7  8  9 10 11 12 13 14
## 
## -- Imputation 18 --
## 
##   1  2  3  4  5  6  7  8
## 
## -- Imputation 19 --
## 
##   1  2  3  4  5  6  7  8  9 10 11
## 
## -- Imputation 20 --
## 
##   1  2  3  4  5  6  7  8  9 10
## 
## -- Imputation 21 --
## 
##   1  2  3  4  5  6  7  8  9 10
## 
## -- Imputation 22 --
## 
##   1  2  3  4  5
## 
## -- Imputation 23 --
## 
##   1  2  3  4  5  6  7
## 
## -- 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
## 
## -- Imputation 29 --
## 
##   1  2  3  4  5  6  7  8
## 
## -- Imputation 30 --
## 
##   1  2  3  4  5  6  7  8  9
## 
## -- Imputation 31 --
## 
##   1  2  3  4  5  6  7  8
## 
## -- Imputation 32 --
## 
##   1  2  3  4  5  6  7  8  9
## 
## -- Imputation 33 --
## 
##   1  2  3  4  5  6  7  8
## 
## -- Imputation 34 --
## 
##   1  2  3  4  5  6  7  8
## 
## -- 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 10 11
## 
## -- 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  6  7  8
## 
## -- Imputation 42 --
## 
##   1  2  3  4  5  6  7  8  9 10 11 12
## 
## -- Imputation 43 --
## 
##   1  2  3  4  5  6  7  8  9
## 
## -- Imputation 44 --
## 
##   1  2  3  4  5  6  7
## 
## -- Imputation 45 --
## 
##   1  2  3  4  5  6  7
## 
## -- Imputation 46 --
## 
##   1  2  3  4  5  6  7
## 
## -- Imputation 47 --
## 
##   1  2  3  4  5  6  7  8
## 
## -- Imputation 48 --
## 
##   1  2  3  4  5  6  7  8  9 10
## 
## -- Imputation 49 --
## 
##   1  2  3  4  5  6  7  8  9 10 11
## 
## -- Imputation 50 --
## 
##   1  2  3  4  5  6  7

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(meanHAPPI ~ 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)  2.3239147 0.42966677  1.4781974 3.1696320     42 %
## GROUP1       0.2627454 0.38208462 -0.4877024 1.0131932     29 %
## WAVE        -0.1759976 0.19988559 -0.5701146 0.2181195     50 %
## BASELINE     0.6011878 0.07127906  0.4603353 0.7420404     58 %
## GROUP1:WAVE  0.2212059 0.25363472 -0.2777693 0.7201811     39 %
summary
##                results         se     (lower    upper) missInfo
## (Intercept)  2.3239147 0.42966677  1.4781974 3.1696320     42 %
## GROUP1       0.2627454 0.38208462 -0.4877024 1.0131932     29 %
## WAVE        -0.1759976 0.19988559 -0.5701146 0.2181195     50 %
## BASELINE     0.6011878 0.07127906  0.4603353 0.7420404     58 %
## GROUP1:WAVE  0.2212059 0.25363472 -0.2777693 0.7201811     39 %
library(pander)

Table

  results se (lower upper) missInfo
(Intercept) 2.324 0.4297 1.478 3.17 42 %
GROUP1 0.2627 0.3821 -0.4877 1.013 29 %
WAVE -0.176 0.1999 -0.5701 0.2181 50 %
BASELINE 0.6012 0.07128 0.4603 0.742 58 %
GROUP1:WAVE 0.2212 0.2536 -0.2778 0.7202 39 %

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(meanHAPPI ~ 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)  2.3148642  0.4312215  5.3681556 1.547380e-07
## GROUP1       0.2631203  0.4117904  0.6389666 5.230332e-01
## WAVE        -0.1759976  0.2141366 -0.8218940 4.118692e-01
## BASELINE     0.6029868  0.0682246  8.8382604 4.495575e-15
## GROUP1:WAVE  0.2211906  0.2752026  0.8037374 4.219727e-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) 2.315 0.4312 5.368 1.547e-07
GROUP1 0.2631 0.4118 0.639 0.523
WAVE -0.176 0.2141 -0.8219 0.4119
BASELINE 0.603 0.06822 8.838 4.496e-15
GROUP1:WAVE 0.2212 0.2752 0.8037 0.422

Check assumptions with Random Computations. Zailig fit with just one of the imputed data sets.

data1=MI$imputations[[1]]
library("Zelig")
zelig.fitdata1 <- zelig(meanHAPPI ~ 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.1636 -0.4236  0.1184  0.4608  1.7227 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  2.10382    0.31667   6.644 3.82e-10 ***
## GROUP1       0.39732    0.33566   1.184    0.238    
## WAVE         0.07394    0.15260   0.485    0.629    
## BASELINE     0.58239    0.04077  14.286  < 2e-16 ***
## GROUP1:WAVE  0.09896    0.21226   0.466    0.642    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.7076 on 173 degrees of freedom
## Multiple R-squared:  0.5615, Adjusted R-squared:  0.5514 
## F-statistic: 55.39 on 4 and 173 DF,  p-value: < 2.2e-16

Describe the meanHAPPI 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 5.03 1.32   5.33    5.16 0.99   2   7     5 -0.92
## meanHAPPI    2 86 5.14 1.20   5.33    5.25 0.99   1   7     6 -0.94
##           kurtosis   se
## BASELINE     -0.11 0.14
## meanHAPPI     0.91 0.13
## -------------------------------------------------------- 
## group: 1
##           vars  n mean   sd median trimmed  mad  min  max range  skew
## BASELINE     1 92 4.83 1.30   4.67    4.84 1.48 1.67 7.00  5.33 -0.09
## meanHAPPI    2 92 5.57 0.86   5.41    5.57 0.87 3.73 7.18  3.45  0.04
##           kurtosis   se
## BASELINE     -0.62 0.14
## meanHAPPI    -0.76 0.09

Create a plot that visualizes meanHAPPI 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(meanHAPPI ~ 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$meanHAPPI)) 
sel2 <- which(!is.na(data1$BASELINE))
data1$residual[intersect(sel1,sel2)] <- residual
qplot(GROUP, meanHAPPI, 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 measures on dataset Randomly Selected Imputed Data ======================================================== Graphing the Two-Way Interaction. Both meanHAPPI and the Residuals

# nlme package
with(data1, boxplot(meanHAPPI ~ WAVE + GROUP))

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

Comparing Basline to Wave 2 and 3 by Group.

fullModeldata1 <- lme(meanHAPPI ~ 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
##   386.4674 408.7399 -186.2337
## 
## Random effects:
##  Formula: ~1 | ID
##         (Intercept) Residual
## StdDev:   0.3279103 0.615683
## 
## Fixed effects: meanHAPPI ~ GROUP * WAVE + BASELINE 
##                 Value Std.Error DF   t-value p-value
## (Intercept) 2.1082930 0.3149563 87  6.693923  0.0000
## GROUP1      0.3971870 0.3045747 87  1.304071  0.1956
## WAVE        0.0739440 0.1346867 86  0.549007  0.5844
## BASELINE    0.5814968 0.0450081 86 12.919827  0.0000
## GROUP1:WAVE 0.0989299 0.1873527 86  0.528041  0.5988
##  Correlation: 
##             (Intr) GROUP1 WAVE   BASELI
## GROUP1      -0.515                     
## WAVE        -0.641  0.663              
## BASELINE    -0.719  0.022  0.000       
## GROUP1:WAVE  0.454 -0.922 -0.719  0.009
## 
## Standardized Within-Group Residuals:
##        Min         Q1        Med         Q3        Max 
## -4.0319571 -0.5406369  0.1254922  0.5620975  2.2726667 
## 
## Number of Observations: 178
## Number of Groups: 89

Another random selected imputation

data10=MI$imputations[[10]]
library("Zelig")
zelig.fitdata10 <- zelig(meanHAPPI ~ 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.0940 -0.3050  0.1150  0.3963  1.7364 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  2.313898   0.304369   7.602 1.77e-12 ***
## GROUP1       0.662858   0.322918   2.053   0.0416 *  
## WAVE         0.008335   0.146809   0.057   0.9548    
## BASELINE     0.552680   0.039131  14.124  < 2e-16 ***
## GROUP1:WAVE -0.062363   0.204220  -0.305   0.7605    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.6807 on 173 degrees of freedom
## Multiple R-squared:  0.5589, Adjusted R-squared:  0.5487 
## F-statistic:  54.8 on 4 and 173 DF,  p-value: < 2.2e-16

Describe the meanHAPPI 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 5.03 1.32   5.33    5.16 0.99   2   7     5 -0.92
## meanHAPPI    2 86 5.11 1.15   5.33    5.23 0.99   1   7     6 -1.12
##           kurtosis   se
## BASELINE     -0.11 0.14
## meanHAPPI     1.51 0.12
## -------------------------------------------------------- 
## group: 1
##           vars  n mean   sd median trimmed  mad  min  max range  skew
## BASELINE     1 92 4.81 1.31   4.67    4.82 1.48 1.67 7.00  5.33 -0.08
## meanHAPPI    2 92 5.56 0.81   5.39    5.55 0.91 3.65 7.06  3.41  0.06
##           kurtosis   se
## BASELINE     -0.65 0.14
## meanHAPPI    -0.60 0.08

Create a plot that visualizes meanHAPPI variable by the GROUP variable

library(ggplot2)
library(influence.ME)

Take a look at the residuals

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

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

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

Comparing Basline to Wave 2 and 3 by Group.

fullModeldata10 <- lme(meanHAPPI ~ 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
##   368.7387 391.0112 -177.3693
## 
## Random effects:
##  Formula: ~1 | ID
##         (Intercept)  Residual
## StdDev:   0.3679299 0.5612719
## 
## Fixed effects: meanHAPPI ~ GROUP * WAVE + BASELINE 
##                  Value  Std.Error DF   t-value p-value
## (Intercept)  2.3336566 0.30174177 87  7.733953  0.0000
## GROUP1       0.6623523 0.28146262 87  2.353251  0.0209
## WAVE         0.0083347 0.12278379 86  0.067881  0.9460
## BASELINE     0.5487524 0.04449871 86 12.331873  0.0000
## GROUP1:WAVE -0.0625964 0.17080839 86 -0.366472  0.7149
##  Correlation: 
##             (Intr) GROUP1 WAVE   BASELI
## GROUP1      -0.497                     
## WAVE        -0.610  0.654              
## BASELINE    -0.742  0.020  0.000       
## GROUP1:WAVE  0.427 -0.910 -0.719  0.015
## 
## Standardized Within-Group Residuals:
##        Min         Q1        Med         Q3        Max 
## -4.0153136 -0.4485420  0.1514783  0.5797057  2.2625328 
## 
## Number of Observations: 178
## Number of Groups: 89

Another random selected imputation

data15=MI$imputations[[15]]
library("Zelig")
zelig.fitdata15 <- zelig(meanHAPPI ~ 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.0005 -0.4127  0.1257  0.4281  1.9227 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  2.51343    0.34063   7.379  6.4e-12 ***
## GROUP1      -0.01613    0.36105  -0.045   0.9644    
## WAVE        -0.22471    0.16412  -1.369   0.1727    
## BASELINE     0.58729    0.04386  13.390  < 2e-16 ***
## GROUP1:WAVE  0.39339    0.22829   1.723   0.0866 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.761 on 173 degrees of freedom
## Multiple R-squared:  0.5337, Adjusted R-squared:  0.5229 
## F-statistic: 49.51 on 4 and 173 DF,  p-value: < 2.2e-16

Describe the meanHAPPI 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 5.03 1.32   5.33    5.16 0.99   2 7.00  5.00 -0.92
## meanHAPPI    2 86 5.13 1.21   5.33    5.23 0.99   1 7.15  6.15 -0.90
##           kurtosis   se
## BASELINE     -0.11 0.14
## meanHAPPI     0.91 0.13
## -------------------------------------------------------- 
## group: 1
##           vars  n mean   sd median trimmed  mad  min  max range  skew
## BASELINE     1 92 4.83 1.30   4.67    4.84 1.48 1.67 7.00  5.33 -0.09
## meanHAPPI    2 92 5.58 0.95   5.38    5.58 0.97 3.54 7.63  4.09  0.12
##           kurtosis   se
## BASELINE     -0.61 0.14
## meanHAPPI    -0.81 0.10

Create a plot that visualizes meanHAPPI variable by the GROUP variable

library(ggplot2)
library(influence.ME)

Take a look at the residuals

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

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

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

Comparing Basline to Wave 2 and 3 by Group.

fullModeldata15 <- lme(meanHAPPI ~ 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
##   409.6419 431.9144 -197.821
## 
## Random effects:
##  Formula: ~1 | ID
##         (Intercept)  Residual
## StdDev:   0.3960648 0.6371681
## 
## Fixed effects: meanHAPPI ~ GROUP * WAVE + BASELINE 
##                  Value Std.Error DF   t-value p-value
## (Intercept)  2.5203592 0.3382201 87  7.451832  0.0000
## GROUP1      -0.0163989 0.3183226 87 -0.051517  0.9590
## WAVE        -0.2247107 0.1393868 86 -1.612137  0.1106
## BASELINE     0.5859158 0.0495196 86 11.831989  0.0000
## GROUP1:WAVE  0.3933751 0.1938826 86  2.028935  0.0456
##  Correlation: 
##             (Intr) GROUP1 WAVE   BASELI
## GROUP1      -0.508                     
## WAVE        -0.618  0.657              
## BASELINE    -0.737  0.030  0.000       
## GROUP1:WAVE  0.443 -0.914 -0.719  0.002
## 
## Standardized Within-Group Residuals:
##        Min         Q1        Med         Q3        Max 
## -3.6638640 -0.4335494  0.1131387  0.5181568  2.5121361 
## 
## Number of Observations: 178
## Number of Groups: 89

Another randomly selected imputation

data25=MI$imputations[[25]]

library("Zelig")
zelig.fitdata25 <- zelig(meanHAPPI ~ GROUP * WAVE + BASELINE, random = ~1 | ID, data = data25,  model = "ls", cite = FALSE)
summary(zelig.fitdata25)
## 
## Call:
## lm(formula = formula, weights = weights, model = F, data = data)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -2.72035 -0.50978  0.09159  0.53515  2.59831 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  1.83165    0.36265   5.051 1.11e-06 ***
## GROUP1       0.10259    0.38442   0.267   0.7899    
## WAVE        -0.41689    0.17474  -2.386   0.0181 *  
## BASELINE     0.74450    0.04669  15.946  < 2e-16 ***
## GROUP1:WAVE  0.34682    0.24306   1.427   0.1554    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.8102 on 173 degrees of freedom
## Multiple R-squared:  0.6149, Adjusted R-squared:  0.606 
## F-statistic: 69.07 on 4 and 173 DF,  p-value: < 2.2e-16

Describe the meanHAPPI 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 5.03 1.32   5.33    5.16 0.99 2.00   7  5.00 -0.92
## meanHAPPI    2 86 4.95 1.39   5.26    5.09 1.10 0.36   7  6.64 -0.99
##           kurtosis   se
## BASELINE     -0.11 0.14
## meanHAPPI     0.72 0.15
## -------------------------------------------------------- 
## group: 1
##           vars  n mean   sd median trimmed  mad  min  max range  skew
## BASELINE     1 92 4.85 1.30   4.67    4.87 1.48 1.67 7.00  5.33 -0.14
## meanHAPPI    2 92 5.44 1.15   5.33    5.49 0.99 1.96 7.56  5.59 -0.45
##           kurtosis   se
## BASELINE     -0.61 0.14
## meanHAPPI     0.18 0.12

Create a plot that visualizes meanHAPPI variable by the GROUP variable

library(ggplot2)
library(influence.ME)

Take a look at the residuals

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

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

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

Comparing Basline to Wave 2 and 3 by Group.

fullModeldata25 <- lme(meanHAPPI ~ 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
##   439.0017 461.2742 -212.5009
## 
## Random effects:
##  Formula: ~1 | ID
##         (Intercept)  Residual
## StdDev:   0.1645619 0.7816507
## 
## Fixed effects: meanHAPPI ~ GROUP * WAVE + BASELINE 
##                  Value Std.Error DF   t-value p-value
## (Intercept)  1.8321890 0.3622907 87  5.057234  0.0000
## GROUP1       0.1025730 0.3778482 87  0.271466  0.7867
## WAVE        -0.4168945 0.1709938 86 -2.438067  0.0168
## BASELINE     0.7443911 0.0476666 86 15.616616  0.0000
## GROUP1:WAVE  0.3468178 0.2378470 86  1.458155  0.1484
##  Correlation: 
##             (Intr) GROUP1 WAVE   BASELI
## GROUP1      -0.555                     
## WAVE        -0.708  0.679              
## BASELINE    -0.662  0.025  0.000       
## GROUP1:WAVE  0.510 -0.944 -0.719 -0.002
## 
## Standardized Within-Group Residuals:
##        Min         Q1        Med         Q3        Max 
## -3.3036088 -0.6343760  0.1140090  0.6789153  3.2434702 
## 
## Number of Observations: 178
## Number of Groups: 89

Check assumptions on model without any imputations

Describe the meanHAPPI 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 5.03 1.32   5.33    5.16 0.99   2   7     5 -0.92
## meanHAPPI    2 59 5.23 1.24   5.33    5.37 0.99   1   7     6 -1.31
##           kurtosis   se
## BASELINE     -0.11 0.14
## meanHAPPI     1.72 0.16
## -------------------------------------------------------- 
## group: 1
##           vars  n mean   sd median trimmed  mad  min max range skew
## BASELINE     1 88 4.83 1.32   4.67    4.85 1.48 1.67   7  5.33 -0.1
## meanHAPPI    2 54 5.72 0.81   5.67    5.73 0.99 4.00   7  3.00  0.0
##           kurtosis   se
## BASELINE     -0.68 0.14
## meanHAPPI    -0.88 0.11

Create a plot that visualizes meanHAPPI variable by the GROUP variable

library(ggplot2)
library(influence.ME)

Take a look at the residuals

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

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

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

Comparing Basline to Wave 2 and 3 by Group.

fullModel <- lme(meanHAPPI ~ 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
##   242.5344 261.3738 -114.2672
## 
## Random effects:
##  Formula: ~1 | ID
##         (Intercept)  Residual
## StdDev:    0.635113 0.4393518
## 
## Fixed effects: meanHAPPI ~ GROUP * WAVE + BASELINE 
##                  Value Std.Error DF   t-value p-value
## (Intercept)  2.4085903 0.4119506 66  5.846794  0.0000
## GROUP1       0.2790397 0.3210326 66  0.869194  0.3879
## WAVE        -0.1457906 0.1318740 38 -1.105529  0.2759
## BASELINE     0.5707310 0.0689540 66  8.276985  0.0000
## GROUP1:WAVE  0.2483011 0.1923138 38  1.291125  0.2045
##  Correlation: 
##             (Intr) GROUP1 WAVE   BASELI
## GROUP1      -0.369                     
## WAVE        -0.414  0.562              
## BASELINE    -0.849  0.012 -0.029       
## GROUP1:WAVE  0.288 -0.825 -0.686  0.014
## 
## Standardized Within-Group Residuals:
##         Min          Q1         Med          Q3         Max 
## -2.31795119 -0.35922277  0.02430816  0.41902062  1.90237827 
## 
## Number of Observations: 109
## Number of Groups: 69

Table with confidence intervals

Table with confidence intervals

  lower est. upper
(Intercept) 1.605 2.409 3.212
GROUP1 -0.347 0.279 0.9051
WAVE -0.4066 -0.1458 0.115
BASELINE 0.4363 0.5707 0.7052
GROUP1:WAVE -0.132 0.2483 0.6286

Table with p-values

  Value Std.Error DF t-value p-value
(Intercept) 2.409 0.412 66 5.847 1.702e-07
GROUP1 0.279 0.321 66 0.8692 0.3879
WAVE -0.1458 0.1319 38 -1.106 0.2759
BASELINE 0.5707 0.06895 66 8.277 8.458e-12
GROUP1:WAVE 0.2483 0.1923 38 1.291 0.2045