Dataset 13
1 Data Preparation
#Plot Affect
ggplot(daily, aes(na)) +
geom_histogram(binwidth = 0.1) +
ggtitle("Distribution of Negative Emotion Study 13") +
xlab("Negative Emotion") + ylab("Frequency") +
scale_x_continuous(breaks = seq(1, 5, 1)) +
scale_y_continuous(breaks = seq(0, 1000, 200)) +
theme(text = element_text(size = 10, family = "Times New Roman"),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
panel.background = element_blank(),
axis.line = element_line(colour = "black")) +
coord_cartesian(xlim = c(1, 5), ylim = c(0, 1000))## Warning: Removed 43 rows containing non-finite values (`stat_bin()`).
ggsave("NAS13.png", width = 4, height = 3)## Warning: Removed 43 rows containing non-finite values (`stat_bin()`).ggplot(daily, aes(pa)) +
geom_histogram(binwidth = 0.1) +
ggtitle("Distribution of Positive Emotion Study 13") +
xlab("Positive Emotion") + ylab("Frequency") +
scale_x_continuous(breaks = seq(1, 5, 1)) +
scale_y_continuous(breaks = seq(0, 1500, 200)) +
theme(text = element_text(size = 10, family = "Times New Roman"),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
panel.background = element_blank(),
axis.line = element_line(colour = "black")) +
coord_cartesian(xlim = c(1, 5), ylim = c(0, 1200))## Warning: Removed 43 rows containing non-finite values (`stat_bin()`).
ggsave("PAS13.png", width = 4, height = 3)## Warning: Removed 43 rows containing non-finite values (`stat_bin()`).1.1 Censoring
daily$Acens <- case_when(daily$na == 1 ~ "left",
daily$na == 5 ~ "right",
TRUE ~ "none")
table(daily$Acens)##
## left none right
## 577 5527 2daily$Acensp <- case_when(daily$pa == 1 ~ "left",
daily$pa == 5 ~ "right",
TRUE ~ "none")
table(daily$Acensp)##
## left none right
## 50 6025 312 BCLSM Negative Emotion
na_model_neuro3 <- brm(bf(na | cens(Acens) ~ neurot + (1|session),
sigma ~ neurot + (1|session)), data = daily,
iter = 7000, warmup = 2000,chains = 8, init = 0.1,
control = list(adapt_delta = .99),
file = "na_model_neuro3")## Warning: Rows containing NAs were excluded from the model.print(na_model_neuro3)## Family: gaussian
## Links: mu = identity; sigma = log
## Formula: na | cens(Acens) ~ neurot + (1 | session)
## sigma ~ neurot + (1 | session)
## Data: daily (Number of observations: 6063)
## Draws: 8 chains, each with iter = 7000; warmup = 2000; thin = 1;
## total post-warmup draws = 40000
##
## Group-Level Effects:
## ~session (Number of levels: 444)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept) 0.43 0.02 0.40 0.46 1.00 5197 10825
## sd(sigma_Intercept) 0.34 0.02 0.31 0.38 1.00 15232 25577
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept 0.65 0.08 0.50 0.80 1.00 2772 6056
## sigma_Intercept -1.26 0.07 -1.40 -1.12 1.00 12628 21519
## neurot 0.41 0.02 0.37 0.46 1.00 3083 6494
## sigma_neurot 0.11 0.02 0.07 0.16 1.00 12370 21374
##
## Draws were sampled using sample(hmc). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).plot(na_model_neuro3)
prior_summary(na_model_neuro3)## prior class coef group resp dpar nlpar lb ub
## (flat) b
## (flat) b neurot
## (flat) b sigma
## (flat) b neurot sigma
## student_t(3, 1.8, 2.5) Intercept
## student_t(3, 0, 2.5) Intercept sigma
## student_t(3, 0, 2.5) sd 0
## student_t(3, 0, 2.5) sd sigma 0
## student_t(3, 0, 2.5) sd session 0
## student_t(3, 0, 2.5) sd Intercept session 0
## student_t(3, 0, 2.5) sd session sigma 0
## student_t(3, 0, 2.5) sd Intercept session sigma 0
## source
## default
## (vectorized)
## default
## (vectorized)
## default
## default
## default
## default
## (vectorized)
## (vectorized)
## (vectorized)
## (vectorized)2.1 Model comparison
2.1.1 scale vs. no scale parameter
na_model_neuro2 <- brm(na | cens(Acens) ~ neurot + (1|session), data = daily,
iter = 3000, warmup = 1000, chains = 8,
file = "na_model_neuro2")## Warning: Rows containing NAs were excluded from the model.print(na_model_neuro2)## Family: gaussian
## Links: mu = identity; sigma = identity
## Formula: na | cens(Acens) ~ neurot + (1 | session)
## Data: daily (Number of observations: 6063)
## Draws: 8 chains, each with iter = 3000; warmup = 1000; thin = 1;
## total post-warmup draws = 16000
##
## Group-Level Effects:
## ~session (Number of levels: 444)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept) 0.44 0.02 0.41 0.47 1.00 1923 3980
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept 0.64 0.08 0.49 0.79 1.01 1221 2457
## neurot 0.42 0.03 0.37 0.47 1.01 1239 2484
##
## Family Specific Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma 0.44 0.00 0.43 0.45 1.00 20507 11459
##
## Draws were sampled using sample(hmc). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).plot(na_model_neuro2)
modelA <- na_model_neuro2
modelB <- na_model_neuro3
modelA <- add_criterion(modelA, "loo")
modelB <- add_criterion(modelB, "loo")
loo <- loo_compare(modelA,modelB, criterion = "loo")
loo <- as.data.frame(loo)
loo$Dataset <- "13"
loo <- tibble::rownames_to_column(loo, "model")
library("writexl")
write_xlsx(loo,"~/looDataset13.xlsx")
kable(loo)| model | elpd_diff | se_diff | elpd_loo | se_elpd_loo | p_loo | se_p_loo | looic | se_looic | Dataset |
|---|---|---|---|---|---|---|---|---|---|
| modelB | 0.0000 | 0.0000 | -3511.970 | 67.36137 | 717.2413 | 17.748407 | 7023.941 | 134.7227 | 13 |
| modelA | -423.5521 | 35.9603 | -3935.523 | 69.29213 | 398.8505 | 8.936909 | 7871.045 | 138.5843 | 13 |
2.1.2 censoring vs. no censoring
na_model_neuro4 <- brm(bf(na ~ neurot + (1|session),
sigma ~ neurot + (1|session)), data = daily,
iter = 7000, warmup = 2000,chains = 4, init = 0.1,
control = list(adapt_delta = .99),
file = "na_model_neuro4")## Warning: Rows containing NAs were excluded from the model.print(na_model_neuro4)## Family: gaussian
## Links: mu = identity; sigma = log
## Formula: na ~ neurot + (1 | session)
## sigma ~ neurot + (1 | session)
## Data: daily (Number of observations: 6063)
## Draws: 4 chains, each with iter = 7000; warmup = 2000; thin = 1;
## total post-warmup draws = 20000
##
## Group-Level Effects:
## ~session (Number of levels: 444)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept) 0.40 0.01 0.37 0.43 1.00 1990 4231
## sd(sigma_Intercept) 0.37 0.02 0.34 0.41 1.00 6589 11188
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept 0.77 0.07 0.63 0.91 1.00 1121 2615
## sigma_Intercept -1.70 0.07 -1.85 -1.56 1.00 6311 10516
## neurot 0.38 0.02 0.34 0.43 1.00 1338 3009
## sigma_neurot 0.23 0.02 0.18 0.28 1.00 6337 11139
##
## Draws were sampled using sample(hmc). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).extract_param <- function(model, parameter) {
ci <- posterior_summary(model, variable = parameter)
est <- sprintf("%.2f %.2f [%.2f;%.2f]", ci[,"Estimate"],ci[,"Est.Error"], ci[,"Q2.5"], ci[,"Q97.5"])
est
}
results_LScens <- data.frame(matrix(nrow = 2, # Modelle & RVI
ncol = 6+1)) # 2 Affekte a 3 Spalten
names(results_LScens) <- c("model", "negemo_b_neuro", "negemo_b_neuro_sigma", "negemo_sigma",
"posemo_b_neuro", "posemo_b_neuro_sigma", "posemo_sigma"
)
results_LScens$model <- c("modelCensoring", "modelnoCensoring")
#NA
results_LScens[1, "negemo_b_neuro"] <- extract_param(na_model_neuro3, "b_neurot")
results_LScens[1, "negemo_b_neuro_sigma"] <- extract_param(na_model_neuro3, "b_sigma_neurot")
results_LScens[1, "negemo_sigma"] <- extract_param(na_model_neuro3, "b_sigma_Intercept")
results_LScens[2, "negemo_b_neuro"] <- extract_param(na_model_neuro4, "b_neurot")
results_LScens[2, "negemo_b_neuro_sigma"] <- extract_param(na_model_neuro4, "b_sigma_neurot")
results_LScens[2, "negemo_sigma"] <- extract_param(na_model_neuro4, "b_sigma_Intercept")2.1.3 BCLSM vs. model C (two-part model)
daily <- daily %>% left_join(daily %>% distinct(session, neurot) %>% mutate(neuro_Q =Hmisc::cut2(neurot, g = 4)), by = c("session", "neurot"))
na_model_neuro_jinxed <- brm(bf(na | cens(Acens) ~ neurot + (1|gr(session, by = neuro_Q)),
sigma ~ neurot + (1|session)), data = daily,
iter = 5000, warmup = 2000, chains = 4,
control = list(adapt_delta = .95), init = 0.1,
file = "na_model_neuro3_jinxed")## Warning: Rows containing NAs were excluded from the model.print(na_model_neuro_jinxed)## Family: gaussian
## Links: mu = identity; sigma = log
## Formula: na | cens(Acens) ~ neurot + (1 | gr(session, by = neuro_Q))
## sigma ~ neurot + (1 | session)
## Data: daily (Number of observations: 6063)
## Draws: 4 chains, each with iter = 5000; warmup = 2000; thin = 1;
## total post-warmup draws = 12000
##
## Group-Level Effects:
## ~session (Number of levels: 444)
## Estimate Est.Error l-95% CI u-95% CI Rhat
## sd(Intercept:neuro_Q[1.00,2.50)) 0.38 0.02 0.34 0.44 1.00
## sd(Intercept:neuro_Q[2.50,3.00)) 0.46 0.04 0.39 0.55 1.00
## sd(Intercept:neuro_Q[3.00,3.67)) 0.41 0.03 0.36 0.47 1.00
## sd(Intercept:neuro_Q[3.67,4.83]) 0.52 0.04 0.44 0.60 1.00
## sd(sigma_Intercept) 0.34 0.02 0.31 0.38 1.00
## Bulk_ESS Tail_ESS
## sd(Intercept:neuro_Q[1.00,2.50)) 1708 3294
## sd(Intercept:neuro_Q[2.50,3.00)) 1357 2967
## sd(Intercept:neuro_Q[3.00,3.67)) 1256 2831
## sd(Intercept:neuro_Q[3.67,4.83]) 1386 2344
## sd(sigma_Intercept) 3769 6329
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept 0.65 0.07 0.50 0.80 1.01 609 1142
## sigma_Intercept -1.26 0.07 -1.40 -1.12 1.00 3291 5767
## neurot 0.42 0.03 0.36 0.47 1.00 625 1308
## sigma_neurot 0.11 0.02 0.07 0.16 1.00 3306 5693
##
## Draws were sampled using sample(hmc). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).modelA <- na_model_neuro3
modelB <- na_model_neuro_jinxed
modelA <- add_criterion(modelA, "loo")
modelB <- add_criterion(modelB, "loo")
loo_c <- loo_compare(modelA,modelB, criterion = "loo")
loo_c <- as.data.frame(loo)
loo_c$Dataset <- "13"
#loo_c <- tibble::rownames_to_column(loo, "model")
library("writexl")
write_xlsx(loo_c,"~/Study13loojinxed.xlsx")
kable(loo_c)| model | elpd_diff | se_diff | elpd_loo | se_elpd_loo | p_loo | se_p_loo | looic | se_looic | Dataset |
|---|---|---|---|---|---|---|---|---|---|
| modelB | 0.0000 | 0.0000 | -3511.970 | 67.36137 | 717.2413 | 17.748407 | 7023.941 | 134.7227 | 13 |
| modelA | -423.5521 | 35.9603 | -3935.523 | 69.29213 | 398.8505 | 8.936909 | 7871.045 | 138.5843 | 13 |
2.2 control for gender
table(Start$sex)##
## 1 2 3
## 383 189 5#1=female, 2=male, 3=diverse
gender <- select(Start, session, sex)
# bins sex to daily
daily2 <- merge(daily, gender, all.x = T, by= "session")
daily2$sex <- as.numeric(daily2$sex)
## not enough diverse people to compare them (excluded from model)
daily2$sex[daily2$sex == 3] <- NA
daily2$sex[daily2$sex == 2] <- 0
daily2$sex <- as.factor(daily2$sex)
table(daily2$sex, exclude = NULL)##
## 0 1 <NA>
## 2080 3969 57m3n_gender <- brm(bf(na | cens(Acens) ~ neurot+ sex + (1|session),
sigma ~ neurot+sex + (1|session)), data = daily2,
iter = 5000, warmup = 2000,chains = 4,
control = list(adapt_delta = .95), inits = 0.1,
file = "m3n_gender")## Warning: Argument 'inits' is deprecated. Please use argument 'init' instead.## Warning: Rows containing NAs were excluded from the model.print(m3n_gender)## Family: gaussian
## Links: mu = identity; sigma = log
## Formula: na | cens(Acens) ~ neurot + sex + (1 | session)
## sigma ~ neurot + sex + (1 | session)
## Data: daily2 (Number of observations: 6006)
## Draws: 4 chains, each with iter = 5000; warmup = 2000; thin = 1;
## total post-warmup draws = 12000
##
## Group-Level Effects:
## ~session (Number of levels: 440)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept) 0.43 0.02 0.40 0.46 1.00 1246 3193
## sd(sigma_Intercept) 0.34 0.02 0.31 0.38 1.00 3762 6305
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept 0.65 0.08 0.50 0.81 1.01 563 1129
## sigma_Intercept -1.29 0.07 -1.43 -1.15 1.00 3538 6321
## neurot 0.41 0.03 0.36 0.46 1.01 550 1356
## sex1 0.01 0.05 -0.08 0.10 1.01 607 1275
## sigma_neurot 0.10 0.02 0.05 0.15 1.00 3257 5750
## sigma_sex1 0.08 0.04 -0.00 0.17 1.00 3027 4645
##
## Draws were sampled using sample(hmc). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).3 BCLSM Positive Emotion
pa_model_neuro3 <- brm(bf(pa | cens(Acensp) ~ neurot + (1|session),
sigma ~ neurot + (1|session)), data = daily,
control = list(adapt_delta = .99),chains = 8,
iter = 5000, warmup = 2000,
file = "pa_model_neuro3")## Warning: Rows containing NAs were excluded from the model.print(pa_model_neuro3)## Family: gaussian
## Links: mu = identity; sigma = log
## Formula: pa | cens(Acensp) ~ neurot + (1 | session)
## sigma ~ neurot + (1 | session)
## Data: daily (Number of observations: 6063)
## Draws: 8 chains, each with iter = 5000; warmup = 2000; thin = 1;
## total post-warmup draws = 24000
##
## Group-Level Effects:
## ~session (Number of levels: 444)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept) 0.52 0.02 0.49 0.56 1.00 3063 6621
## sd(sigma_Intercept) 0.37 0.02 0.33 0.43 1.00 4278 7931
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept 4.25 0.09 4.07 4.42 1.00 1755 4122
## sigma_Intercept -1.34 0.07 -1.48 -1.19 1.00 8360 13296
## neurot -0.36 0.03 -0.41 -0.30 1.00 1910 4139
## sigma_neurot 0.17 0.02 0.12 0.22 1.00 8582 13387
##
## Draws were sampled using sample(hmc). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).plot(pa_model_neuro3)
prior_summary(pa_model_neuro3)## prior class coef group resp dpar nlpar lb ub
## (flat) b
## (flat) b neurot
## (flat) b sigma
## (flat) b neurot sigma
## student_t(3, 3.3, 2.5) Intercept
## student_t(3, 0, 2.5) Intercept sigma
## student_t(3, 0, 2.5) sd 0
## student_t(3, 0, 2.5) sd sigma 0
## student_t(3, 0, 2.5) sd session 0
## student_t(3, 0, 2.5) sd Intercept session 0
## student_t(3, 0, 2.5) sd session sigma 0
## student_t(3, 0, 2.5) sd Intercept session sigma 0
## source
## default
## (vectorized)
## default
## (vectorized)
## default
## default
## default
## default
## (vectorized)
## (vectorized)
## (vectorized)
## (vectorized)pp_check(pa_model_neuro3)## Using 10 posterior draws for ppc type 'dens_overlay' by default.## Warning: Censored responses are not shown in 'pp_check'.
3.1 Model comparison
3.1.1 scale vs. no scale parameter
pa_model_neuro2 <- brm(pa | cens(Acensp) ~ neurot + (1|session), data = daily,
iter = 4000, warmup = 1000, chains = 8,
file = "pa_model_neuro2")## Warning: Rows containing NAs were excluded from the model.print(pa_model_neuro2)## Family: gaussian
## Links: mu = identity; sigma = identity
## Formula: pa | cens(Acensp) ~ neurot + (1 | session)
## Data: daily (Number of observations: 6063)
## Draws: 8 chains, each with iter = 4000; warmup = 1000; thin = 1;
## total post-warmup draws = 24000
##
## Group-Level Effects:
## ~session (Number of levels: 444)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept) 0.52 0.02 0.49 0.56 1.00 2605 5220
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept 4.24 0.09 4.06 4.42 1.00 1507 2977
## neurot -0.35 0.03 -0.41 -0.29 1.00 1540 3126
##
## Family Specific Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma 0.47 0.00 0.47 0.48 1.00 28578 18787
##
## Draws were sampled using sample(hmc). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).plot(pa_model_neuro2)
modelAp <- pa_model_neuro2
modelBp <- pa_model_neuro3
modelAp <- add_criterion(modelAp, "loo")
modelBp <- add_criterion(modelBp, "loo")
looP <- loo_compare(modelAp,modelBp, criterion = "loo")
looP <- as.data.frame(looP)
looP$Dataset <- "13"
#looP <- tibble::rownames_to_column(looP, "model")
library("writexl")
write_xlsx(looP,"~/looPDataset13.xlsx")
kable(looP)| elpd_diff | se_diff | elpd_loo | se_elpd_loo | p_loo | se_p_loo | looic | se_looic | Dataset | |
|---|---|---|---|---|---|---|---|---|---|
| modelBp | 0.000 | 0.00000 | -3880.188 | 68.97765 | 748.2984 | 19.478463 | 7760.377 | 137.9553 | 13 |
| modelAp | -453.445 | 39.04761 | -4333.633 | 66.75055 | 404.8004 | 8.476593 | 8667.267 | 133.5011 | 13 |
3.1.2 censoring vs. no censoring
pa_model_neuro4 <- brm(bf(pa ~ neurot + (1|session),
sigma ~ neurot + (1|session)), data = daily,
control = list(adapt_delta = .99),chains = 4,
iter = 5000, warmup = 2000,
file = "pa_model_neuro4")## Warning: Rows containing NAs were excluded from the model.print(pa_model_neuro4)## Family: gaussian
## Links: mu = identity; sigma = log
## Formula: pa ~ neurot + (1 | session)
## sigma ~ neurot + (1 | session)
## Data: daily (Number of observations: 6063)
## Draws: 4 chains, each with iter = 5000; warmup = 2000; thin = 1;
## total post-warmup draws = 12000
##
## Group-Level Effects:
## ~session (Number of levels: 444)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept) 0.51 0.02 0.48 0.55 1.01 977 2107
## sd(sigma_Intercept) 0.37 0.02 0.32 0.42 1.00 2332 4516
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept 4.24 0.09 4.06 4.41 1.01 541 1019
## sigma_Intercept -1.33 0.07 -1.48 -1.19 1.00 2526 4900
## neurot -0.35 0.03 -0.41 -0.29 1.01 590 1280
## sigma_neurot 0.16 0.02 0.11 0.21 1.00 2640 5095
##
## Draws were sampled using sample(hmc). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).results_LScens[1, "posemo_b_neuro"] <- extract_param(pa_model_neuro3, "b_neurot")
results_LScens[1, "posemo_b_neuro_sigma"] <- extract_param(pa_model_neuro3, "b_sigma_neurot")
results_LScens[1, "posemo_sigma"] <- extract_param(pa_model_neuro3, "b_sigma_Intercept")
results_LScens[2, "posemo_b_neuro"] <- extract_param(pa_model_neuro4, "b_neurot")
results_LScens[2, "posemo_b_neuro_sigma"] <- extract_param(pa_model_neuro4, "b_sigma_neurot")
results_LScens[2, "posemo_sigma"] <- extract_param(pa_model_neuro4, "b_sigma_Intercept")
library("writexl")
write_xlsx(results_LScens,"~/results_LSCENS.xlsx")3.1.3 BCLSM vs. model C (two-part model)
pa_model_neuro_jinxed <- brm(bf(pa | cens(Acens) ~ neurot + (1|gr(session, by = neuro_Q)),
sigma ~ neurot + (1|session)), data = daily,
iter = 5000, warmup = 2000, chains = 4,
control = list(adapt_delta = .95), init = 0.1,
file = "pa_model_neuro3_jinxed")## Warning: Rows containing NAs were excluded from the model.print(pa_model_neuro_jinxed)## Family: gaussian
## Links: mu = identity; sigma = log
## Formula: pa | cens(Acens) ~ neurot + (1 | gr(session, by = neuro_Q))
## sigma ~ neurot + (1 | session)
## Data: daily (Number of observations: 6063)
## Draws: 4 chains, each with iter = 5000; warmup = 2000; thin = 1;
## total post-warmup draws = 12000
##
## Group-Level Effects:
## ~session (Number of levels: 444)
## Estimate Est.Error l-95% CI u-95% CI Rhat
## sd(Intercept:neuro_Q[1.00,2.50)) 0.51 0.03 0.45 0.58 1.00
## sd(Intercept:neuro_Q[2.50,3.00)) 0.46 0.04 0.39 0.55 1.00
## sd(Intercept:neuro_Q[3.00,3.67)) 0.49 0.03 0.43 0.55 1.00
## sd(Intercept:neuro_Q[3.67,4.83]) 0.57 0.05 0.49 0.67 1.00
## sd(sigma_Intercept) 0.35 0.02 0.31 0.39 1.00
## Bulk_ESS Tail_ESS
## sd(Intercept:neuro_Q[1.00,2.50)) 1178 2596
## sd(Intercept:neuro_Q[2.50,3.00)) 1848 3974
## sd(Intercept:neuro_Q[3.00,3.67)) 1727 3750
## sd(Intercept:neuro_Q[3.67,4.83]) 1827 3942
## sd(sigma_Intercept) 2762 5100
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept 4.02 0.09 3.84 4.21 1.01 671 1238
## sigma_Intercept -1.27 0.07 -1.42 -1.13 1.00 4525 6857
## neurot -0.30 0.03 -0.36 -0.23 1.01 684 1326
## sigma_neurot 0.14 0.02 0.10 0.19 1.00 4533 6528
##
## Draws were sampled using sample(hmc). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).modelA <- pa_model_neuro3
modelB <- pa_model_neuro_jinxed
modelA <- add_criterion(modelA, "loo")
modelB <- add_criterion(modelB, "loo")
loo_cP <- loo_compare(modelA,modelB, criterion = "loo")## Warning: Not all models have the same y variable. ('yhash' attributes do not
## match)loo_cP <- as.data.frame(loo)
loo_cP$Dataset <- "13"
#loo_cP <- tibble::rownames_to_column(loo_cP, "model")
library("writexl")
write_xlsx(loo_cP,"~/Study13POloojinxed.xlsx")
kable(loo_cP)| model | elpd_diff | se_diff | elpd_loo | se_elpd_loo | p_loo | se_p_loo | looic | se_looic | Dataset |
|---|---|---|---|---|---|---|---|---|---|
| modelB | 0.0000 | 0.0000 | -3511.970 | 67.36137 | 717.2413 | 17.748407 | 7023.941 | 134.7227 | 13 |
| modelA | -423.5521 | 35.9603 | -3935.523 | 69.29213 | 398.8505 | 8.936909 | 7871.045 | 138.5843 | 13 |
extract_param <- function(model, parameter) {
ci <- posterior_summary(model, variable = parameter)
est <- sprintf("%.2f %.2f [%.2f;%.2f]", ci[,"Estimate"],ci[,"Est.Error"], ci[,"Q2.5"], ci[,"Q97.5"])
est
}
results_LS <- data.frame(matrix(nrow = 7, # Modelle & RVI
ncol = 8+1)) # 2 Affekte a 3 Spalten
names(results_LS) <-c("model", "negemo_b_neuro", "negemo_b_neuro_sigma", "negemo_sigma", "b_neg_sigma_sex",
"posemo_b_neuro", "posemo_b_neuro_sigma", "posemo_sigma", "b_pos_sigma_sex"
)
results_LS$model <- c("model1", "model2", "model3",
"RSD", "RSD_weight", "SD", "gender")
#NA
results_LS[2, "negemo_b_neuro"] <- extract_param(na_model_neuro2, "b_neurot")
results_LS[2, "negemo_sigma"] <- extract_param(na_model_neuro2, "sigma")
results_LS[3, "negemo_b_neuro"] <- extract_param(na_model_neuro3, "b_neurot")
results_LS[3, "negemo_b_neuro_sigma"] <- extract_param(na_model_neuro3, "b_sigma_neurot")
results_LS[3, "negemo_sigma"] <- extract_param(na_model_neuro3, "b_sigma_Intercept")
#pa
results_LS[2, "posemo_b_neuro"] <- extract_param(pa_model_neuro2, "b_neurot")
results_LS[2, "posemo_sigma"] <- extract_param(pa_model_neuro2, "sigma")
results_LS[3, "posemo_b_neuro"] <- extract_param(pa_model_neuro3, "b_neurot")
results_LS[3, "posemo_b_neuro_sigma"] <- extract_param(pa_model_neuro3, "b_sigma_neurot")
results_LS[3, "posemo_sigma"] <- extract_param(pa_model_neuro3, "b_sigma_Intercept")
#gender
results_LS[7, "negemo_b_neuro"] <- extract_param(m3n_gender, "b_neurot")
results_LS[7, "negemo_b_neuro_sigma"] <- extract_param(m3n_gender, "b_sigma_neurot")
results_LS[7, "negemo_sigma"] <- extract_param(m3n_gender, "b_sigma_Intercept")
results_LS[7, "b_neg_sigma_sex"] <- extract_param(m3n_gender, "b_sigma_sex1")4 RVI (Relative-Variability-Index)
4.1 Unweighted RVI
#id <- unique(daily$session)
#id <- as.data.frame(id)
Start$RSD_NA <- NA
for (i in 1:nrow(Start)) {
# aktuelle Zeilennummer
if (Start$session[i] %in% Daily$session) {
if (
sum(!is.na(Daily$na[Daily$session == Start$session[i]])) >= 5
) {
Start$RSD_NA[i] <- relativeSD(Daily$na[Daily$session == Start$session[i]],
1, 5)
}
}
}
Start$logrsd_n <- log(Start$RSD_NA)
m_rsd_na <- brm(logrsd_n ~ Neurot, data= Start,
file = "na_logrsd_uw")## Warning: Rows containing NAs were excluded from the model.print(m_rsd_na)## Family: gaussian
## Links: mu = identity; sigma = identity
## Formula: logrsd_n ~ Neurot
## Data: Start (Number of observations: 441)
## Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
## total post-warmup draws = 4000
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept -1.26 0.07 -1.39 -1.13 1.00 4553 3078
## Neurot -0.05 0.02 -0.10 -0.01 1.00 4551 2971
##
## Family Specific Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma 0.40 0.01 0.38 0.43 1.00 3944 3089
##
## Draws were sampled using sample(hmc). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).results_LS[4,3] <- extract_param(m_rsd_na, "b_Neurot")
Start$RSD_PA <- NA
for (i in 1:nrow(Start)) {
# aktuelle Zeilennummer
if (Start$session[i] %in% Daily$session) {
if (
sum(!is.na(Daily$pa[Daily$session == Start$session[i]])) >= 5
) {
Start$RSD_PA[i] <- relativeSD(Daily$pa[Daily$session == Start$session[i]],
1, 5)
}
}
}
Start$RSD_PA[Start$RSD_PA == 0] <- NA # sonst log(0) = - Inf
Start$logrsd_p <- log(Start$RSD_PA)
range(Start$logrsd_p, na.rm = T)## [1] -3.144725e+00 -1.110223e-16m_rsd_pa <- brm(logrsd_p ~ Neurot, data= Start,
file = "pa_logrsd_uw")## Warning: Rows containing NAs were excluded from the model.print(m_rsd_pa)## Family: gaussian
## Links: mu = identity; sigma = identity
## Formula: logrsd_p ~ Neurot
## Data: Start (Number of observations: 437)
## Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
## total post-warmup draws = 4000
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept -1.90 0.06 -2.02 -1.78 1.00 3573 2813
## Neurot 0.13 0.02 0.09 0.17 1.00 3605 2933
##
## Family Specific Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma 0.37 0.01 0.35 0.40 1.00 3940 2651
##
## Draws were sampled using sample(hmc). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).results_LS[4,6] <- extract_param(m_rsd_pa, "b_Neurot")4.2 Weighted RVI
Start$mean_NA <- NA
for (i in 1:nrow(Start)) {
if (Start$session[i] %in% Daily$session) {
if (
sum(!is.na(Daily$na[Daily$session == Start$session[i]])) >= 5
) {
Start$mean_NA[i] <- mean(Daily$na[Daily$session == Start$session[i]],
na.rm = T)
}
}
}
Start$mean_PA <- NA
for (i in 1:nrow(Start)) {
if (Start$session[i] %in% Daily$session) {
if (
sum(!is.na(Daily$pa[Daily$session == Start$session[i]])) >= 5
) {
Start$mean_PA[i] <- mean(Daily$pa[Daily$session == Start$session[i]],
na.rm = T)
}
}
}Start$weight_NA <- NA
for (i in 1:nrow(Start)) {
if (Start$session[i] %in% Daily$session) {
if (!is.na(Start$mean_NA[i])) {
# If mean is not missing
Start$weight_NA[i] <- maximumSD(Start$mean_NA[i], # mean
1, # min
5, # max
sum(!is.na(Daily$na[Daily$session == Start$session[i]]))
)
# W as reported in paper
Start$weight_NA[i] <- Start$weight_NA[i]^2
}
}
}
m_rsd_na_w <- brm(logrsd_n| weights(weight_NA) ~ Neurot, data= Start,
file = "na_logrsd_w")## Warning: Rows containing NAs were excluded from the model.print(m_rsd_na_w)## Family: gaussian
## Links: mu = identity; sigma = identity
## Formula: logrsd_n | weights(weight_NA) ~ Neurot
## Data: Start (Number of observations: 441)
## Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
## total post-warmup draws = 4000
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept -1.49 0.05 -1.58 -1.39 1.00 3023 2871
## Neurot 0.00 0.01 -0.02 0.03 1.00 3438 3188
##
## Family Specific Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma 0.38 0.01 0.37 0.40 1.00 3483 2857
##
## Draws were sampled using sample(hmc). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).results_LS[5,3] <- extract_param(m_rsd_na_w, "b_Neurot")
Start$weight_PA <- NA
for (i in 1:nrow(Start)) {
if (Start$session[i] %in% Daily$session) {
if (!is.na(Start$mean_PA[i])) {
Start$weight_PA[i] <- maximumSD(Start$mean_PA[i],
1,
5,
sum(!is.na(Daily$pa[Daily$session == Start$session[i]]))
)
Start$weight_PA[i] <- Start$weight_PA[i]^2
}
}
}
m_rsd_Pa_w <- brm(logrsd_p| weights(weight_PA) ~ Neurot, data= Start,
file = "pa_logrsd_w")## Warning: Rows containing NAs were excluded from the model.print(m_rsd_Pa_w)## Family: gaussian
## Links: mu = identity; sigma = identity
## Formula: logrsd_p | weights(weight_PA) ~ Neurot
## Data: Start (Number of observations: 437)
## Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
## total post-warmup draws = 4000
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept -1.93 0.03 -2.00 -1.87 1.00 4205 2698
## Neurot 0.14 0.01 0.12 0.16 1.00 4284 2896
##
## Family Specific Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma 0.36 0.01 0.34 0.37 1.00 3632 2833
##
## Draws were sampled using sample(hmc). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).results_LS[5,6] <- extract_param(m_rsd_Pa_w, "b_Neurot")5 SD
Start$sd_NA <- NA
for (i in 1:nrow(Start)) {
if (Start$session[i] %in% Daily$session) {
if (
sum(!is.na(Daily$na[Daily$session == Start$session[i]])) >= 5
) {
Start$sd_NA[i] <- sd(Daily$na[Daily$session == Start$session[i]],
na.rm = T)
}
}
}
Start$sd_PA <- NA
for (i in 1:nrow(Start)) {
if (Start$session[i] %in% Daily$session) {
if (
sum(!is.na(Daily$pa[Daily$session == Start$session[i]])) >= 5
) {
Start$sd_PA[i] <- sd(Daily$pa[Daily$session == Start$session[i]],
na.rm = T)
}
}
}
mean(Start$sd_NA, na.rm = T)## [1] 0.3728279mean(Start$sd_PA, na.rm= T)## [1] 0.437996Start$sd_PA[Start$sd_PA == 0] <- NA
Start$logsd_NA <- log(Start$sd_NA)
Start$logsd_PA <- log(Start$sd_PA)
m_sd_na <- brm(logsd_NA ~ Neurot, data= Start,
file = "na_logsd")## Warning: Rows containing NAs were excluded from the model.m_sd_na## Family: gaussian
## Links: mu = identity; sigma = identity
## Formula: logsd_NA ~ Neurot
## Data: Start (Number of observations: 441)
## Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
## total post-warmup draws = 4000
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept -1.76 0.07 -1.90 -1.62 1.00 3989 3128
## Neurot 0.23 0.02 0.18 0.28 1.00 3829 2911
##
## Family Specific Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma 0.43 0.01 0.40 0.46 1.00 3480 2888
##
## Draws were sampled using sample(hmc). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).results_LS[6,3] <- extract_param(m_sd_na, "b_Neurot")
m_sd_pa <- brm(logsd_PA ~ Neurot, data= Start,
file = "pa_logsd")## Warning: Rows containing NAs were excluded from the model.m_sd_pa## Family: gaussian
## Links: mu = identity; sigma = identity
## Formula: logsd_PA ~ Neurot
## Data: Start (Number of observations: 437)
## Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
## total post-warmup draws = 4000
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept -1.30 0.06 -1.42 -1.18 1.00 4939 3146
## Neurot 0.14 0.02 0.10 0.18 1.00 5066 3060
##
## Family Specific Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma 0.37 0.01 0.35 0.40 1.00 4081 2874
##
## Draws were sampled using sample(hmc). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).results_LS[6,6] <- extract_param(m_sd_pa, "b_Neurot")6 Incremental Validity of SD
na_noneurot <- brm(bf(na | cens(Acens) ~ (1|session),
sigma ~ (1|session)), data = daily,
iter = 7000, warmup = 2000,chains = 4,
control = list(adapt_delta = .99), init = 0.1,
file = "na_noneurot13")
print(na_noneurot)## Family: gaussian
## Links: mu = identity; sigma = log
## Formula: na | cens(Acens) ~ (1 | session)
## sigma ~ (1 | session)
## Data: daily (Number of observations: 6063)
## Draws: 4 chains, each with iter = 7000; warmup = 2000; thin = 1;
## total post-warmup draws = 20000
##
## Group-Level Effects:
## ~session (Number of levels: 444)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept) 0.56 0.02 0.52 0.60 1.00 1808 4000
## sd(sigma_Intercept) 0.35 0.02 0.32 0.39 1.00 6791 11798
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept 1.85 0.03 1.79 1.90 1.00 802 1756
## sigma_Intercept -0.93 0.02 -0.97 -0.89 1.00 5971 10402
##
## Draws were sampled using sample(hmc). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).rans <- coef(na_noneurot, summary = T)
rans_i <- as.data.frame(rans$session[,,"Intercept"]) %>% tibble::rownames_to_column("session")
rans_s <- as.data.frame(rans$session[,,"sigma_Intercept"]) %>% tibble::rownames_to_column("session")
nrow(rans_s)## [1] 444nrow(rans_i)## [1] 444nrow(Start)## [1] 2126dat <- merge(rans_s, rans_i, all = T, by= "session")
dat <- merge(dat, Start, all = T, by= "session")
names(dat)[2] <- "Est.SD"
names(dat)[6] <- "Est.M"
fit1 <- lm(Neurot ~ Est.SD + Est.M , data=dat)
summary(fit1)##
## Call:
## lm(formula = Neurot ~ Est.SD + Est.M, data = dat)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.93705 -0.48706 0.00357 0.48500 2.41832
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.41720 0.16705 8.484 3.33e-16 ***
## Est.SD 0.27820 0.10914 2.549 0.0111 *
## Est.M 0.93622 0.05915 15.828 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.665 on 441 degrees of freedom
## (1682 observations deleted due to missingness)
## Multiple R-squared: 0.3935, Adjusted R-squared: 0.3907
## F-statistic: 143 on 2 and 441 DF, p-value: < 2.2e-16fit1.2 <- lm(Neurot ~ Est.M , data=dat)
summary(fit1.2)##
## Call:
## lm(formula = Neurot ~ Est.M, data = dat)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.84242 -0.49728 -0.00481 0.49628 2.55625
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.09986 0.11207 9.814 <2e-16 ***
## Est.M 0.96751 0.05822 16.618 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.6691 on 442 degrees of freedom
## (1682 observations deleted due to missingness)
## Multiple R-squared: 0.3845, Adjusted R-squared: 0.3831
## F-statistic: 276.1 on 1 and 442 DF, p-value: < 2.2e-16aov <- anova(fit1.2, fit1)
aov## Analysis of Variance Table
##
## Model 1: Neurot ~ Est.M
## Model 2: Neurot ~ Est.SD + Est.M
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 442 197.87
## 2 441 195.00 1 2.8731 6.4977 0.01114 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1summary(fit1)$r.squared-summary(fit1.2)$r.squared## [1] 0.008936786results_SDin <- data.frame(matrix(nrow = 1, ncol = 9))
names(results_SDin) <- c("Dataset","b_SD","Err.SD","p(b_SD)","b_M","Err.M","p(b_M)","ΔR²", "p")
results_SDin$Dataset <- "13"
results_SDin$`ΔR²` <- summary(fit1)$r.squared-summary(fit1.2)$r.squared
results_SDin$`p` <- aov$`Pr(>F)`[2]
results_SDin$Err.SD <- summary(fit1)$coefficients[2,2]
results_SDin$b_SD <- fit1$coefficients[2]
results_SDin$`p(b_SD)` <- summary(fit1)$coefficients[2,4]
results_SDin$b_M <- fit1$coefficients[3]
results_SDin$`p(b_M)` <- summary(fit1)$coefficients[3,4]
results_SDin$Err.M <- summary(fit1)$coefficients[3,2]
library("writexl")
write_xlsx(results_SDin,"~/results_SD13.xlsx")