Loading the dataset

data.test4 <- read.csv("/Volumes/TOSHIBA EXT/Dropbox/ADULT STUDY/adult_study011615.csv")
# Load the psych package
library(psych)
data.test4$HAPPIPERMA17 <- apply(data.test4[, c ("HAPPI1" ,"HAPPI2", "HAPPI3", "PERMA17")], 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 HAPPIPERMA17 and ID Group and wave from data.test4 and create a new #dataset with only those variables.
data <- data.test4[,c("ID", "GROUP", "wave", "HAPPIPERMA17")]
#Use dcast to cnage from long-format data to wide format data
data <- dcast(data, ID + GROUP ~ wave, mean, value.var = "HAPPIPERMA17")
# 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 HAPPIPERMA17 and wave 2 and 3 of HAPPIPERMA17. 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", "HAPPIPERMA17", "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  8  9 10 11
## 
## -- Imputation 2 --
## 
##   1  2  3  4  5  6  7
## 
## -- Imputation 3 --
## 
##   1  2  3  4  5  6  7  8  9 10 11 12 13 14
## 
## -- Imputation 4 --
## 
##   1  2  3  4  5  6  7  8  9
## 
## -- Imputation 5 --
## 
##   1  2  3  4  5  6  7
## 
## -- Imputation 6 --
## 
##   1  2  3  4  5  6  7  8
## 
## -- Imputation 7 --
## 
##   1  2  3  4  5  6  7
## 
## -- Imputation 8 --
## 
##   1  2  3  4  5
## 
## -- Imputation 9 --
## 
##   1  2  3  4  5  6  7  8  9
## 
## -- Imputation 10 --
## 
##   1  2  3  4  5  6
## 
## -- Imputation 11 --
## 
##   1  2  3  4  5  6  7  8  9
## 
## -- Imputation 12 --
## 
##   1  2  3  4  5
## 
## -- Imputation 13 --
## 
##   1  2  3  4  5  6  7  8  9 10
## 
## -- Imputation 14 --
## 
##   1  2  3  4  5  6  7  8
## 
## -- Imputation 15 --
## 
##   1  2  3  4  5  6  7
## 
## -- Imputation 16 --
## 
##   1  2  3  4  5  6  7  8
## 
## -- Imputation 17 --
## 
##   1  2  3  4  5  6  7  8
## 
## -- 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
## 
## -- Imputation 22 --
## 
##   1  2  3  4  5  6  7  8
## 
## -- Imputation 23 --
## 
##   1  2  3  4  5  6  7  8  9
## 
## -- Imputation 24 --
## 
##   1  2  3  4  5  6  7  8  9 10
## 
## -- Imputation 25 --
## 
##   1  2  3  4  5  6  7  8  9 10 11 12 13
## 
## -- 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  8  9
## 
## -- Imputation 29 --
## 
##   1  2  3  4  5  6
## 
## -- Imputation 30 --
## 
##   1  2  3  4  5  6  7  8
## 
## -- Imputation 31 --
## 
##   1  2  3  4  5  6  7  8  9 10 11
## 
## -- Imputation 32 --
## 
##   1  2  3  4  5  6  7  8
## 
## -- Imputation 33 --
## 
##   1  2  3  4  5  6  7  8  9 10 11
## 
## -- Imputation 34 --
## 
##   1  2  3  4  5  6  7
## 
## -- Imputation 35 --
## 
##   1  2  3  4  5  6  7  8
## 
## -- Imputation 36 --
## 
##   1  2  3  4  5  6  7
## 
## -- Imputation 37 --
## 
##   1  2  3  4  5  6  7  8  9
## 
## -- 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  9 10
## 
## -- 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
## 
## -- 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  8  9
## 
## -- Imputation 46 --
## 
##   1  2  3  4  5  6
## 
## -- Imputation 47 --
## 
##   1  2  3  4  5  6  7  8  9 10 11 12
## 
## -- Imputation 48 --
## 
##   1  2  3  4  5  6  7
## 
## -- Imputation 49 --
## 
##   1  2  3  4  5  6
## 
## -- Imputation 50 --
## 
##   1  2  3  4  5  6  7  8  9 10 11 12 13

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(HAPPIPERMA17 ~ 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)  2.48198450 0.46543611  1.5668592 3.3971098     36 %
## GROUP1       0.37459855 0.41461270 -0.4394285 1.1886256     27 %
## WAVE        -0.08768055 0.20589283 -0.4929526 0.3175915     42 %
## BASELINE     0.58851954 0.07106886  0.4483176 0.7287214     52 %
## GROUP1:WAVE  0.15113732 0.25878170 -0.3570797 0.6593543     29 %

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(HAPPIPERMA17 ~ 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)  2.46960634 0.46778945  5.2793118 2.084649e-07
## GROUP1       0.37497965 0.45036052  0.8326211 4.052618e-01
## WAVE        -0.08768055 0.22394418 -0.3915286 6.956167e-01
## BASELINE     0.59077967 0.06744318  8.7596648 2.050128e-15
## GROUP1:WAVE  0.15113244 0.28630014  0.5278811 5.977106e-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(HAPPIPERMA17 ~ 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.7183 -0.5046  0.1132  0.4774  2.4616 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  2.73563    0.43220   6.330 2.04e-09 ***
## GROUP1       0.77719    0.45173   1.720   0.0871 .  
## WAVE         0.07238    0.20538   0.352   0.7250    
## BASELINE     0.51598    0.05208   9.908  < 2e-16 ***
## GROUP1:WAVE -0.18317    0.28570  -0.641   0.5223    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.9523 on 173 degrees of freedom
## Multiple R-squared:  0.3819, Adjusted R-squared:  0.3676 
## F-statistic: 26.72 on 4 and 173 DF,  p-value: < 2.2e-16

Describe the HAPPIPERMA17 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.48 1.38   5.75    5.62 1.11 2.25 7.75  5.50 -0.89
## HAPPIPERMA17    2 86 5.67 1.34   5.97    5.81 1.16 0.75 7.80  7.05 -1.16
##              kurtosis   se
## BASELINE         0.00 0.15
## HAPPIPERMA17     1.65 0.14
## -------------------------------------------------------- 
## group: 1
##              vars  n mean   sd median trimmed  mad  min  max range  skew
## BASELINE        1 92 5.29 1.38   5.25    5.32 1.48 1.75 7.75  6.00 -0.28
## HAPPIPERMA17    2 92 6.08 1.01   6.00    6.10 0.74 3.59 9.01  5.43 -0.13
##              kurtosis   se
## BASELINE        -0.25 0.14
## HAPPIPERMA17     0.44 0.11

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

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

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

Comparing Basline to Wave 2 and 3 by Group.

fullModeldata1 <- lme(HAPPIPERMA17 ~ 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
##   494.3811 516.6536 -240.1905
## 
## Random effects:
##  Formula: ~1 | ID
##         (Intercept)  Residual
## StdDev:   0.3753402 0.8605729
## 
## Fixed effects: HAPPIPERMA17 ~ GROUP * WAVE + BASELINE 
##                  Value Std.Error DF   t-value p-value
## (Intercept)  2.7510576 0.4311952 87  6.380075  0.0000
## GROUP1       0.7769417 0.4218697 87  1.841663  0.0689
## WAVE         0.0723763 0.1882588 86  0.384451  0.7016
## BASELINE     0.5131618 0.0559677 86  9.168898  0.0000
## GROUP1:WAVE -0.1833580 0.2618890 86 -0.700136  0.4857
##  Correlation: 
##             (Intr) GROUP1 WAVE   BASELI
## GROUP1      -0.514                     
## WAVE        -0.655  0.669              
## BASELINE    -0.711  0.011  0.000       
## GROUP1:WAVE  0.460 -0.931 -0.719  0.015
## 
## Standardized Within-Group Residuals:
##         Min          Q1         Med          Q3         Max 
## -3.92664844 -0.45276149  0.08007613  0.47010014  2.54713066 
## 
## Number of Observations: 178
## Number of Groups: 89

Another random imputation

data10=MI$imputations[[10]]
library("Zelig")
zelig.fitdata10 <- zelig(HAPPIPERMA17 ~ 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.6494 -0.3322  0.0970  0.4359  2.0699 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  2.28721    0.36625   6.245 3.18e-09 ***
## GROUP1       0.37829    0.38308   0.987    0.325    
## WAVE         0.01534    0.17413   0.088    0.930    
## BASELINE     0.60573    0.04410  13.735  < 2e-16 ***
## GROUP1:WAVE  0.10623    0.24221   0.439    0.662    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.8074 on 173 degrees of freedom
## Multiple R-squared:  0.5401, Adjusted R-squared:  0.5295 
## F-statistic: 50.79 on 4 and 173 DF,  p-value: < 2.2e-16

Describe the HAPPIPERMA17 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.48 1.38   5.75    5.62 1.11 2.25 7.75   5.5 -0.89
## HAPPIPERMA17    2 86 5.63 1.32   5.75    5.76 1.27 0.75 7.75   7.0 -1.13
##              kurtosis   se
## BASELINE         0.00 0.15
## HAPPIPERMA17     1.68 0.14
## -------------------------------------------------------- 
## group: 1
##              vars  n mean   sd median trimmed  mad  min  max range  skew
## BASELINE        1 92 5.33 1.38    5.5    5.38 1.48 1.75 7.75  6.00 -0.35
## HAPPIPERMA17    2 92 6.08 0.98    6.0    6.11 0.74 2.73 8.49  5.75 -0.39
##              kurtosis   se
## BASELINE        -0.23 0.14
## HAPPIPERMA17     0.71 0.10

Create a plot that visualizes HAPPIPERMA17 variable by the GROUP variable

library(ggplot2)
library(influence.ME)

Take a look at the residuals

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

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

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

Comparing Basline to Wave 2 and 3 by Group.

fullModeldata10 <- lme(HAPPIPERMA17 ~ 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
##   425.1098 447.3823 -205.5549
## 
## Random effects:
##  Formula: ~1 | ID
##         (Intercept)  Residual
## StdDev:   0.4817449 0.6336993
## 
## Fixed effects: HAPPIPERMA17 ~ GROUP * WAVE + BASELINE 
##                 Value Std.Error DF   t-value p-value
## (Intercept) 2.3259144 0.3637763 87  6.393803  0.0000
## GROUP1      0.3767049 0.3222281 87  1.169063  0.2456
## WAVE        0.0153392 0.1386280 86  0.110650  0.9122
## BASELINE    0.5986597 0.0512346 86 11.684670  0.0000
## GROUP1:WAVE 0.1065849 0.1928442 86  0.552700  0.5819
##  Correlation: 
##             (Intr) GROUP1 WAVE   BASELI
## GROUP1      -0.485                     
## WAVE        -0.572  0.645              
## BASELINE    -0.771  0.036  0.000       
## GROUP1:WAVE  0.421 -0.898 -0.719 -0.014
## 
## Standardized Within-Group Residuals:
##         Min          Q1         Med          Q3         Max 
## -4.03175710 -0.43331841  0.09074686  0.49022247  1.99579838 
## 
## Number of Observations: 178
## Number of Groups: 89

Another random imputation

data15=MI$imputations[[15]]
library("Zelig")
zelig.fitdata15 <- zelig(HAPPIPERMA17 ~ 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.3830 -0.4156  0.1329  0.4937  2.1319 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  2.71850    0.37759   7.200 1.77e-11 ***
## GROUP1       0.48162    0.39382   1.223   0.2230    
## WAVE        -0.31472    0.17902  -1.758   0.0805 .  
## BASELINE     0.59890    0.04563  13.125  < 2e-16 ***
## GROUP1:WAVE  0.15076    0.24901   0.605   0.5457    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.8301 on 173 degrees of freedom
## Multiple R-squared:  0.5362, Adjusted R-squared:  0.5255 
## F-statistic:    50 on 4 and 173 DF,  p-value: < 2.2e-16

Describe the HAPPIPERMA17 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.48 1.38   5.75    5.62 1.11 2.25 7.75   5.5 -0.89
## HAPPIPERMA17    2 86 5.53 1.38   5.75    5.67 1.11 0.75 7.75   7.0 -1.09
##              kurtosis   se
## BASELINE         0.00 0.15
## HAPPIPERMA17     1.09 0.15
## -------------------------------------------------------- 
## group: 1
##              vars  n mean   sd median trimmed  mad  min  max range  skew
## BASELINE        1 92 5.31 1.36   5.38    5.35 1.30 1.75 7.75  6.00 -0.34
## HAPPIPERMA17    2 92 6.13 0.93   6.25    6.17 0.79 3.68 8.59  4.91 -0.29
##              kurtosis   se
## BASELINE        -0.14 0.14
## HAPPIPERMA17     0.07 0.10

Create a plot that visualizes HAPPIPERMA17 variable by the GROUP variable

library(ggplot2)
library(influence.ME)

Take a look at the residuals

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

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

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

Comparing Basline to Wave 2 and 3 by Group.

fullModeldata15 <- lme(HAPPIPERMA17 ~ 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
##   437.8522 460.1247 -211.9261
## 
## Random effects:
##  Formula: ~1 | ID
##         (Intercept)  Residual
## StdDev:   0.4663966 0.6724331
## 
## Fixed effects: HAPPIPERMA17 ~ GROUP * WAVE + BASELINE 
##                  Value Std.Error DF   t-value p-value
## (Intercept)  2.7176408 0.3768178 87  7.212081  0.0000
## GROUP1       0.4816453 0.3388803 87  1.421285  0.1588
## WAVE        -0.3147201 0.1471014 86 -2.139477  0.0352
## BASELINE     0.5990581 0.0525053 86 11.409477  0.0000
## GROUP1:WAVE  0.1507582 0.2046151 86  0.736789  0.4633
##  Correlation: 
##             (Intr) GROUP1 WAVE   BASELI
## GROUP1      -0.487                     
## WAVE        -0.586  0.651              
## BASELINE    -0.763  0.030  0.000       
## GROUP1:WAVE  0.425 -0.906 -0.719 -0.005
## 
## Standardized Within-Group Residuals:
##         Min          Q1         Med          Q3         Max 
## -4.10510975 -0.51939733  0.06502394  0.55460482  2.18045198 
## 
## Number of Observations: 178
## Number of Groups: 89

Another random imputation

data25=MI$imputations[[25]]

library("Zelig")
zelig.fitdata25 <- zelig(HAPPIPERMA17 ~ 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.3649 -0.4546  0.0694  0.5456  2.2968 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  2.05485    0.38227   5.375 2.45e-07 ***
## GROUP1       0.21619    0.39940   0.541    0.589    
## WAVE        -0.17408    0.18160  -0.959    0.339    
## BASELINE     0.66513    0.04608  14.435  < 2e-16 ***
## GROUP1:WAVE  0.35772    0.25262   1.416    0.159    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.842 on 173 degrees of freedom
## Multiple R-squared:  0.5768, Adjusted R-squared:  0.567 
## F-statistic: 58.94 on 4 and 173 DF,  p-value: < 2.2e-16

Describe the HAPPIPERMA17 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.48 1.38   5.75    5.62 1.11 2.25 7.75   5.5 -0.89
## HAPPIPERMA17    2 86 5.44 1.42   5.75    5.58 1.11 0.75 7.75   7.0 -1.01
##              kurtosis   se
## BASELINE         0.00 0.15
## HAPPIPERMA17     0.65 0.15
## -------------------------------------------------------- 
## group: 1
##              vars  n mean   sd median trimmed  mad  min  max range  skew
## BASELINE        1 92 5.31 1.38   5.49    5.35 1.48 1.75 7.75  6.00 -0.33
## HAPPIPERMA17    2 92 6.08 1.04   6.00    6.10 1.11 3.91 8.47  4.55 -0.11
##              kurtosis   se
## BASELINE        -0.21 0.14
## HAPPIPERMA17    -0.64 0.11

Create a plot that visualizes HAPPIPERMA17 variable by the GROUP variable

library(ggplot2)
library(influence.ME)

Take a look at the residuals

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

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

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

Comparing Basline to Wave 2 and 3 by Group.

fullModeldata25 <- lme(HAPPIPERMA17 ~ 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
##   447.1819 469.4544 -216.5909
## 
## Random effects:
##  Formula: ~1 | ID
##         (Intercept)  Residual
## StdDev:   0.4140574 0.7195042
## 
## Fixed effects: HAPPIPERMA17 ~ GROUP * WAVE + BASELINE 
##                  Value Std.Error DF   t-value p-value
## (Intercept)  2.0744863 0.3810316 87  5.444395  0.0000
## GROUP1       0.2159692 0.3574629 87  0.604172  0.5473
## WAVE        -0.1740841 0.1573987 86 -1.106007  0.2718
## BASELINE     0.6615459 0.0513682 86 12.878505  0.0000
## GROUP1:WAVE  0.3574705 0.2189648 86  1.632548  0.1062
##  Correlation: 
##             (Intr) GROUP1 WAVE   BASELI
## GROUP1      -0.491                     
## WAVE        -0.620  0.660              
## BASELINE    -0.738  0.009  0.000       
## GROUP1:WAVE  0.433 -0.918 -0.719  0.016
## 
## Standardized Within-Group Residuals:
##        Min         Q1        Med         Q3        Max 
## -3.3781718 -0.4689206  0.1057731  0.5436523  2.6747643 
## 
## Number of Observations: 178
## Number of Groups: 89

Check assumptions on model without any imputations

Describe the HAPPIPERMA17 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.48 1.38   5.75    5.62 1.11 2.25 7.75   5.5 -0.89
## HAPPIPERMA17    2 59 5.69 1.40   6.00    5.87 1.11 0.75 7.75   7.0 -1.48
##              kurtosis   se
## BASELINE         0.00 0.15
## HAPPIPERMA17     2.39 0.18
## -------------------------------------------------------- 
## group: 1
##              vars  n mean   sd median trimmed  mad  min  max range  skew
## BASELINE        1 88 5.30 1.39   5.38    5.35 1.48 1.75 7.75   6.0 -0.32
## HAPPIPERMA17    2 54 6.22 0.80   6.25    6.24 0.74 4.25 7.75   3.5 -0.26
##              kurtosis   se
## BASELINE        -0.23 0.15
## HAPPIPERMA17     0.10 0.11

Create a plot that visualizes HAPPIPERMA17 variable by the GROUP variable

library(ggplot2)
library(influence.ME)

Take a look at the residuals

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

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

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

Comparing Basline to Wave 2 and 3 by Group.

fullModel <- lme(HAPPIPERMA17 ~ 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
##   267.5962 286.4356 -126.7981
## 
## Random effects:
##  Formula: ~1 | ID
##         (Intercept)  Residual
## StdDev:   0.7307607 0.4814194
## 
## Fixed effects: HAPPIPERMA17 ~ GROUP * WAVE + BASELINE 
##                  Value Std.Error DF   t-value p-value
## (Intercept)  2.5293563 0.4798912 66  5.270687  0.0000
## GROUP1       0.4675688 0.3566186 66  1.311117  0.1944
## WAVE        -0.0954297 0.1449653 38 -0.658293  0.5143
## BASELINE     0.5673585 0.0751682 66  7.547850  0.0000
## GROUP1:WAVE  0.1536837 0.2113500 38  0.727153  0.4716
##  Correlation: 
##             (Intr) GROUP1 WAVE   BASELI
## GROUP1      -0.344                     
## WAVE        -0.390  0.556              
## BASELINE    -0.864  0.002 -0.027       
## GROUP1:WAVE  0.270 -0.816 -0.686  0.015
## 
## Standardized Within-Group Residuals:
##         Min          Q1         Med          Q3         Max 
## -2.48680420 -0.30247786  0.07239194  0.38319867  1.98864961 
## 
## Number of Observations: 109
## Number of Groups: 69
Table with P-value

|             |       Value|  Std.Error|  DF|     t-value|    p-value|
|:------------|-----------:|----------:|---:|-----------:|----------:|
|(Intercept)  |   2.5293563|  0.4798912|  66|   5.2706870|  0.0000016|
|GROUP1       |   0.4675688|  0.3566186|  66|   1.3111174|  0.1943615|
|WAVE         |  -0.0954297|  0.1449653|  38|  -0.6582935|  0.5143184|
|BASELINE     |   0.5673585|  0.0751682|  66|   7.5478501|  0.0000000|
|GROUP1:WAVE  |   0.1536837|  0.2113500|  38|   0.7271526|  0.4715888|

Table with confidence intervals

est. lower upper
(Intercept) 2.5293563 1.5934561 3.4652565
GROUP1 0.4675688 -0.2279209 1.1630586
WAVE -0.0954297 -0.3820867 0.1912273
BASELINE 0.5673585 0.4207628 0.7139541
GROUP1:WAVE 0.1536837 -0.2642436 0.5716110