Loading data, setting up
library(tidyverse)
## ── Attaching packages ───────────────────────────────────── tidyverse 1.2.1 ──
## ✔ ggplot2 3.2.0 ✔ purrr 0.3.2
## ✔ tibble 2.1.3 ✔ dplyr 0.8.1
## ✔ tidyr 0.8.3 ✔ stringr 1.4.0
## ✔ readr 1.3.1 ✔ forcats 0.3.0
## ── Conflicts ──────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag() masks stats::lag()
library(lme4)
## Loading required package: Matrix
##
## Attaching package: 'Matrix'
## The following object is masked from 'package:tidyr':
##
## expand
library(brms)
## Loading required package: Rcpp
## Loading 'brms' package (version 2.8.0). Useful instructions
## can be found by typing help('brms'). A more detailed introduction
## to the package is available through vignette('brms_overview').
##
## Attaching package: 'brms'
## The following object is masked from 'package:lme4':
##
## ngrps
library(sjstats)
library(tidybayes)
##
## Attaching package: 'tidybayes'
## The following object is masked from 'package:sjstats':
##
## hdi
library(gt)
d1 <- read_csv("imuscle-esm.csv")
## Parsed with column specification:
## cols(
## .default = col_double(),
## Date = col_date(format = ""),
## beep_ID = col_character()
## )
## See spec(...) for full column specifications.
d2 <- read_csv("scimo-esm.csv")
## Parsed with column specification:
## cols(
## .default = col_double(),
## stud_ID = col_character(),
## response_date = col_date(format = ""),
## act_re = col_character(),
## ThirdQuar_2011 = col_logical(),
## participant_ID = col_character(),
## is_outlier = col_logical(),
## participant_ID1 = col_character(),
## beep_ID = col_character()
## )
## See spec(...) for full column specifications.
## Warning: 40 parsing failures.
## row col expected actual file
## 1727 ThirdQuar_2010 a double B 'scimo-esm.csv'
## 1727 FourthQuar_2010 a double B 'scimo-esm.csv'
## 1728 ThirdQuar_2010 a double B 'scimo-esm.csv'
## 1728 FourthQuar_2010 a double B 'scimo-esm.csv'
## 1729 ThirdQuar_2010 a double B 'scimo-esm.csv'
## .... ............... ........ ...... ...............
## See problems(...) for more details.
d3 <- read_csv("stemie-esm.csv")
## Parsed with column specification:
## cols(
## .default = col_double(),
## beep_ID = col_character(),
## additional_comments = col_character(),
## response_date = col_date(format = ""),
## time_after_signal_seconds = col_time(format = ""),
## sociedad_class = col_character(),
## beep_ID_new = col_character(),
## liked_best = col_character(),
## change_this_program = col_character(),
## why_recommend_to_friends = col_character(),
## what_will_study_comments = col_character(),
## job_30_years_old_comments = col_character(),
## race = col_character(),
## gender = col_character(),
## English = col_character(),
## Mathematics = col_character(),
## Science = col_character(),
## STAR_BOY_Math_ScaledScore = col_logical(),
## STAR_MOY_Math_ScaledScore = col_logical(),
## STAR_EOY_Math_ScaledScore = col_logical(),
## standardized_bos = col_character()
## # ... with 8 more columns
## )
## See spec(...) for full column specifications.
## Warning: 4324 parsing failures.
## row col expected actual file
## 1186 DaysEnrolled14/15 1/0/T/F/TRUE/FALSE 180 'stemie-esm.csv'
## 1186 English15/16 1/0/T/F/TRUE/FALSE B 'stemie-esm.csv'
## 1186 Math15/16 1/0/T/F/TRUE/FALSE B 'stemie-esm.csv'
## 1186 Science15/16 1/0/T/F/TRUE/FALSE B 'stemie-esm.csv'
## 1186 Standardized_Bos 15/16 1/0/T/F/TRUE/FALSE 726 'stemie-esm.csv'
## .... ...................... .................. ...... ................
## See problems(...) for more details.
# STEM-IE
d3 <- d3 %>%
mutate(program_ID = paste0(program_ID, " ", sociedad_class)) %>%
mutate(program_ID = str_replace(program_ID, " NA", ""))
# Overall
d1s <- d1 %>% dplyr::mutate(participant_ID = STUDID,
program_ID = str_c(d1$TEACHID, " - ", d1$PERIOD)) %>%
select(participant_ID,
beep_ID,
program_ID,
aff_enga,
beh_enga,
cog_enga) %>%
mutate(dataset = "imuscle",
participant_ID = as.character(participant_ID),
program_ID = as.character(program_ID))
d2s <- d2 %>%
select(participant_ID = stud_ID,
program_ID = teacher_ID,
beep_ID,
aff_enga,
beh_enga,
cog_enga) %>%
mutate(dataset = "scimo",
program_ID = as.character(program_ID))
subtact_one <- function(x) x - 1
d3s <- d3 %>%
mutate(program_ID = paste0(program_ID, " ", sociedad_class)) %>%
mutate(program_ID = str_replace(program_ID, " NA", "")) %>%
select(beep_ID = beep_ID_new,
cog_enga = cognitive_engagement,
beh_enga = behavioral_engagement,
aff_enga = affective_engagement,
participant_ID,
program_ID) %>%
mutate(dataset = "stemie") %>%
mutate(participant_ID = as.character(participant_ID)) %>%
mutate_at(vars(contains("enga")), subtact_one)
dd <- bind_rows(d1s, d2s, d3s)
dd <- dd %>%
group_by(participant_ID, program_ID, beep_ID) %>%
gather(key, val, -participant_ID, -beep_ID, -program_ID, -dataset) %>%
summarize(mean_enga = mean(val, na.rm = TRUE)) %>%
right_join(dd)
## Joining, by = c("participant_ID", "program_ID", "beep_ID")
dd
## # A tibble: 13,715 x 8
## # Groups: participant_ID, program_ID [1,315]
## participant_ID program_ID beep_ID mean_enga aff_enga beh_enga cog_enga
## <chr> <chr> <chr> <dbl> <dbl> <dbl> <dbl>
## 1 1160043 10 - 4 432012… 1.5 1 2 1.5
## 2 1170131 10 - 4 432012… 1.33 1 3 0
## 3 1170180 10 - 4 432012… 0.667 1 0 1
## 4 1170252 10 - 4 432012… 2.67 2 3 3
## 5 1170256 10 - 4 432012… 2 2 2 2
## 6 1170744 10 - 4 432012… 1.33 2 2 0
## 7 1170836 10 - 4 432012… 2.33 2 3 2
## 8 1170846 10 - 4 432012… 1.33 1 1.5 1.5
## 9 1170859 10 - 4 432012… 3 3 3 3
## 10 1170861 10 - 4 432012… 1.33 1 2.5 0.5
## # … with 13,705 more rows, and 1 more variable: dataset <chr>
Model
mood <- brm(mvbind(aff_enga, beh_enga, cog_enga) ~ -1 +
(1| a | participant_ID) + (1| b | beep_ID) + (1| c | program_ID) +
dataset,
data = dd,
iter = 2000, chains = 4,
cores = 3)
write_rds(mood, "final-model.rds")
mood <- read_rds("final-model.rds")
mood
## Family: MV(gaussian, gaussian, gaussian)
## Links: mu = identity; sigma = identity
## mu = identity; sigma = identity
## mu = identity; sigma = identity
## Formula: aff_enga ~ -1 + (1 | a | participant_ID) + (1 | b | beep_ID) + (1 | c | program_ID) + dataset
## beh_enga ~ -1 + (1 | a | participant_ID) + (1 | b | beep_ID) + (1 | c | program_ID) + dataset
## cog_enga ~ -1 + (1 | a | participant_ID) + (1 | b | beep_ID) + (1 | c | program_ID) + dataset
## Data: dd (Number of observations: 13556)
## Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
## total post-warmup samples = 4000
##
## Group-Level Effects:
## ~beep_ID (Number of levels: 1047)
## Estimate Est.Error l-95% CI
## sd(affenga_Intercept) 0.30 0.01 0.28
## sd(behenga_Intercept) 0.21 0.01 0.19
## sd(cogenga_Intercept) 0.18 0.01 0.16
## cor(affenga_Intercept,behenga_Intercept) 0.58 0.04 0.51
## cor(affenga_Intercept,cogenga_Intercept) 0.14 0.06 0.03
## cor(behenga_Intercept,cogenga_Intercept) 0.67 0.04 0.58
## u-95% CI Eff.Sample Rhat
## sd(affenga_Intercept) 0.31 1823 1.00
## sd(behenga_Intercept) 0.22 1618 1.00
## sd(cogenga_Intercept) 0.19 1974 1.00
## cor(affenga_Intercept,behenga_Intercept) 0.65 1534 1.00
## cor(affenga_Intercept,cogenga_Intercept) 0.25 1865 1.00
## cor(behenga_Intercept,cogenga_Intercept) 0.75 1416 1.00
##
## ~participant_ID (Number of levels: 1186)
## Estimate Est.Error l-95% CI
## sd(affenga_Intercept) 0.55 0.01 0.53
## sd(behenga_Intercept) 0.52 0.01 0.49
## sd(cogenga_Intercept) 0.67 0.02 0.64
## cor(affenga_Intercept,behenga_Intercept) 0.78 0.02 0.74
## cor(affenga_Intercept,cogenga_Intercept) 0.75 0.02 0.72
## cor(behenga_Intercept,cogenga_Intercept) 0.63 0.02 0.59
## u-95% CI Eff.Sample Rhat
## sd(affenga_Intercept) 0.58 825 1.00
## sd(behenga_Intercept) 0.54 766 1.00
## sd(cogenga_Intercept) 0.70 822 1.00
## cor(affenga_Intercept,behenga_Intercept) 0.80 988 1.00
## cor(affenga_Intercept,cogenga_Intercept) 0.78 848 1.00
## cor(behenga_Intercept,cogenga_Intercept) 0.67 833 1.00
##
## ~program_ID (Number of levels: 57)
## Estimate Est.Error l-95% CI
## sd(affenga_Intercept) 0.20 0.03 0.15
## sd(behenga_Intercept) 0.12 0.03 0.07
## sd(cogenga_Intercept) 0.22 0.03 0.16
## cor(affenga_Intercept,behenga_Intercept) 0.73 0.13 0.43
## cor(affenga_Intercept,cogenga_Intercept) 0.72 0.10 0.48
## cor(behenga_Intercept,cogenga_Intercept) 0.82 0.11 0.56
## u-95% CI Eff.Sample Rhat
## sd(affenga_Intercept) 0.27 838 1.00
## sd(behenga_Intercept) 0.17 467 1.00
## sd(cogenga_Intercept) 0.29 948 1.00
## cor(affenga_Intercept,behenga_Intercept) 0.92 935 1.00
## cor(affenga_Intercept,cogenga_Intercept) 0.88 830 1.01
## cor(behenga_Intercept,cogenga_Intercept) 0.98 491 1.00
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Eff.Sample
## affenga_datasetimuscle 1.69 0.05 1.60 1.78 1213
## affenga_datasetscimo 1.23 0.07 1.09 1.36 1195
## affenga_datasetstemie 1.87 0.07 1.72 2.01 925
## behenga_datasetimuscle 2.08 0.03 2.01 2.14 1361
## behenga_datasetscimo 1.80 0.05 1.70 1.89 1133
## behenga_datasetstemie 1.90 0.05 1.80 2.01 1006
## cogenga_datasetimuscle 1.53 0.05 1.44 1.63 1468
## cogenga_datasetscimo 1.12 0.07 0.97 1.27 1335
## cogenga_datasetstemie 1.64 0.08 1.49 1.80 1041
## Rhat
## affenga_datasetimuscle 1.00
## affenga_datasetscimo 1.00
## affenga_datasetstemie 1.00
## behenga_datasetimuscle 1.00
## behenga_datasetscimo 1.00
## behenga_datasetstemie 1.00
## cogenga_datasetimuscle 1.00
## cogenga_datasetscimo 1.00
## cogenga_datasetstemie 1.00
##
## Family Specific Parameters:
## Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
## sigma_affenga 0.65 0.00 0.64 0.66 5237 1.00
## sigma_behenga 0.66 0.00 0.65 0.66 5053 1.00
## sigma_cogenga 0.66 0.00 0.65 0.67 4731 1.00
##
## Residual Correlations:
## Estimate Est.Error l-95% CI u-95% CI Eff.Sample
## rescor(affenga,behenga) 0.48 0.01 0.47 0.50 5192
## rescor(affenga,cogenga) 0.42 0.01 0.40 0.43 4963
## rescor(behenga,cogenga) 0.34 0.01 0.33 0.36 4949
## Rhat
## rescor(affenga,behenga) 1.00
## rescor(affenga,cogenga) 1.00
## rescor(behenga,cogenga) 1.00
##
## Samples were drawn using sampling(NUTS). For each parameter, Eff.Sample
## is a crude measure of effective sample size, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
Summary table
d1 <- summary(mood)$fixed
d2 <- summary(mood)$random
d2a <- d2$beep_ID
d2b <- d2$participant_ID
d2c <- d2$program_ID
d3 <- summary(mood)$spec_pars
d4 <- summary(mood)$rescor_pars
d <- rbind(d1, d2a, d2b, d2c, d3, d4) %>%
as.data.frame() %>%
rownames_to_column("Var") %>%
as_tibble() %>%
select(Var, Estimate, Est_Error = Est.Error, Lower_CI = `l-95% CI`, Upper_CI = `u-95% CI`)
t <- d %>%
mutate(category = c(rep("Response-Dataset Intercept", 9),
rep("Situation Var", 3),
rep("Individual Var", 3),
rep("Program Var", 3),
rep("Ranef Corr.", 9),
rep("Residual", 3),
rep("Residual Corr.", 3))) %>%
group_by(category) %>%
gt() %>%
gt::fmt_number(columns = vars(Estimate, Est_Error, Lower_CI, Upper_CI)) %>%
gt::tab_row_group()
t
| Var |
Estimate |
Est_Error |
Lower_CI |
Upper_CI |
| Response-Dataset Intercept |
| affenga_datasetimuscle |
1.69 |
0.05 |
1.60 |
1.78 |
| affenga_datasetscimo |
1.23 |
0.07 |
1.09 |
1.36 |
| affenga_datasetstemie |
1.87 |
0.07 |
1.72 |
2.01 |
| behenga_datasetimuscle |
2.08 |
0.03 |
2.01 |
2.14 |
| behenga_datasetscimo |
1.80 |
0.05 |
1.70 |
1.89 |
| behenga_datasetstemie |
1.90 |
0.05 |
1.80 |
2.01 |
| cogenga_datasetimuscle |
1.53 |
0.05 |
1.44 |
1.63 |
| cogenga_datasetscimo |
1.12 |
0.07 |
0.97 |
1.27 |
| cogenga_datasetstemie |
1.64 |
0.08 |
1.49 |
1.80 |
| Situation Var |
| sd.affenga_Intercept. |
0.30 |
0.01 |
0.28 |
0.31 |
| sd.behenga_Intercept. |
0.21 |
0.01 |
0.19 |
0.22 |
| sd.cogenga_Intercept. |
0.18 |
0.01 |
0.16 |
0.19 |
| Individual Var |
| cor.affenga_Intercept.behenga_Intercept. |
0.58 |
0.04 |
0.51 |
0.65 |
| cor.affenga_Intercept.cogenga_Intercept. |
0.14 |
0.06 |
0.03 |
0.25 |
| cor.behenga_Intercept.cogenga_Intercept. |
0.67 |
0.04 |
0.58 |
0.75 |
| Program Var |
| sd.affenga_Intercept..1 |
0.55 |
0.01 |
0.53 |
0.58 |
| sd.behenga_Intercept..1 |
0.52 |
0.01 |
0.49 |
0.54 |
| sd.cogenga_Intercept..1 |
0.67 |
0.02 |
0.64 |
0.70 |
| Ranef Corr. |
| cor.affenga_Intercept.behenga_Intercept..1 |
0.78 |
0.02 |
0.74 |
0.80 |
| cor.affenga_Intercept.cogenga_Intercept..1 |
0.75 |
0.02 |
0.72 |
0.78 |
| cor.behenga_Intercept.cogenga_Intercept..1 |
0.63 |
0.02 |
0.59 |
0.67 |
| sd.affenga_Intercept..2 |
0.20 |
0.03 |
0.15 |
0.27 |
| sd.behenga_Intercept..2 |
0.12 |
0.03 |
0.07 |
0.17 |
| sd.cogenga_Intercept..2 |
0.22 |
0.03 |
0.16 |
0.29 |
| cor.affenga_Intercept.behenga_Intercept..2 |
0.73 |
0.13 |
0.43 |
0.92 |
| cor.affenga_Intercept.cogenga_Intercept..2 |
0.72 |
0.10 |
0.48 |
0.88 |
| cor.behenga_Intercept.cogenga_Intercept..2 |
0.82 |
0.11 |
0.56 |
0.98 |
| Residual |
| sigma_affenga |
0.65 |
0.00 |
0.64 |
0.66 |
| sigma_behenga |
0.66 |
0.00 |
0.65 |
0.66 |
| sigma_cogenga |
0.66 |
0.00 |
0.65 |
0.67 |
| Residual Corr. |
| rescor.affenga.behenga. |
0.48 |
0.01 |
0.47 |
0.50 |
| rescor.affenga.cogenga. |
0.42 |
0.01 |
0.40 |
0.43 |
| rescor.behenga.cogenga. |
0.34 |
0.01 |
0.33 |
0.36 |
t %>%
gt::gtsave("esm-icc-estimates.rtf")
# parnames(mood)[1:33]
# var component within outcome
## aff
hypothesis(mood,
c("sd_beep_ID__affenga_Intercept=sd_participant_ID__affenga_Intercept",
"sd_participant_ID__affenga_Intercept=sd_program_ID__affenga_Intercept",
"sd_beep_ID__affenga_Intercept=sd_program_ID__affenga_Intercept"),
class = NULL)
## Hypothesis Tests for class :
## Hypothesis Estimate Est.Error CI.Lower CI.Upper Evid.Ratio
## 1 (sd_beep_ID__affe... = 0 -0.26 0.02 -0.29 -0.23 NA
## 2 (sd_participant_I... = 0 0.35 0.03 0.28 0.41 NA
## 3 (sd_beep_ID__affe... = 0 0.09 0.03 0.03 0.15 NA
## Post.Prob Star
## 1 NA *
## 2 NA *
## 3 NA *
## ---
## '*': The expected value under the hypothesis lies outside the 95%-CI.
## Posterior probabilities of point hypotheses assume equal prior probabilities.
## beh
hypothesis(mood,
c("sd_beep_ID__behenga_Intercept=sd_participant_ID__behenga_Intercept",
"sd_participant_ID__behenga_Intercept=sd_program_ID__behenga_Intercept",
"sd_beep_ID__behenga_Intercept=sd_program_ID__behenga_Intercept"),
class = NULL)
## Hypothesis Tests for class :
## Hypothesis Estimate Est.Error CI.Lower CI.Upper Evid.Ratio
## 1 (sd_beep_ID__behe... = 0 -0.31 0.02 -0.34 -0.28 NA
## 2 (sd_participant_I... = 0 0.40 0.03 0.34 0.45 NA
## 3 (sd_beep_ID__behe... = 0 0.09 0.03 0.03 0.14 NA
## Post.Prob Star
## 1 NA *
## 2 NA *
## 3 NA *
## ---
## '*': The expected value under the hypothesis lies outside the 95%-CI.
## Posterior probabilities of point hypotheses assume equal prior probabilities.
## cog
hypothesis(mood,
c("sd_beep_ID__cogenga_Intercept=sd_participant_ID__cogenga_Intercept",
"sd_participant_ID__cogenga_Intercept=sd_program_ID__cogenga_Intercept",
"sd_beep_ID__cogenga_Intercept=sd_program_ID__cogenga_Intercept"),
class = NULL)
## Hypothesis Tests for class :
## Hypothesis Estimate Est.Error CI.Lower CI.Upper Evid.Ratio
## 1 (sd_beep_ID__coge... = 0 -0.49 0.02 -0.53 -0.46 NA
## 2 (sd_participant_I... = 0 0.45 0.04 0.37 0.52 NA
## 3 (sd_beep_ID__coge... = 0 -0.04 0.03 -0.11 0.02 NA
## Post.Prob Star
## 1 NA *
## 2 NA *
## 3 NA
## ---
## '*': The expected value under the hypothesis lies outside the 95%-CI.
## Posterior probabilities of point hypotheses assume equal prior probabilities.
# var component by outcome
## situation
hypothesis(mood,
c("sd_beep_ID__affenga_Intercept=sd_beep_ID__behenga_Intercept",
"sd_beep_ID__behenga_Intercept=sd_beep_ID__cogenga_Intercept",
"sd_beep_ID__affenga_Intercept=sd_beep_ID__cogenga_Intercept"),
class = NULL)
## Hypothesis Tests for class :
## Hypothesis Estimate Est.Error CI.Lower CI.Upper Evid.Ratio
## 1 (sd_beep_ID__affe... = 0 0.09 0.01 0.07 0.11 NA
## 2 (sd_beep_ID__behe... = 0 0.03 0.01 0.01 0.05 NA
## 3 (sd_beep_ID__affe... = 0 0.12 0.01 0.09 0.14 NA
## Post.Prob Star
## 1 NA *
## 2 NA *
## 3 NA *
## ---
## '*': The expected value under the hypothesis lies outside the 95%-CI.
## Posterior probabilities of point hypotheses assume equal prior probabilities.
## situation
hypothesis(mood,
c("sd_participant_ID__affenga_Intercept=sd_participant_ID__behenga_Intercept",
"sd_participant_ID__behenga_Intercept=sd_participant_ID__cogenga_Intercept",
"sd_participant_ID__affenga_Intercept=sd_participant_ID__cogenga_Intercept"),
class = NULL)
## Hypothesis Tests for class :
## Hypothesis Estimate Est.Error CI.Lower CI.Upper Evid.Ratio
## 1 (sd_participant_I... = 0 0.04 0.01 0.01 0.06 NA
## 2 (sd_participant_I... = 0 -0.15 0.02 -0.19 -0.12 NA
## 3 (sd_participant_I... = 0 -0.12 0.01 -0.15 -0.09 NA
## Post.Prob Star
## 1 NA *
## 2 NA *
## 3 NA *
## ---
## '*': The expected value under the hypothesis lies outside the 95%-CI.
## Posterior probabilities of point hypotheses assume equal prior probabilities.
## program
hypothesis(mood,
c("sd_program_ID__affenga_Intercept=sd_program_ID__behenga_Intercept",
"sd_program_ID__behenga_Intercept=sd_program_ID__cogenga_Intercept",
"sd_program_ID__affenga_Intercept=sd_program_ID__cogenga_Intercept"),
class = NULL)
## Hypothesis Tests for class :
## Hypothesis Estimate Est.Error CI.Lower CI.Upper Evid.Ratio
## 1 (sd_program_ID__a... = 0 0.09 0.03 0.03 0.14 NA
## 2 (sd_program_ID__b... = 0 -0.10 0.03 -0.16 -0.04 NA
## 3 (sd_program_ID__a... = 0 -0.01 0.03 -0.08 0.05 NA
## Post.Prob Star
## 1 NA *
## 2 NA *
## 3 NA
## ---
## '*': The expected value under the hypothesis lies outside the 95%-CI.
## Posterior probabilities of point hypotheses assume equal prior probabilities.
# Dataset
hypothesis(mood,
c("b_affenga_datasetstemie=b_affenga_datasetscimo",
"b_affenga_datasetstemie=b_affenga_datasetimuscle",
"b_affenga_datasetscimo=b_affenga_datasetimuscle"),
class = NULL)
## Hypothesis Tests for class :
## Hypothesis Estimate Est.Error CI.Lower CI.Upper Evid.Ratio
## 1 (b_affenga_datase... = 0 0.63 0.09 0.46 0.83 NA
## 2 (b_affenga_datase... = 0 0.18 0.08 0.01 0.34 NA
## 3 (b_affenga_datase... = 0 -0.46 0.08 -0.62 -0.29 NA
## Post.Prob Star
## 1 NA *
## 2 NA *
## 3 NA *
## ---
## '*': The expected value under the hypothesis lies outside the 95%-CI.
## Posterior probabilities of point hypotheses assume equal prior probabilities.
hypothesis(mood,
c("b_behenga_datasetstemie=b_behenga_datasetscimo",
"b_behenga_datasetstemie=b_behenga_datasetimuscle",
"b_behenga_datasetscimo=b_behenga_datasetimuscle"),
class = NULL)
## Hypothesis Tests for class :
## Hypothesis Estimate Est.Error CI.Lower CI.Upper Evid.Ratio
## 1 (b_behenga_datase... = 0 0.11 0.07 -0.03 0.25 NA
## 2 (b_behenga_datase... = 0 -0.17 0.06 -0.29 -0.05 NA
## 3 (b_behenga_datase... = 0 -0.28 0.06 -0.39 -0.16 NA
## Post.Prob Star
## 1 NA
## 2 NA *
## 3 NA *
## ---
## '*': The expected value under the hypothesis lies outside the 95%-CI.
## Posterior probabilities of point hypotheses assume equal prior probabilities.
hypothesis(mood,
c("b_cogenga_datasetstemie=b_cogenga_datasetscimo",
"b_cogenga_datasetstemie=b_cogenga_datasetimuscle",
"b_cogenga_datasetscimo=b_cogenga_datasetimuscle"),
class = NULL)
## Hypothesis Tests for class :
## Hypothesis Estimate Est.Error CI.Lower CI.Upper Evid.Ratio
## 1 (b_cogenga_datase... = 0 0.52 0.11 0.32 0.74 NA
## 2 (b_cogenga_datase... = 0 0.11 0.09 -0.07 0.29 NA
## 3 (b_cogenga_datase... = 0 -0.41 0.09 -0.59 -0.24 NA
## Post.Prob Star
## 1 NA *
## 2 NA
## 3 NA *
## ---
## '*': The expected value under the hypothesis lies outside the 95%-CI.
## Posterior probabilities of point hypotheses assume equal prior probabilities.
# Not 0
## aff
hypothesis(mood,
c("sd_beep_ID__affenga_Intercept>.02",
"sd_participant_ID__affenga_Intercept>0.02",
"sd_program_ID__affenga_Intercept>0.02"),
class = NULL)
## Hypothesis Tests for class :
## Hypothesis Estimate Est.Error CI.Lower CI.Upper Evid.Ratio
## 1 (sd_beep_ID__affe... > 0 0.28 0.01 0.26 Inf Inf
## 2 (sd_participant_I... > 0 0.53 0.01 0.51 Inf Inf
## 3 (sd_program_ID__a... > 0 0.18 0.03 0.14 Inf Inf
## Post.Prob Star
## 1 1 *
## 2 1 *
## 3 1 *
## ---
## '*': The expected value under the hypothesis lies outside the 95%-CI.
## Posterior probabilities of point hypotheses assume equal prior probabilities.
## beh
hypothesis(mood,
c("sd_beep_ID__behenga_Intercept>.02",
"sd_participant_ID__behenga_Intercept>0.02",
"sd_program_ID__behenga_Intercept>0.02"),
class = NULL)
## Hypothesis Tests for class :
## Hypothesis Estimate Est.Error CI.Lower CI.Upper Evid.Ratio
## 1 (sd_beep_ID__behe... > 0 0.19 0.01 0.17 Inf Inf
## 2 (sd_participant_I... > 0 0.50 0.01 0.47 Inf Inf
## 3 (sd_program_ID__b... > 0 0.10 0.03 0.06 Inf Inf
## Post.Prob Star
## 1 1 *
## 2 1 *
## 3 1 *
## ---
## '*': The expected value under the hypothesis lies outside the 95%-CI.
## Posterior probabilities of point hypotheses assume equal prior probabilities.
## cog
hypothesis(mood,
c("sd_beep_ID__cogenga_Intercept>.02",
"sd_participant_ID__cogenga_Intercept>0.02",
"sd_program_ID__cogenga_Intercept>0.02"),
class = NULL)
## Hypothesis Tests for class :
## Hypothesis Estimate Est.Error CI.Lower CI.Upper Evid.Ratio
## 1 (sd_beep_ID__coge... > 0 0.16 0.01 0.14 Inf Inf
## 2 (sd_participant_I... > 0 0.65 0.02 0.62 Inf Inf
## 3 (sd_program_ID__c... > 0 0.20 0.03 0.15 Inf Inf
## Post.Prob Star
## 1 1 *
## 2 1 *
## 3 1 *
## ---
## '*': The expected value under the hypothesis lies outside the 95%-CI.
## Posterior probabilities of point hypotheses assume equal prior probabilities.
## res correlations
## at program level
hypothesis(mood,
c("cor_program_ID__affenga_Intercept__behenga_Intercept=cor_program_ID__affenga_Intercept__cogenga_Intercept",
"cor_program_ID__behenga_Intercept__cogenga_Intercept=cor_program_ID__affenga_Intercept__cogenga_Intercept",
"cor_program_ID__affenga_Intercept__behenga_Intercept=cor_program_ID__behenga_Intercept__cogenga_Intercept"),
class = NULL)
## Hypothesis Tests for class :
## Hypothesis Estimate Est.Error CI.Lower CI.Upper Evid.Ratio
## 1 (cor_program_ID__... = 0 0.01 0.13 -0.27 0.26 NA
## 2 (cor_program_ID__... = 0 0.10 0.13 -0.14 0.35 NA
## 3 (cor_program_ID__... = 0 -0.10 0.14 -0.41 0.17 NA
## Post.Prob Star
## 1 NA
## 2 NA
## 3 NA
## ---
## '*': The expected value under the hypothesis lies outside the 95%-CI.
## Posterior probabilities of point hypotheses assume equal prior probabilities.
## at participant level
hypothesis(mood,
c("cor_participant_ID__affenga_Intercept__behenga_Intercept=cor_participant_ID__affenga_Intercept__cogenga_Intercept",
"cor_participant_ID__behenga_Intercept__cogenga_Intercept=cor_participant_ID__affenga_Intercept__cogenga_Intercept",
"cor_participant_ID__affenga_Intercept__behenga_Intercept=cor_participant_ID__behenga_Intercept__cogenga_Intercept"),
class = NULL)
## Hypothesis Tests for class :
## Hypothesis Estimate Est.Error CI.Lower CI.Upper Evid.Ratio
## 1 (cor_participant_... = 0 0.03 0.02 -0.01 0.06 NA
## 2 (cor_participant_... = 0 -0.12 0.02 -0.16 -0.08 NA
## 3 (cor_participant_... = 0 0.15 0.02 0.11 0.18 NA
## Post.Prob Star
## 1 NA
## 2 NA *
## 3 NA *
## ---
## '*': The expected value under the hypothesis lies outside the 95%-CI.
## Posterior probabilities of point hypotheses assume equal prior probabilities.
## aff and cog, at different levels
hypothesis(mood,
c("cor_participant_ID__affenga_Intercept__cogenga_Intercept=cor_beep_ID__affenga_Intercept__cogenga_Intercept",
"cor_participant_ID__affenga_Intercept__cogenga_Intercept=cor_program_ID__affenga_Intercept__cogenga_Intercept",
"cor_program_ID__affenga_Intercept__cogenga_Intercept=cor_beep_ID__affenga_Intercept__cogenga_Intercept"),
class = NULL)
## Hypothesis Tests for class :
## Hypothesis Estimate Est.Error CI.Lower CI.Upper Evid.Ratio
## 1 (cor_participant_... = 0 0.61 0.06 0.49 0.72 NA
## 2 (cor_participant_... = 0 0.03 0.10 -0.13 0.27 NA
## 3 (cor_program_ID__... = 0 0.58 0.12 0.32 0.78 NA
## Post.Prob Star
## 1 NA *
## 2 NA
## 3 NA *
## ---
## '*': The expected value under the hypothesis lies outside the 95%-CI.
## Posterior probabilities of point hypotheses assume equal prior probabilities.
Plot
sjstats::get_re_var(mood)
# AFF
resid <- mood %>%
spread_draws(sigma_affenga) %>%
mutate(est = sigma_affenga ^ 2) %>%
mutate(var = "residual") %>%
select(-.chain, -sigma_affenga) %>%
mutate(ICC = est /tot_var$total_var[1])
beeps <- mood%>%
spread_draws(r_beep_ID__affenga[condition, term]) %>%
group_by(.chain, .iteration) %>%
summarize(r_beep_ID = var(r_beep_ID__affenga)) %>%
ungroup() %>%
mutate(ICC = r_beep_ID/tot_var$total_var[1]) %>%
mutate(var = 'situation') %>%
rename(est = r_beep_ID)
peeps <- mood %>%
spread_draws(r_participant_ID__affenga[condition, term]) %>%
group_by(.chain, .iteration) %>%
summarize(r_participant_ID = var(r_participant_ID__affenga)) %>%
ungroup() %>%
mutate(ICC = r_participant_ID/tot_var$total_var[1]) %>%
mutate(var = 'individual') %>%
rename(est = r_participant_ID)
progs <- mood%>%
spread_draws(r_program_ID__affenga[condition, term]) %>%
group_by(.chain, .iteration) %>%
summarize(r_program_ID = var(r_program_ID__affenga)) %>%
ungroup() %>%
mutate(ICC = r_program_ID/tot_var$total_var[1]) %>%
mutate(var = 'program') %>%
rename(est = r_program_ID)
dda <- bind_rows(beeps, peeps, progs, resid)
gma <- dda %>%
group_by(var) %>%
summarize(ICC_mean = mean(ICC)) %>%
mutate(ICC = ICC_mean)
# BEH
resid <- mood %>%
spread_draws(sigma_behenga) %>%
mutate(est = sigma_behenga ^ 2) %>%
mutate(var = "residual") %>%
select(-.chain, -sigma_behenga) %>%
mutate(ICC = est /tot_var$total_var[1])
beeps <- mood%>%
spread_draws(r_beep_ID__behenga[condition, term]) %>%
group_by(.chain, .iteration) %>%
summarize(r_beep_ID = var(r_beep_ID__behenga)) %>%
ungroup() %>%
mutate(ICC = r_beep_ID/tot_var$total_var[1]) %>%
mutate(var = 'situation') %>%
rename(est = r_beep_ID)
peeps <- mood %>%
spread_draws(r_participant_ID__behenga[condition, term]) %>%
group_by(.chain, .iteration) %>%
summarize(r_participant_ID = var(r_participant_ID__behenga)) %>%
ungroup() %>%
mutate(ICC = r_participant_ID/tot_var$total_var[1]) %>%
mutate(var = 'individual') %>%
rename(est = r_participant_ID)
progs <- mood%>%
spread_draws(r_program_ID__behenga[condition, term]) %>%
group_by(.chain, .iteration) %>%
summarize(r_program_ID = var(r_program_ID__behenga)) %>%
ungroup() %>%
mutate(ICC = r_program_ID/tot_var$total_var[1]) %>%
mutate(var = 'program') %>%
rename(est = r_program_ID)
ddb <- bind_rows(beeps, peeps, progs, resid)
gmb <- ddb %>%
group_by(var) %>%
summarize(ICC_mean = mean(ICC)) %>%
mutate(ICC = ICC_mean)
# COG
resid <- mood %>%
spread_draws(sigma_cogenga) %>%
mutate(est = sigma_cogenga ^ 2) %>%
mutate(var = "residual") %>%
select(-.chain, -sigma_cogenga) %>%
mutate(ICC = est /tot_var$total_var[1])
beeps <- mood%>%
spread_draws(r_beep_ID__cogenga[condition, term]) %>%
group_by(.chain, .iteration) %>%
summarize(r_beep_ID = var(r_beep_ID__cogenga)) %>%
ungroup() %>%
mutate(ICC = r_beep_ID/tot_var$total_var[1]) %>%
mutate(var = 'situation') %>%
rename(est = r_beep_ID)
peeps <- mood %>%
spread_draws(r_participant_ID__cogenga[condition, term]) %>%
group_by(.chain, .iteration) %>%
summarize(r_participant_ID = var(r_participant_ID__cogenga)) %>%
ungroup() %>%
mutate(ICC = r_participant_ID/tot_var$total_var[1]) %>%
mutate(var = 'individual') %>%
rename(est = r_participant_ID)
progs <- mood%>%
spread_draws(r_program_ID__cogenga[condition, term]) %>%
group_by(.chain, .iteration) %>%
summarize(r_program_ID = var(r_program_ID__cogenga)) %>%
ungroup() %>%
mutate(ICC = r_program_ID/tot_var$total_var[1]) %>%
mutate(var = 'program') %>%
rename(est = r_program_ID)
ddc <- bind_rows(beeps, peeps, progs, resid)
gmc <- ddc %>%
group_by(var) %>%
summarize(ICC_mean = mean(ICC)) %>%
mutate(ICC = ICC_mean)
dda <- dda %>% mutate(out = "affective")
ddb <- ddb %>% mutate(out = "behavioral")
ddc <- ddc %>% mutate(out = "cognitive")
dddd <- dda %>%
bind_rows(ddb) %>%
bind_rows(ddc) %>%
select(-.draw)
gma <- gma %>% mutate(out = "affective")
gmb <- gmb %>% mutate(out = "behavioral")
gmc <- gmc %>% mutate(out = "cognitive")
gmd <- gma %>%
bind_rows(gmb) %>%
bind_rows(gmc)
dddd <- dddd %>% mutate(outcome = tools::toTitleCase(out))
gmd <- gmd %>% mutate(outcome = tools::toTitleCase(out),
var = tools::toTitleCase(var)) %>%
filter(var!="residual")
p <- dddd %>%
mutate(outcome = tools::toTitleCase(out),
var = tools::toTitleCase(var)) %>%
filter(var != "Residual") %>%
ggplot(aes(x = ICC, fill = var)) +
geom_density(alpha = .5) +
facet_grid(rows = vars(outcome)) +
geom_vline(aes(xintercept = ICC_mean, color = var), linetype = 4, size = .75, data = gmd) +
ggrepel::geom_label_repel(aes(
x = ICC_mean,
label = str_c(tools::toTitleCase(gmd$var), " ICC = ", round(gmd$ICC_mean, 3)),
y = 25),
data = gmd)
p +
# geom_density(aes(x = ICC), data = alldo) +
theme_bw() +
scale_fill_manual(values = c( '#fb8072', '#8dd3c7','#ffffb3')) +
scale_color_manual(values = c('#fb8072', '#8dd3c7', '#ffffb3')) +
ylab(NULL) +
xlab("Intra-class Correlation") +
theme(legend.position = "none") +
theme(axis.title.y=element_blank(),
axis.text.y=element_blank(),
axis.ticks.y=element_blank())
ggsave("2019-06-17-esm-iccs.png", width = 9, height = 9)