library(MCMCglmm)
## Loading required package: Matrix
## Loading required package: coda
## Loading required package: ape
library(lme4)
library(tidyverse)
## ── Attaching packages ───────────────────────────────────────────────────────────────── tidyverse 1.2.1 ──
## ✔ ggplot2 2.2.1.9000 ✔ purrr 0.2.4
## ✔ tibble 1.4.2 ✔ dplyr 0.7.4
## ✔ tidyr 0.8.0 ✔ stringr 1.3.0
## ✔ readr 1.1.1 ✔ forcats 0.3.0
## ── Conflicts ──────────────────────────────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ tidyr::expand() masks Matrix::expand()
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag() masks stats::lag()
library(tidyverse)
library(lme4)
library(corrr)
library(jmRtools)
library(tidyLPA)
library(kableExtra)
library(sjPlot)
library(broom)
library(broom.mixed)
library(konfound)
Engagement
lme4
m0i <- lmer(rm_engagement ~ 1 + pre_interest + gender_female + (1 | participant_ID), data = d_red)
d_BLUP_2 <- ranef(m0i) %>%
pluck(1) %>%
rownames_to_column("participant_ID") %>%
mutate(participant_ID = as.integer(participant_ID)) %>%
rename(rm_engagement_BLUP = `(Intercept)`)
d_ind_level_2 <- distinct(d, participant_ID, post_interest, program_ID, .keep_all = TRUE)
d_for_m2 <- left_join(d_ind_level_2, d_BLUP_2, by = "participant_ID")
d_for_m2 <- filter(d_for_m2, !is.na(post_interest) & !is.na(gender_female) & !is.na(urm) & !is.na(rm_engagement_BLUP))
m0ii <- lm(post_interest ~ 1 + rm_engagement_BLUP + gender_female + pre_interest, data = d_for_m2)
summary(m0i)
r2beta(m0i)
summary(m0ii)
r2beta(m0ii)
sd(d_red$rm_engagement)
sd(d_red$post_interest, na.rm = T)
## Linear mixed model fit by REML ['lmerMod']
## Formula:
## rm_engagement ~ 1 + pre_interest + gender_female + (1 | participant_ID)
## Data: d_red
##
## REML criterion at convergence: 5728.7
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -4.3115 -0.5279 0.1058 0.5844 3.7095
##
## Random effects:
## Groups Name Variance Std.Dev.
## participant_ID (Intercept) 0.3336 0.5776
## Residual 0.4083 0.6390
## Number of obs: 2714, groups: participant_ID, 180
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 2.55029 0.17009 14.994
## pre_interest 0.10847 0.05056 2.145
## gender_female -0.06856 0.09075 -0.756
##
## Correlation of Fixed Effects:
## (Intr) pr_ntr
## pre_interst -0.925
## gender_feml -0.335 0.069
## Effect Rsq upper.CL lower.CL
## 1 Model 0.021 0.034 0.012
## 2 pre_interest 0.018 0.029 0.009
## 3 gender_female 0.002 0.008 0.000
##
## Call:
## lm(formula = post_interest ~ 1 + rm_engagement_BLUP + gender_female +
## pre_interest, data = d_for_m2)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.97365 -0.35127 0.03355 0.42297 1.79135
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.24076 0.23397 5.303 4.45e-07 ***
## rm_engagement_BLUP 0.46171 0.10915 4.230 4.25e-05 ***
## gender_female -0.14041 0.11879 -1.182 0.239
## pre_interest 0.62788 0.06884 9.121 8.37e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.6976 on 137 degrees of freedom
## Multiple R-squared: 0.4309, Adjusted R-squared: 0.4184
## F-statistic: 34.58 on 3 and 137 DF, p-value: < 2.2e-16
##
## Effect Rsq upper.CL lower.CL
## 1 Model 0.431 0.543 0.325
## 4 pre_interest 0.378 0.492 0.264
## 2 rm_engagement_BLUP 0.116 0.227 0.036
## 3 gender_female 0.010 0.069 0.000
## [1] 0.8643042
## [1] 0.9360477
MCMCglmm
prior = list(R = list(V = diag(2), nu = 0.002),
G = list(G1 = list(V = diag(2), nu = 2, alpha.mu = rep(0,2),
alpha.V = diag(25^2,2,2))))
m2 <- MCMCglmm(fixed = cbind(rm_engagement, post_interest) ~ trait - 1 +
trait:gender_female +
trait:pre_interest,
random=~us(trait):participant_ID,
rcov=~idh(trait):units, # this denotes response vars varying across response rows
family = rep("gaussian",2),
data=as.data.frame(d_red),
prior=prior,
burnin = 3000,
nitt = 23000,
thin = 20,
verbose=TRUE)
summary(m2)
m2_cor <- m2$VCV[,"traitpost_interest:traitrm_engagement.participant_ID"]/
(sqrt(m2$VCV[,"traitpost_interest:traitpost_interest.participant_ID"])*
sqrt(m2$VCV[,"traitrm_engagement:traitrm_engagement.participant_ID"]))
broom::tidy(m2)
posterior.mode(m2_cor)
plot(m2_cor)
HPDinterval(m2_cor)
##
## MCMC iteration = 0
##
## MCMC iteration = 1000
##
## MCMC iteration = 2000
##
## MCMC iteration = 3000
##
## MCMC iteration = 4000
##
## MCMC iteration = 5000
##
## MCMC iteration = 6000
##
## MCMC iteration = 7000
##
## MCMC iteration = 8000
##
## MCMC iteration = 9000
##
## MCMC iteration = 10000
##
## MCMC iteration = 11000
##
## MCMC iteration = 12000
##
## MCMC iteration = 13000
##
## MCMC iteration = 14000
##
## MCMC iteration = 15000
##
## MCMC iteration = 16000
##
## MCMC iteration = 17000
##
## MCMC iteration = 18000
##
## MCMC iteration = 19000
##
## MCMC iteration = 20000
##
## MCMC iteration = 21000
##
## MCMC iteration = 22000
##
## MCMC iteration = 23000
##
## Iterations = 3001:22981
## Thinning interval = 20
## Sample size = 1000
##
## DIC: -22307.52
##
## G-structure: ~us(trait):participant_ID
##
## post.mean l-95% CI
## traitrm_engagement:traitrm_engagement.participant_ID 0.3422 0.2599
## traitpost_interest:traitrm_engagement.participant_ID 0.1355 0.0515
## traitrm_engagement:traitpost_interest.participant_ID 0.1355 0.0515
## traitpost_interest:traitpost_interest.participant_ID 0.6191 0.4157
## u-95% CI eff.samp
## traitrm_engagement:traitrm_engagement.participant_ID 0.4194 652.9
## traitpost_interest:traitrm_engagement.participant_ID 0.2163 1000.0
## traitrm_engagement:traitpost_interest.participant_ID 0.2163 1000.0
## traitpost_interest:traitpost_interest.participant_ID 0.8670 1000.0
##
## R-structure: ~idh(trait):units
##
## post.mean l-95% CI u-95% CI eff.samp
## traitrm_engagement.units 4.092e-01 3.881e-01 4.321e-01 1000
## traitpost_interest.units 9.159e-07 8.607e-07 9.626e-07 1178
##
## Location effects: cbind(rm_engagement, post_interest) ~ trait - 1 + trait:gender_female + trait:pre_interest
##
## post.mean l-95% CI u-95% CI eff.samp
## traitrm_engagement 2.55048 2.00351 3.01261 1000.0
## traitpost_interest 1.77651 -0.27295 3.40217 919.5
## traitrm_engagement:gender_female -0.06622 -0.25922 0.12726 1160.9
## traitpost_interest:gender_female -0.11897 -0.46559 0.28071 1000.0
## traitrm_engagement:pre_interest 0.10781 -0.04526 0.25560 836.2
## traitpost_interest:pre_interest 0.45492 -0.03864 1.04628 1277.6
## pMCMC
## traitrm_engagement <0.001 ***
## traitpost_interest 0.066 .
## traitrm_engagement:gender_female 0.472
## traitpost_interest:gender_female 0.450
## traitrm_engagement:pre_interest 0.144
## traitpost_interest:pre_interest 0.068 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## effect term estimate std.error
## 1 fixed traitrm_engagement 2.55047805 0.26258836
## 2 fixed traitpost_interest 1.77650556 0.94205882
## 3 fixed traitrm_engagement:gender_female -0.06622246 0.09689615
## 4 fixed traitpost_interest:gender_female -0.11897292 0.18484084
## 5 fixed traitrm_engagement:pre_interest 0.10780600 0.07790819
## 6 fixed traitpost_interest:pre_interest 0.45492379 0.27754389
## var1
## 0.2976175
## lower upper
## var1 0.1466879 0.4360675
## attr(,"Probability")
## [1] 0.95
Negative affect
lme4
m0i <- lmer(negative_affect ~ 1 + pre_interest + gender_female + (1 | participant_ID), data = d_red)
d_BLUP_2 <- ranef(m0i) %>%
pluck(1) %>%
rownames_to_column("participant_ID") %>%
mutate(participant_ID = as.integer(participant_ID)) %>%
rename(negative_affect_BLUP = `(Intercept)`)
d_ind_level_2 <- distinct(d, participant_ID, post_interest, program_ID, .keep_all = TRUE)
d_for_m2 <- left_join(d_ind_level_2, d_BLUP_2, by = "participant_ID")
d_for_m2 <- filter(d_for_m2, !is.na(post_interest) & !is.na(gender_female) & !is.na(urm) & !is.na(negative_affect_BLUP))
m0ii <- lm(post_interest ~ 1 + negative_affect_BLUP + gender_female + pre_interest, data = d_for_m2)
summary(m0i)
r2beta(m0i)
summary(m0ii)
r2beta(m0ii)
sqrt(.007)
## Linear mixed model fit by REML ['lmerMod']
## Formula:
## negative_affect ~ 1 + pre_interest + gender_female + (1 | participant_ID)
## Data: d_red
##
## REML criterion at convergence: 6364.5
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.6000 -0.4527 -0.1262 0.3479 3.9044
##
## Random effects:
## Groups Name Variance Std.Dev.
## participant_ID (Intercept) 0.6154 0.7845
## Residual 0.5039 0.7099
## Number of obs: 2714, groups: participant_ID, 180
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 1.97608 0.22729 8.694
## pre_interest -0.01052 0.06755 -0.156
## gender_female -0.14074 0.12138 -1.160
##
## Correlation of Fixed Effects:
## (Intr) pr_ntr
## pre_interst -0.925
## gender_feml -0.337 0.070
## Effect Rsq upper.CL lower.CL
## 1 Model 0.008 0.016 0.003
## 3 gender_female 0.008 0.016 0.002
## 2 pre_interest 0.000 0.003 0.000
##
## Call:
## lm(formula = post_interest ~ 1 + negative_affect_BLUP + gender_female +
## pre_interest, data = d_for_m2)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.7796 -0.4793 0.1402 0.3830 2.0170
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.31997 0.24706 5.343 3.72e-07 ***
## negative_affect_BLUP -0.08470 0.08580 -0.987 0.325
## gender_female -0.18527 0.12533 -1.478 0.142
## pre_interest 0.61314 0.07285 8.417 4.58e-14 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.7391 on 137 degrees of freedom
## Multiple R-squared: 0.3611, Adjusted R-squared: 0.3471
## F-statistic: 25.81 on 3 and 137 DF, p-value: 2.683e-13
##
## Effect Rsq upper.CL lower.CL
## 1 Model 0.361 0.482 0.253
## 4 pre_interest 0.341 0.458 0.226
## 3 gender_female 0.016 0.082 0.000
## 2 negative_affect_BLUP 0.007 0.061 0.000
## [1] 0.083666
MCMCglmm
prior = list(R = list(V = diag(2), nu = 0.002),
G = list(G1 = list(V = diag(2), nu = 2, alpha.mu = rep(0,2),
alpha.V = diag(25^2,2,2))))
m3 <- MCMCglmm(fixed = cbind(negative_affect, post_interest) ~ trait - 1 +
trait:gender_female +
trait:pre_interest,
random=~us(trait):participant_ID,
rcov=~corg(trait):units, # this denotes response vars varying across response rows
family = rep("gaussian",2),
data=as.data.frame(d_red),
prior=prior,
burnin = 3000,
nitt = 23000,
thin = 20,
verbose=TRUE)
##
## MCMC iteration = 0
##
## Acceptance ratio for liability set 1 = 0.000668
##
## MCMC iteration = 1000
##
## Acceptance ratio for liability set 1 = 0.323320
##
## MCMC iteration = 2000
##
## Acceptance ratio for liability set 1 = 0.315046
##
## MCMC iteration = 3000
##
## Acceptance ratio for liability set 1 = 0.306995
##
## MCMC iteration = 4000
##
## Acceptance ratio for liability set 1 = 0.330169
##
## MCMC iteration = 5000
##
## Acceptance ratio for liability set 1 = 0.332455
##
## MCMC iteration = 6000
##
## Acceptance ratio for liability set 1 = 0.332655
##
## MCMC iteration = 7000
##
## Acceptance ratio for liability set 1 = 0.331803
##
## MCMC iteration = 8000
##
## Acceptance ratio for liability set 1 = 0.332125
##
## MCMC iteration = 9000
##
## Acceptance ratio for liability set 1 = 0.331918
##
## MCMC iteration = 10000
##
## Acceptance ratio for liability set 1 = 0.331368
##
## MCMC iteration = 11000
##
## Acceptance ratio for liability set 1 = 0.330844
##
## MCMC iteration = 12000
##
## Acceptance ratio for liability set 1 = 0.331611
##
## MCMC iteration = 13000
##
## Acceptance ratio for liability set 1 = 0.331069
##
## MCMC iteration = 14000
##
## Acceptance ratio for liability set 1 = 0.331918
##
## MCMC iteration = 15000
##
## Acceptance ratio for liability set 1 = 0.330220
##
## MCMC iteration = 16000
##
## Acceptance ratio for liability set 1 = 0.330187
##
## MCMC iteration = 17000
##
## Acceptance ratio for liability set 1 = 0.330182
##
## MCMC iteration = 18000
##
## Acceptance ratio for liability set 1 = 0.330834
##
## MCMC iteration = 19000
##
## Acceptance ratio for liability set 1 = 0.331090
##
## MCMC iteration = 20000
##
## Acceptance ratio for liability set 1 = 0.331419
##
## MCMC iteration = 21000
##
## Acceptance ratio for liability set 1 = 0.330990
##
## MCMC iteration = 22000
##
## Acceptance ratio for liability set 1 = 0.330624
##
## MCMC iteration = 23000
##
## Acceptance ratio for liability set 1 = 0.331005
summary(m3)
m3_cor <- m3$VCV[,"traitpost_interest:traitnegative_affect.participant_ID"]/
(sqrt(m3$VCV[,"traitpost_interest:traitpost_interest.participant_ID"])*
sqrt(m3$VCV[,"traitnegative_affect:traitnegative_affect.participant_ID"]))
broom::tidy(m3)
posterior.mode(m3_cor)
HPDinterval(m3_cor)
##
## Iterations = 3001:22981
## Thinning interval = 20
## Sample size = 1000
##
## DIC: 11131.63
##
## G-structure: ~us(trait):participant_ID
##
## post.mean
## traitnegative_affect:traitnegative_affect.participant_ID 0.58350
## traitpost_interest:traitnegative_affect.participant_ID -0.04849
## traitnegative_affect:traitpost_interest.participant_ID -0.04849
## traitpost_interest:traitpost_interest.participant_ID 0.50841
## l-95% CI u-95% CI
## traitnegative_affect:traitnegative_affect.participant_ID 0.4465 0.72829
## traitpost_interest:traitnegative_affect.participant_ID -0.1564 0.04868
## traitnegative_affect:traitpost_interest.participant_ID -0.1564 0.04868
## traitpost_interest:traitpost_interest.participant_ID 0.3802 0.64832
## eff.samp
## traitnegative_affect:traitnegative_affect.participant_ID 1000.0
## traitpost_interest:traitnegative_affect.participant_ID 716.2
## traitnegative_affect:traitpost_interest.participant_ID 716.2
## traitpost_interest:traitpost_interest.participant_ID 823.6
##
## R-structure: ~corg(trait):units
##
## post.mean l-95% CI
## traitnegative_affect:traitnegative_affect.units 1.000000 1.00000
## traitpost_interest:traitnegative_affect.units -0.001201 -0.06502
## traitnegative_affect:traitpost_interest.units -0.001201 -0.06502
## traitpost_interest:traitpost_interest.units 1.000000 1.00000
## u-95% CI eff.samp
## traitnegative_affect:traitnegative_affect.units 1.00000 0
## traitpost_interest:traitnegative_affect.units 0.06424 1000
## traitnegative_affect:traitpost_interest.units 0.06424 1000
## traitpost_interest:traitpost_interest.units 1.00000 0
##
## Location effects: cbind(negative_affect, post_interest) ~ trait - 1 + trait:gender_female + trait:pre_interest
##
## post.mean l-95% CI u-95% CI eff.samp
## traitnegative_affect 1.97790 1.53342 2.42277 904.9
## traitpost_interest 1.30691 0.81554 1.84677 668.3
## traitnegative_affect:gender_female -0.13491 -0.37820 0.09063 1000.0
## traitpost_interest:gender_female -0.18277 -0.46270 0.05989 838.6
## traitnegative_affect:pre_interest -0.01275 -0.14531 0.11241 871.8
## traitpost_interest:pre_interest 0.61532 0.47221 0.76046 701.2
## pMCMC
## traitnegative_affect <0.001 ***
## traitpost_interest <0.001 ***
## traitnegative_affect:gender_female 0.248
## traitpost_interest:gender_female 0.146
## traitnegative_affect:pre_interest 0.888
## traitpost_interest:pre_interest <0.001 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## effect term estimate std.error
## 1 fixed traitnegative_affect 1.97790069 0.23409588
## 2 fixed traitpost_interest 1.30690663 0.25605558
## 3 fixed traitnegative_affect:gender_female -0.13491374 0.12101202
## 4 fixed traitpost_interest:gender_female -0.18276519 0.13188117
## 5 fixed traitnegative_affect:pre_interest -0.01275496 0.06951439
## 6 fixed traitpost_interest:pre_interest 0.61531586 0.07432096
## var1
## -0.03725614
## lower upper
## var1 -0.2870488 0.08172756
## attr(,"Probability")
## [1] 0.95