MLM - Ex3 - Eldad Aviv
Eldad Aviv
2024-08-16
MLM - Ex3
Eldad Aviv: 206836165
Pre Process
# Load Libraries
library(tidyverse)
library(datawizard)
library(lmerTest)
library(performance)
library(parameters)# Load Data
data <- read.csv("ChangeInRiskyBehaviorAdolescence_Hoffman_Chap.7.csv")
head(data)## PersonID Attitude12 Age12 Age13 Age14 Age15 Age16 Age17
## 1 1 5.000000 11.69729 13.21557 14.04451 15.29666 15.85848 17.04552
## 2 2 3.045675 11.88690 12.75442 13.86911 14.90883 15.78037 16.76942
## 3 3 3.768105 12.32108 12.90831 13.99549 15.06845 16.12700 16.88830
## 4 4 5.000000 12.03376 13.09281 14.16630 14.80565 15.87344 16.79336
## 5 5 4.311243 11.92280 13.03784 13.90019 14.96241 15.94619 17.02968
## 6 6 3.808912 11.77549 12.78307 14.09661 15.21296 16.14227 17.02965
## Age18 Risky12 Risky13 Risky14 Risky15 Risky16 Risky17 Risky18
## 1 17.92486 27.36784 26.71078 25.00833 22.33165 22.93829 22.38455 22.19030
## 2 17.98165 15.04131 17.70607 13.87121 17.51945 21.01569 19.84318 30.03761
## 3 17.78550 19.64890 11.74143 11.53835 15.48973 11.56835 20.27700 20.08044
## 4 17.84713 14.92558 14.90911 19.55625 19.83930 15.06753 20.55906 18.84474
## 5 17.79814 10.24811 19.73088 13.40048 14.61820 16.78653 24.40688 27.58493
## 6 18.08528 14.52408 19.41424 22.92023 20.17371 17.35936 22.97153 24.50401
## Monitor12 Monitor13 Monitor14 Monitor15 Monitor16 Monitor17 Monitor18
## 1 2.689663 2.602434 2.458114 2.837886 2.629586 2.500052 3.041118
## 2 3.816063 4.557213 3.763838 4.122347 3.717130 3.070058 3.980447
## 3 4.855208 5.000000 4.386148 4.436663 3.444854 3.670691 3.705285
## 4 3.425048 2.710094 3.480533 3.759563 3.120889 3.125168 3.449266
## 5 3.251903 3.337643 3.277787 3.677860 3.156562 3.581186 3.865481
## 6 3.513176 3.779001 3.735469 3.738414 3.637157 3.276559 2.961469## Data Prep
# Creating long (stacked) format of the data:
data.long <- data |>
pivot_longer(
cols = starts_with(c("Age", "Risky", "Monitor")),
names_pattern = "(.*)([0-9][0-9])",
names_to = c(".value", "Age_rounded")
) |>
mutate(
# we centered age to 18:
YearsPre18 = Age - 18
)
head(data.long)## # A tibble: 6 × 7
## PersonID Attitude12 Age_rounded Age Risky Monitor YearsPre18
## <int> <dbl> <chr> <dbl> <dbl> <dbl> <dbl>
## 1 1 5 12 11.7 27.4 2.69 -6.30
## 2 1 5 13 13.2 26.7 2.60 -4.78
## 3 1 5 14 14.0 25.0 2.46 -3.96
## 4 1 5 15 15.3 22.3 2.84 -2.70
## 5 1 5 16 15.9 22.9 2.63 -2.14
## 6 1 5 17 17.0 22.4 2.50 -0.9541.For each predictor, examine if it should be treated as a time varrying predictor
Monitor ICC
lmer(Monitor ~ 1 + (1 | PersonID), data = data.long) |>
icc()## # Intraclass Correlation Coefficient
##
## Adjusted ICC: 0.688
## Unadjusted ICC: 0.688Report Monitor ICC:
69% of variance was due to BP mean differences. thus we will address it’s random component in our analysis.
# Visually:
data.long |>
mutate(
PersonID = factor(PersonID) |> fct_reorder(Monitor, .fun = mean)
) |>
ggplot(aes(PersonID, Monitor)) +
geom_point() +
stat_summary(geom = "point", fill = "red", size = 2, shape = 21)## No summary function supplied, defaulting to `mean_se()`
Age ICC
Age_ICC <- lmer(Age ~ 1 + (1 | PersonID), data = data.long)## boundary (singular) fit: see help('isSingular')# Summary of the model
summary(Age_ICC)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: Age ~ 1 + (1 | PersonID)
## Data: data.long
##
## REML criterion at convergence: 5922.7
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.73146 -0.93560 0.00465 0.93874 1.66445
##
## Random effects:
## Groups Name Variance Std.Dev.
## PersonID (Intercept) 0.000 0.000
## Residual 4.017 2.004
## Number of obs: 1400, groups: PersonID, 200
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 1.500e+01 5.356e-02 1.399e+03 280.1 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')Report
The variance of zero in the random effect for age by person suggests that age does not vary significantly between individuals after accounting for the fixed effect of age. This indicates that a simpler model with age as a fixed effect only may be more appropriate for this analysis.
Predictors Configuration
data.long <- data.long |>
group_by(PersonID) |>
mutate(
Monitor_PM = mean(Monitor, na.rm = TRUE), # person mean
Monitor_WP = Monitor - Monitor_PM, # person mean centering!
) |>
select(PersonID, Age, Monitor, Monitor_PM, Monitor_WP, Risky)Modeling
First we will check to proption of Risky variance explained by individual differences.
M0 <- lmer(Risky ~ 1 + (1 | PersonID),
data = data.long)
icc(M0)## # Intraclass Correlation Coefficient
##
## Adjusted ICC: 0.388
## Unadjusted ICC: 0.38839% of the variance is due to BP differences
M1a <- lmer(Risky ~ Monitor_PM + Monitor_WP + # mood new vars as fixed effects
(1 | PersonID),
data = data.long)
model_parameters(M1a, ci_method = "S")## # Fixed Effects
##
## Parameter | Coefficient | SE | 95% CI | t | df | p
## ----------------------------------------------------------------------------
## (Intercept) | 26.96 | 1.36 | [24.27, 29.65] | 19.77 | 198.00 | < .001
## Monitor PM | -2.46 | 0.44 | [-3.32, -1.60] | -5.64 | 198.00 | < .001
## Monitor WP | 1.42 | 0.33 | [ 0.77, 2.07] | 4.29 | 1199.00 | < .001
##
## # Random Effects
##
## Parameter | Coefficient
## --------------------------------------
## SD (Intercept: PersonID) | 3.03
## SD (Residual) | 4.12##
## Uncertainty intervals (equal-tailed) and p-values (two-tailed) computed
## using a Wald t-distribution with Satterthwaite approximation.Monitor_PM: The person mean of monitoring (Monitor_PM) has a negative effect on risky behavior (coef = -2.46, p < 0.001). This suggests that, on average, higher levels of monitoring across the assessment period are associated with lower levels of risky behavior
Monitor_WP: The within-person deviation of monitoring from their personal mean (Monitor_WP) has a positive effect on risky behavior(coef = 1.42, p<0.001). This indicates that, on times when a girl’s monitoring score is higher than her personal average, she tends to exhibit more risky behavior.
M1b <- lmer(Risky ~ Monitor_PM + Monitor_WP +
(Monitor_WP | PersonID), # setting WP mood changes as random + COV with the intercept
data = data.long)
model_parameters(M1b, ci_method = "S")## # Fixed Effects
##
## Parameter | Coefficient | SE | 95% CI | t | df | p
## ---------------------------------------------------------------------------
## (Intercept) | 26.79 | 1.36 | [24.10, 29.48] | 19.65 | 197.98 | < .001
## Monitor PM | -2.41 | 0.44 | [-3.27, -1.55] | -5.52 | 197.96 | < .001
## Monitor WP | 1.55 | 0.44 | [ 0.67, 2.43] | 3.49 | 146.00 | < .001
##
## # Random Effects
##
## Parameter | Coefficient
## --------------------------------------------------
## SD (Intercept: PersonID) | 3.07
## SD (Monitor_WP: PersonID) | 3.98
## Cor (Intercept~Monitor_WP: PersonID) | 0.06
## SD (Residual) | 3.89##
## Uncertainty intervals (equal-tailed) and p-values (two-tailed) computed
## using a Wald t-distribution with Satterthwaite approximation.anova(M1a, M1b)## refitting model(s) with ML (instead of REML)## Data: data.long
## Models:
## M1a: Risky ~ Monitor_PM + Monitor_WP + (1 | PersonID)
## M1b: Risky ~ Monitor_PM + Monitor_WP + (Monitor_WP | PersonID)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## M1a 5 8258.5 8284.7 -4124.2 8248.5
## M1b 7 8231.9 8268.6 -4109.0 8217.9 30.575 2 2.295e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Interpretation:
The model comparison using ANOVA shows that M1b, which includes random slopes for Monitor_WP, significantly improves model fit compared to M1a. The addition of random slopes reduces the AIC and BIC values, and the chi-square test is highly significant (p> 0.001). This indicates that the more complex model (M1b) better captures the variability in risky behavior due to individual differences in how deviations in monitoring impact behavior.
In summary, M1b is preferred over M1a as it accounts for individual differences in the effect of monitoring deviations, providing a more accurate representation of the relationship between monitoring and risky behavior.
Compare Models: Pseudo-R2
VarCorr(M1a)## Groups Name Std.Dev.
## PersonID (Intercept) 3.0251
## Residual 4.1224VarCorr(M1b)## Groups Name Std.Dev. Corr
## PersonID (Intercept) 3.0694
## Monitor_WP 3.9841 0.058
## Residual 3.88751 - (3.8875 / 4.1224)^2## [1] 0.1107159This result indicates that the more complex model (M1b) explains approximately 11% more of the variance in risky behavior compared to the simpler model (M1a) when considering residual variance.