Setup

library(haven)
library(data.table)
library(JWileymisc)
library(ggplot2)
library(ggpubr)
## Loading required package: magrittr
library(visreg)
library(VIM)
## Warning: package 'VIM' was built under R version 3.6.3
## Loading required package: colorspace
## Loading required package: grid
## Registered S3 methods overwritten by 'car':
##   method                          from
##   influence.merMod                lme4
##   cooks.distance.influence.merMod lme4
##   dfbeta.influence.merMod         lme4
##   dfbetas.influence.merMod        lme4
## VIM is ready to use. 
##  Since version 4.0.0 the GUI is in its own package VIMGUI.
## 
##           Please use the package to use the new (and old) GUI.
## Suggestions and bug-reports can be submitted at: https://github.com/alexkowa/VIM/issues
## 
## Attaching package: 'VIM'
## The following object is masked from 'package:datasets':
## 
##     sleep
library(mice)
## 
## Attaching package: 'mice'
## The following objects are masked from 'package:base':
## 
##     cbind, rbind
library(lme4)
## Loading required package: Matrix
library(lmerTest)
## Warning: package 'lmerTest' was built under R version 3.6.3
## 
## Attaching package: 'lmerTest'
## The following object is masked from 'package:lme4':
## 
##     lmer
## The following object is masked from 'package:stats':
## 
##     step
library(multilevelTools)
## Warning: package 'multilevelTools' was built under R version 3.6.3
library(emmeans)
## Welcome to emmeans.
## NOTE -- Important change from versions <= 1.41:
##     Indicator predictors are now treated as 2-level factors by default.
##     To revert to old behavior, use emm_options(cov.keep = character(0))
## read in data
db <- as.data.table(read_sav("B 19032020.sav")) # baseline data
data(aces_daily)
d <- as.data.table(aces_daily)

1. Generalized Linear Models 1

Instructions

  1. Fit one moderated, multiple linear regression. Pick from:
  • Outcome: positive or negative affect at baseline
  • Predictors: female and any other variable of your choosing
  • Make sure that you include the interaction between your predictor and female. Note: you will need to score positive or negative affect. For a refresher, see the content for Data Visualization 1.
  1. Conduct model diagnostics on your regression model. Where applicable, apply appropriate transformations and/or exclude outliers.
  2. Create a summary of your final model, after any relevant transformations or extreme values have been addressed, where applicable.
  • use modelTest() and APAStyler() to get a nice result.
  1. Create a graph to help visualize your results using visreg(). Note, you may want to adjust the fig.width and fig.height options to this R chunk to make the graph have the ratios you like.
  2. Briefly summarize your analysis steps, assumption checks, any changes performed, and interpret the interaction and regression coefficients, referencing your figure visualizing the result where appropriate.
## create positive affect variable
PosAffAVG <- db[, PosAffAVG := rowMeans(.SD, na.rm = TRUE),.SDcols = c("PANAS1", "PANAS3", "PANAS5", "PANAS9", "PANAS10", "PANAS12", "PANAS14", "PANAS14", "PANAS16", "PANAS17", "PANAS19")]

## put your regression model code using lm() here
PAMR <- lm(formula = PosAffAVG ~ female * stress, data = db)

## put your regression diagnostics code including plots using modelDiagnostics() here
PAMRd <- modelDiagnostics(PAMR, ev.perc = .005)
## Warning in .local(x, ...): singularity problem
## Warning in rq.fit.sfn(x, y, tau = tau, rhs = rhs, control = control, ...): tiny diagonals replaced with Inf when calling blkfct
## Warning in .local(x, ...): singularity problem
## Warning in rq.fit.sfn(x, y, tau = tau, rhs = rhs, control = control, ...): tiny diagonals replaced with Inf when calling blkfct
plot(PAMRd, ncol = 2, ask = FALSE)

## put your code to summarize the final model using modelTest() here
PAMRtest <- modelTest(PAMR)
knitr::kable(APAStyler(PAMRtest))
Term Est Type
(Intercept) 3.27 [-0.19, 6.74] Fixed Effects
female 1.62 [-1.96, 5.21] Fixed Effects
stress 0.00 [-0.15, 0.15] Fixed Effects
female:stress -0.07 [-0.22, 0.09] Fixed Effects
N (Observations) 56 Overall Model
logLik DF 5 Overall Model
logLik -50.83 Overall Model
AIC 111.66 Overall Model
BIC 121.78 Overall Model
F2 0.30 Overall Model
R2 0.23 Overall Model
Adj R2 0.19 Overall Model
female f2 = 0.02, p = .368 Effect Sizes
stress f2 = 0.00, p = .958 Effect Sizes
female:stress f2 = 0.01, p = .388 Effect Sizes 
## put your code to visualize your regression model using visreg() here
visreg(PAMR, xvar = "PosAffAVG", by = "stress")

confint(PAMR)
##                    2.5 %     97.5 %
## (Intercept)   -0.1873599 6.73541180
## female        -1.9647229 5.20848661
## stress        -0.1453695 0.15316173
## female:stress -0.2191422 0.08646222
## remove interaction as it was non-significant
PAMRni <- lm(formula = PosAffAVG ~ female + stress , data = db)
PAMRnitest <- modelTest(PAMRni)
knitr::kable(APAStyler(PAMRnitest))
Term Est Type
(Intercept) 4.73*** [ 3.88, 5.58] Fixed Effects
female 0.08 [-0.42, 0.58] Fixed Effects
stress -0.06*** [-0.09, -0.03] Fixed Effects
N (Observations) 56 Overall Model
logLik DF 4 Overall Model
logLik -51.23 Overall Model
AIC 110.47 Overall Model
BIC 118.57 Overall Model
F2 0.28 Overall Model
R2 0.22 Overall Model
Adj R2 0.19 Overall Model
female f2 = 0.00, p = .751 Effect Sizes
stress f2 = 0.26, p < .001 Effect Sizes 
visreg(PAMRni)

confint(PAMRni)
##                   2.5 %      97.5 %
## (Intercept)  3.87550059  5.58456606
## female      -0.42344273  0.58394860
## stress      -0.09124982 -0.02756731

Write Up

A multiple regression was run with an interaction between sex and stress, with the predictor variable being positive affect. As this interaction was non-significant, it was dropped. No outliers in the data were detected. The data was approximately normally distributed, as indicated by the density plots and the QQ deviates plot. Homogeneity of variance was not violated as shown by the spread of residuals on the scatterplot. We assume independence of observations. There was a statistically significant association between stress and positive affect, controlling for sex. A one unit higher stress score was associated with a 0.06 [95% CI -0.09 - -0.03] lower positive affect score, p < .001. There was no statistically significant association between sex and positive affect, controlling for stress (p < .001). This was a moderate effect (f2 = 0.26). Compared to males, females had 0.08 [95% CI -0.42 – 0.58] higher positive affect scores, this was non-significant association. Overall, stress and sex explained 22% of the variance in positive affect.

2. Generalized Linear Models 2

Instructions

  1. Conduct one multiple logistic regression predicting high selfesteem “SelfesteemHigh” as created in the GLM2 lecture content from two predictors of your choosing.
  2. Create a graph of the predicted probabilities of being in the SelfesteemHigh group by each of your two predictors. It should parallel the predicted probability plot here: http://joshuawiley.com/MonashHonoursStatistics/GLM2.html#glms-with-multiple-predictors
  3. Report and interpret the odds ratios, 95% confidence intervals for each predictor and present the average marginal effect for each predictor to interpret the results on the probability scale.
## create self esteem high
db[, SelfesteemHigh := as.integer(selfesteem >= 15)]

## run multiple logistic regression predicting high self esteem from two predictors
mlog <- glm(SelfesteemHigh ~ neuroticism + stress, data = db, family = binomial())
summary(mlog)
## 
## Call:
## glm(formula = SelfesteemHigh ~ neuroticism + stress, family = binomial(), 
##     data = db)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -2.1071  -0.8510   0.5195   0.8917   1.7081  
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)   
## (Intercept)  6.25688    1.92967   3.242  0.00119 **
## neuroticism -0.18238    0.16658  -1.095  0.27357   
## stress      -0.17052    0.06528  -2.612  0.00900 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 78.157  on 56  degrees of freedom
## Residual deviance: 64.775  on 54  degrees of freedom
## AIC: 70.775
## 
## Number of Fisher Scoring iterations: 4
## create graph of predicted probabilities
visreg(mlog, xvar = "stress",
       by = "neuroticism", breaks = c(4, 9),
       scale = "response", overlay = TRUE,
       partial = FALSE, rug = FALSE, gg = TRUE) +
  ylab("predicted probability of having high self esteem") +
  theme_pubr()

## pick h value for difference, store that in variable h
h <- .01

## make the original dataset
originaldata2 <- db[, .(stress, neuroticism)]

## make the incrased self esteem dataset (stress + h)
increasedStress <- db[, .(stress = stress + h, neuroticism)]

## make the incrased neuroticism dataset (neuroticism + h)
increasedNeuro <- db[, .(stress, neuroticism = neuroticism + h)]

## calculate original predicted probabilities
originaldata2$Prob <- predict(mlog, newdata = originaldata2,
                             type = "response")

## calculate increased predicted probabilities for self esteem
increasedStress$Prob <- predict(mlog, newdata = increasedStress,
                              type = "response")

## calculate increased predicted probabilities for neuroticism
increasedNeuro$Prob <- predict(mlog, newdata = increasedNeuro,
                              type = "response")

## calculate the difference in probabilities per unit
stress.diffprob <- (increasedStress$Prob - originaldata2$Prob) / h
neuro.diffprob <- (increasedNeuro$Prob - originaldata2$Prob) / h

## calculate the average marginal effects
mean(stress.diffprob)
## [1] -0.03284269
mean(neuro.diffprob)
## [1] -0.03512734
## report coeffients and confidence intervals
exp(coef(mlog))
## (Intercept) neuroticism      stress 
## 521.5888663   0.8332815   0.8432239
exp(confint(mlog))
## Waiting for profiling to be done...
##                  2.5 %       97.5 %
## (Intercept) 16.3758405 3.503105e+04
## neuroticism  0.5884599 1.143402e+00
## stress       0.7331093 9.504645e-01

Write Up

A logistic regression predicting high self esteem from stress and neuroticism showed that students who had low stress were expected to have 0.84 odds of being in the high self esteem group, [95% CI 0.73 - 9.50]. Students who had low neuroticism were expected to have 0.83 odds of being in the high self esteem group, [95% CI 0.59 - 1.14]. Results show that in the sample, holding neuroticism constant increasing stress by one unit on average is predicted to result in a -0.033 decline in the probability of being in the high self esteem group. Meanwhile, when stress is being held constant, increasing neuroticism by one unit on average is predicted to result in a -0.035 decline in the probability of being in the high self esteem group.

3. Missing Data

Instructions

  1. Using the missing dataset as in lecture, multiply impute the data, and use the with() and pool() functions to run a linear regression with any of the variables from the imputed dataset (i.e., Female, Age, Stress, PosAff, NegAff) and report pooled results.
  2. Code for the data setup is below. Briefly describe and interpret the linear regression results on the multiply imputed data.
options(scipen = 999)

## between person data, no missing
davg <- na.omit(d[, .(
 Female = factor(na.omit(Female)[1], levels = 0:1),
 Age = na.omit(Age)[1],
 STRESS = mean(STRESS, na.rm = TRUE),
 PosAff = mean(PosAff, na.rm = TRUE),
 NegAff = mean(NegAff, na.rm = TRUE)),
 by = UserID])

## create missing data
davgmiss <- copy(davg)
davgmiss[STRESS < 1, NegAff := NA]
davgmiss[STRESS > 4, PosAff := NA]

## random missingness on age
set.seed(1234)
davgmiss[, Age := ifelse(rbinom(
 .N, size = 1, prob = .1) == 1,
 NA, Age)]

## drop unneeded variables to make analysis easier
davgmiss[, UserID := NULL]

##evidence of convergence
mi.1 <- mice(davgmiss, m = 5, maxit = 3, defaultMethod = c("norm", "logreg", "polyreg", "polr"), seed = 1234, printFlag = FALSE)

## run an additional iterations
mi.1 <- mice.mids(mi.1, maxit = 20, printFlag = FALSE)
mi.reg <- with(mi.1, lm(NegAff ~ STRESS + PosAff))
mi.reg
## call :
## with.mids(data = mi.1, expr = lm(NegAff ~ STRESS + PosAff))
## 
## call1 :
## mice.mids(obj = mi.1, maxit = 20, printFlag = FALSE)
## 
## nmis :
## Female    Age STRESS PosAff NegAff 
##      0     19      0     30     40 
## 
## analyses :
## [[1]]
## 
## Call:
## lm(formula = NegAff ~ STRESS + PosAff)
## 
## Coefficients:
## (Intercept)       STRESS       PosAff  
##     0.77786      0.28801      0.01452  
## 
## 
## [[2]]
## 
## Call:
## lm(formula = NegAff ~ STRESS + PosAff)
## 
## Coefficients:
## (Intercept)       STRESS       PosAff  
##     1.08261      0.24576     -0.05112  
## 
## 
## [[3]]
## 
## Call:
## lm(formula = NegAff ~ STRESS + PosAff)
## 
## Coefficients:
## (Intercept)       STRESS       PosAff  
##      0.9174       0.2670      -0.0137  
## 
## 
## [[4]]
## 
## Call:
## lm(formula = NegAff ~ STRESS + PosAff)
## 
## Coefficients:
## (Intercept)       STRESS       PosAff  
##     1.05268      0.25099     -0.04699  
## 
## 
## [[5]]
## 
## Call:
## lm(formula = NegAff ~ STRESS + PosAff)
## 
## Coefficients:
## (Intercept)       STRESS       PosAff  
##     0.69397      0.27431      0.07079
pool(mi.reg)
## Class: mipo    m = 5 
##          term m     estimate         ubar            b            t dfcom
## 1 (Intercept) 5  0.904902312 0.0122777391 0.0285518184 0.0465399211   187
## 2      STRESS 5  0.265205487 0.0003246815 0.0002965485 0.0006805397   187
## 3      PosAff 5 -0.005299381 0.0009591494 0.0025225840 0.0039862503   187
##          df      riv    lambda       fmi
## 1  6.411087 2.790594 0.7361891 0.7922529
## 2 12.549327 1.096022 0.5229059 0.5842711
## 3  6.001400 3.156026 0.7593855 0.8128471
m.mireg <- summary(pool(mi.reg), conf.int=TRUE)
print(m.mireg)
##          term     estimate  std.error   statistic        df
## 1 (Intercept)  0.904902312 0.21573113  4.19458381  6.411087
## 2      STRESS  0.265205487 0.02608716 10.16613247 12.549327
## 3      PosAff -0.005299381 0.06313676 -0.08393496  6.001400
##           p.value      2.5 %    97.5 %
## 1 0.0049390532811  0.3851245 1.4246801
## 2 0.0000002029388  0.2086412 0.3217698
## 3 0.9358379359955 -0.1597807 0.1491820
pool.r.squared(mi.reg)
##           est     lo 95     hi 95 fmi
## R^2 0.5957116 0.4745241 0.6969421 NaN

Write up

For a one unit change in stress, negative affect increases by 0.26 (95% CI 0.21, 0.32) (p < .001). Meanwhile, for a one unit change in positive affect, negative affect decreases slightly by 0.005 (95% CI -0.16, 0.15), but this assoociation was not significant. 59.57% of the variability in negative affect can be explained by positive affect and stress.

4. Linear Mixed Models 1

Instructions

  1. run an intercept only linear mixed model with PosAff ratings as the outcome variable including a fixed and random intercept by person (the id variable is called UserID).
  2. Describe the number of observations and participants included in the analysis
  3. Report & interpret the intraclass correlation coefficient (ICC).
  4. Conduct model diagnostics and briefly comment on whether the homogeneity of variance and normality assumptions appear to be met or not (even if not met, you do not need to address it in this lab report section).
  5. Finally, report the intercept with 95% confidence interval.
## linear mixed model with PosAff ratings as the outcome variable including a fixed and random intercept by person
ri.m <- lmer(PosAff ~ 1 +(1 | UserID), data = d, REML = TRUE)
summary(ri.m)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: PosAff ~ 1 + (1 | UserID)
##    Data: d
## 
## REML criterion at convergence: 14794.6
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -4.3450 -0.6468 -0.0341  0.6169  4.0577 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  UserID   (Intercept) 0.6289   0.7930  
##  Residual             0.5290   0.7273  
## Number of obs: 6399, groups:  UserID, 191
## 
## Fixed effects:
##              Estimate Std. Error        df t value            Pr(>|t|)    
## (Intercept)   2.67866    0.05811 189.83075    46.1 <0.0000000000000002 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## report ICC
iccMixed("PosAff", id = "UserID", data = d)
##         Var     Sigma       ICC
## 1:   UserID 0.6289049 0.5431473
## 2: Residual 0.5289851 0.4568527
## conduct model diagnostics
md <- modelDiagnostics(ri.m, ev.perc = .001)
plot(md, ask = FALSE, ncol = 2, nrow = 2)

## report confidence intervals
ri.ci <- confint(ri.m, method = "profile", oldNames = FALSE)
## Computing profile confidence intervals ...
ri.ci
##                           2.5 %    97.5 %
## sd_(Intercept)|UserID 0.7157462 0.8795099
## sigma                 0.7147055 0.7402970
## (Intercept)           2.5644941 2.7928351

Write up

An intercept only linear mixed model was fit to 6399 positive affect ratings from 191 people. The assumption of normality was met as indicated by the denisty plot. However, the residuals plot reveals that homogeneity of variance was violated at the high and low ends of the data. The intraclass correlation coefficient was 0.54. This means that slightly more of the variation exists between people rather than within people. Therefore, about 54% of the total varation is between person. The remaining 46% variation can be attributed to differences between days. The fixed effect intercept shows that the average [95% CI] posisitive affect was 2.67 [2.56, 2.79].

5. Linear Mixed Models 2

Instructions

  1. Create a between and within version of STRESS using meanDeviations()
  2. Fit a linear mixed model using lmer() with both the between and within version of STRESS included as fixed effects and both a random intercept and random slope for the within version of STRESS by UserID.
  3. Check model diagnostics and where appropriate use transformations or exclude outliers / extreme values. You only need to do this once (i.e., if after removing extreme values, there are still more extreme values, you do not need to address them, just do it once). Also, deal with any singularity / non convergence issues, simplifying the model if needed to achieve convergence. Keep all code to document your process.
  4. Make a summary of the final model (after any cleaning).
  5. Briefly report the results from the LMM (interpret the fixed effects along with confidence intervals / p-values).
## make between person and within person stress
d[, c("Bstress", "Wstress") := meanDeviations(STRESS), by = UserID]

## run linear mixed model with both the between and within version of stress included as fixed effects and both a random intercept and random slope for the within version of stress
lmm <- lmer(SUPPORT ~ Bstress + Wstress + (Wstress | UserID), data = d)
lmm.md <- modelDiagnostics(lmm, ev.perc = .001)
plot(lmm.md, ask = FALSE, ncol = 2, nrow = 2)

## remove outliers 
lmm.md$extremeValues
##         SUPPORT UserID Index                        EffectType
##  1:  0.33744722      9   273                         Residuals
##  2:  0.00000000     12   394                         Residuals
##  3:  0.00000000     12   400                         Residuals
##  4:  1.36898422     19   631                         Residuals
##  5:  1.52065041     19   643                         Residuals
##  6:  0.00000000     75  2588                         Residuals
##  7:  0.76026374     75  2609                         Residuals
##  8: 10.00000000     97  3329                         Residuals
##  9:  0.00000000     98  3373                         Residuals
## 10:  0.44581106    108  3716                         Residuals
## 11:  0.70943453    117  4022                         Residuals
## 12:  9.39224694    135  4647                         Residuals
## 13:  9.39741464    154  5326                         Residuals
## 14:  3.71539797     67  2296    Random Effect UserID : Wstress
## 15:  0.70475312     67  2299    Random Effect UserID : Wstress
## 16:  0.05905785     67  2302    Random Effect UserID : Wstress
## 17:  3.51222852     67  2305    Random Effect UserID : Wstress
## 18:  1.64901116     67  2308    Random Effect UserID : Wstress
## 19:  2.04100785     67  2311    Random Effect UserID : Wstress
## 20:  8.51653426     67  2314    Random Effect UserID : Wstress
## 21:  9.10268704     67  2317    Random Effect UserID : Wstress
## 22:  3.67570322     67  2320    Random Effect UserID : Wstress
## 23:  3.05623782     67  2323    Random Effect UserID : Wstress
## 24:  2.37536564     67  2326    Random Effect UserID : Wstress
## 25:  5.96875113     67  2329    Random Effect UserID : Wstress
## 26:  8.02684440     92  3146    Random Effect UserID : Wstress
## 27:  4.57739883     92  3149    Random Effect UserID : Wstress
## 28:  0.00000000     92  3152    Random Effect UserID : Wstress
## 29:  0.08087353     92  3155    Random Effect UserID : Wstress
## 30:  6.90382678     92  3158    Random Effect UserID : Wstress
## 31:  9.79360987     92  3161    Random Effect UserID : Wstress
## 32:  4.12676114     92  3164    Random Effect UserID : Wstress
## 33:  1.12680047     92  3167    Random Effect UserID : Wstress
## 34:  5.79055586     92  3170    Random Effect UserID : Wstress
## 35:  5.71980549     92  3173    Random Effect UserID : Wstress
## 36:  7.75583221     92  3176    Random Effect UserID : Wstress
## 37:  8.02684440     92  3146 Multivariate Random Effect UserID
## 38:  4.57739883     92  3149 Multivariate Random Effect UserID
## 39:  0.00000000     92  3152 Multivariate Random Effect UserID
## 40:  0.08087353     92  3155 Multivariate Random Effect UserID
## 41:  6.90382678     92  3158 Multivariate Random Effect UserID
## 42:  9.79360987     92  3161 Multivariate Random Effect UserID
## 43:  4.12676114     92  3164 Multivariate Random Effect UserID
## 44:  1.12680047     92  3167 Multivariate Random Effect UserID
## 45:  5.79055586     92  3170 Multivariate Random Effect UserID
## 46:  5.71980549     92  3173 Multivariate Random Effect UserID
## 47:  7.75583221     92  3176 Multivariate Random Effect UserID
##         SUPPORT UserID Index                        EffectType
d.noev <- d[-unique(lmm.md$extremeValues$Index)]

## rerun analysis
lmm.noev <- lmer(SUPPORT ~ Bstress + Wstress + (Wstress | UserID), data = d.noev)
md.noev <- modelDiagnostics(lmm.noev, ev.perc = .001)
plot(md.noev, ask = FALSE, ncol = 2, nrow = 2)

summary(lmm.noev)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: SUPPORT ~ Bstress + Wstress + (Wstress | UserID)
##    Data: d.noev
## 
## REML criterion at convergence: 9293.6
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.6396 -0.6387  0.0627  0.6142  2.8556 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev. Corr 
##  UserID   (Intercept) 3.49668  1.8699        
##           Wstress     0.03714  0.1927   -0.25
##  Residual             3.49728  1.8701        
## Number of obs: 2139, groups:  UserID, 189
## 
## Fixed effects:
##              Estimate Std. Error        df t value             Pr(>|t|)
## (Intercept)   5.77813    0.26583 187.58976  21.736 < 0.0000000000000002
## Bstress      -0.14960    0.09600 188.37212  -1.558                0.121
## Wstress      -0.22128    0.02626 141.10200  -8.428   0.0000000000000367
##                
## (Intercept) ***
## Bstress        
## Wstress     ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##         (Intr) Bstrss
## Bstress -0.845       
## Wstress -0.033 -0.040
confint(lmm.noev, oldNames = FALSE)
## Computing profile confidence intervals ...
##                                     2.5 %      97.5 %
## sd_(Intercept)|UserID           1.6691191  2.08178532
## cor_Wstress.(Intercept)|UserID -0.5063808  0.03139872
## sd_Wstress|UserID               0.1285737  0.25398664
## sigma                           1.8106827  1.93279612
## (Intercept)                     5.2560741  6.29979788
## Bstress                        -0.3380796  0.03899403
## Wstress                        -0.2733782 -0.16980607

Write Up

A linear mixed model was fit to 2139 support scores from 189 people. The fixed effect intercept revealed that the estimated average [95% CI] support score was 5.78 [5.26, 6.30], when Bstress and Wstress are zero (p < .001). There was also a non-significant fixed effect of average stress on support, such that each one unit higher average stress people had was associated with -0.14 95% CI = [-0.34, -0.04] (p = .12). There was also a significant fixed effect of within person stress on support, such that on days where people were one unit higher on stress than their own average, people were expected to have -0.22 95% CI = [-0.27, -0.16] lower support that same day, on average (p < .001). Finally, there was a negative correlation between the random intercept and slope, -0.25 indicating that compared to the population averages, people who had a relatively higher support when support was 0 also tended to have a more negative within person association between stress and support.

6. Interactions and Moderation for LMMs

Instructions

  1. Pick a continuous, within person variable as the outcome (any variable measured more than once per person). Also choose at least two variables that you will use as predictors/explanatory variables (may be between and/or within). At least one of your two predictors/explanatory variables must be a continuous variable (e.g., not Female, EDU).
  2. Fit a linear mixed model using lmer() and test whether your predictors interact with each other by including at least a fixed effect of the interaction term. You may include random slopes or not. Ensure that you have a model that converges and does not have a singularity warning (if you’re unable to solve these, you may need to pick different variables). Note, model diagnostics are not required for this lab report section, only convergence and no singularity issues.
  3. Regardless of whether the interaction is statistically significant or not, create a graph of the interaction. Create a graph of the interaction using visreg().
  4. Regardless of whether the interaction is statistically significant or not, calculate the simple slopes, confidence intervals, and p-values for the simple slopes (e.g., using emtrends()). Briefly report the simple slopes.
## fit a linear mixed model
m <- lmer(STRESS ~ SUPPORT * PosAff + (1 | UserID), data = d)

## calculate simple slopes
mem <- emtrends(m, var = "PosAff", at = list(support = c(6.05 - 2.41, 6.05 + 2.41)), lmer.df = "satterthwaite")
summary(mem, infer=TRUE)
##  SUPPORT PosAff PosAff.trend     SE   df lower.CL upper.CL t.ratio p.value
##     5.38   2.67        -1.22 0.0573 2121    -1.33    -1.11 -21.330 <.0001 
## 
## Degrees-of-freedom method: satterthwaite 
## Confidence level used: 0.95
## create simple slopes graph
egltable("SUPPORT", data = d[!duplicated(UserID)])
##                 M (SD)
## 1: SUPPORT 6.05 (2.41)
visreg(m, xvar = "PosAff", by = "SUPPORT", overlay=TRUE, breaks = c(6.05 - 2.41, 6.05 + 2.41), partial = FALSE, rug = FALSE)

mem <- emtrends(m, var = "PosAff", at = list(SUPPORT = c(6.05 - 2.41, 6.05 + 2.41)), lmer.df = "satterthwaite")
summary(mem, infer=TRUE)
##  SUPPORT PosAff PosAff.trend     SE   df lower.CL upper.CL t.ratio p.value
##     3.64   2.67        -1.24 0.0673 2148    -1.37    -1.11 -18.434 <.0001 
##     8.46   2.67        -1.19 0.0702 2162    -1.33    -1.05 -16.921 <.0001 
## 
## Degrees-of-freedom method: satterthwaite 
## Confidence level used: 0.95

Write up

The simple slope of positive affect on stress was tested and graphed at low (M - 1SD) and high (M + 1SD) values of support. The results showed that there was a significant negative association between positive affect and stress at both high and low levels of support.

7. Model Comparisons for LMMs

Instructions

  1. Pick a continuous, within person variable as the outcome (any variable measured more than once per person). Also choose a continuous, within person predictor/explanatory variable.
  2. Fit a linear mixed model using lmer() and fit models that have an intercept only, linear effect of your predictor, quadratic effect of your predictor, and cubic effect of your predictor. Use AIC and BIC values to compare models. Which model is the best? You may have your explanatory variable as a fixed only or fixed and random effect. Ensure that you have models that converge and do not have a singularity warning (e.g., dropping random effects if needed; if you’re unable to solve these, you may need to pick different variables). Note, model diagnostics are not required for this lab report section, only convergence and no singularity issues.
  3. Report and interpret the conditional and marginal f2 effect size(s) for your predictor from your final model. If your best model is an intercept only model, then pick a different model (e.g., linear effect of your predictor model) to report and interpret the conditional and marginal f2 effect sizes for your predictor.
## run models that have an intercept only, linear effect of your predictor, quadratic effect of your predictor, and cubic effect of your predictor
m0 <- lmer(PosAff ~ 1 + (1 | UserID), data = d, REML = FALSE)
m1 <- lmer(PosAff ~ poly(STRESS, 1) + (1 | UserID), data = d[!is.na(STRESS)], REML = FALSE)
m2 <- lmer(PosAff ~ poly(STRESS, 2) + (1 | UserID), data = d[!is.na(STRESS)], REML = FALSE)
m3 <- lmer(PosAff ~ poly(STRESS, 3) + (1 | UserID), data = d[!is.na(STRESS)], REML = FALSE)

## check number of observations are equivalent
nobs(m0)
## [1] 6399
nobs(m1)
## [1] 6399
nobs(m2)
## [1] 6399
nobs(m3)
## [1] 6399
## calculate AIC and BIC
AIC(m0, m1, m2, m3)
##    df      AIC
## m0  3 14796.71
## m1  4 13372.31
## m2  5 13340.58
## m3  6 13341.93
BIC(m0, m1, m2, m3)
##    df      BIC
## m0  3 14817.00
## m1  4 13399.36
## m2  5 13374.40
## m3  6 13382.51
## report conditional and marginal f2
modelPerformance(m2)
## $Performance
##     Model Estimator N_Obs     N_Groups      AIC     BIC        LL LLDF
## 1: merMod        ML  6399 UserID (191) 13340.58 13374.4 -6665.292    5
##        Sigma MarginalR2 ConditionalR2 MarginalF2 ConditionalF2
## 1: 0.6484209  0.1483331     0.6203202   0.174168      1.633799
## 
## attr(,"class")
## [1] "modelPerformance.merMod" "modelPerformance"

Write up

Linear mixed models that have an intercept only, linear effect of stress, quadratic effect of stress, and cubic effect of stress were fit. AIC and BIC show that of the four models, the quadratic model is the best fit. The marginal effect is medium (f2 = 0.17) and the conditional effect is large (f2 = 1.63).