title: “EG_Happiness” output: html_document
Load Packages
library(arm); library(lmerTest); library(psych)
## Loading required package: MASS
## Loading required package: Matrix
## Loading required package: lme4
## Loading required package: Rcpp
##
## arm (Version 1.7-03, built: 2014-4-27)
##
## Working directory is /private/var/folders/9c/88_yf73x0ys02zvn5hkvvq6r0000gn/T/RtmpmxWTU5
##
## KernSmooth 2.23 loaded
## Copyright M. P. Wand 1997-2009
##
## Attaching package: 'lmerTest'
##
## The following object is masked from 'package:lme4':
##
## lmer
##
## The following object is masked from 'package:stats':
##
## step
##
##
## Attaching package: 'psych'
##
## The following objects are masked from 'package:arm':
##
## logit, rescale, sim
Loading Data
setwd("/Volumes/TOSHIBA EXT/Dropbox/Schools Study Data/Emily Griffith Data and R Scripts")
data <- read.csv("EmilyGriffith_all.csv")
Creating School ID as ID
data$ID <- data$Q1
#Create scale scores
data$meanHAPPI <- apply(data[, c("HAPPI.1_1", "HAPPI.2_1", "HAPPI.3_1", "HAPPI.4_1")], 1, mean, na.rm = TRUE)
#Means or plotting
data$baseline <- ifelse(data$Time<4,0,1)
pdata <- tapply(data[,"meanHAPPI"], data[,3], mean, na.rm=TRUE)
plot(pdata, type="l")
Baseline Model
M0 <- lmer(meanHAPPI ~ 1 + (1|ID), data=data)
fixef(M0)
## (Intercept)
## 4.904
Confidence Intervals
confint(M0)
## Computing profile confidence intervals ...
## 2.5 % 97.5 %
## .sig01 0.0000 0.5517
## .sigma 0.8642 1.1717
## (Intercept) 4.7073 5.1018
Model 1
M1 <- update(M0, .~. + Time, REML=FALSE)
fixef(M1)
## (Intercept) Time
## 5.09514 -0.08583
Confidence Intervals for Model 1
confint(M1)
## Computing profile confidence intervals ...
## Warning: convergence code 3 from bobyqa: bobyqa -- a trust region step
## failed to reduce q
## 2.5 % 97.5 %
## .sig01 0.0000 0.5508
## .sigma 0.8612 1.1678
## (Intercept) 4.6210 5.5694
## Time -0.2801 0.1083
Model 2
M2 <- update(M1, .~. + baseline)
fixef(M2)
## (Intercept) Time baseline
## 5.4187 -0.2730 0.8164
confint(M2)
## Computing profile confidence intervals ...
## 2.5 % 97.5 %
## .sig01 0.00000 0.54475
## .sigma 0.84421 1.14613
## (Intercept) 4.86356 5.97426
## Time -0.53250 -0.01427
## baseline 0.05271 1.58209
Model 3
M3 <- update(M2, .~. + I(Time^2))
fixef(M3)
## (Intercept) Time baseline I(Time^2)
## 7.8677 -3.1113 -1.6384 0.7103
Confidence Intervals for Model 3
confint(M3)
## Computing profile confidence intervals ...
## 2.5 % 97.5 %
## .sig01 0.0000 0.6077
## .sigma 0.7620 1.0497
## (Intercept) 6.5266 9.1947
## Time -4.5463 -1.6586
## baseline -3.0464 -0.2130
## I(Time^2) 0.3519 1.0641
Pdata <- tapply(data[,"meanHAPPI"], data[,3], mean, na.rm=TRUE)
# Add random noise to time to better see the points of interest
data$TimeJIT <- data$Time+runif(126, min=-.1, max=.1)
with(data, plot(TimeJIT, meanHAPPI, col="grey", pch="*"))
lines(pdata, col="red", lwd=2)
Intervention was given between 3 and 4 (Time 1 is not relevant because the data is not clean)
boxplot(data$meanHAPPI~Time, data=data, notch=F, col=(c("red","blue", "green")), main="HAPPINESS", xlab="HAPPINESS", xlim = c(), ylim = c(2, 9 ), yaxs = "i")
```
Note that the echo = FALSE parameter was added to the code chunk to prevent printing of the R code that generated the plot.