Main Simulation
1 Simulate Data
Create a sample with neuroticism scores (normal distributed)
set.seed(040519)
n <- 200
days_per_person <- 30
n_days <- n*days_per_person
people <- tibble(
id = 1:n,
neurot = rnorm(n),
latentaffect = rnorm(n),
latentaffectsd = rnorm(n),
)
sd(people$neurot)## [1] 1.020802This is how we simulate measurement. (measure: censored 1-5)
measure <- function(x) {
x[x < -2] <- -2
x[x > 2] <- 2
round(x,1) +3
}1.1 3 true models
model 1: association with the mean (diary1) model 2: association with variability (diary2) model 3: both(diary3)
diary1 <- people %>% full_join(
tibble(
id = rep(1:n, each = days_per_person),
), by = 'id') %>%
mutate(
Aff1 = -1.4 + 0.5* neurot + 0.3* latentaffect +
rnorm(n_days, mean =0 , sd = exp(-1.1 + 0 * neurot + 0.3* latentaffectsd))
)
qplot(diary1$neurot, diary1$Aff1, alpha = I(0.1)) + geom_hline(yintercept = -2, linetype = "dashed")## Warning: `qplot()` was deprecated in ggplot2 3.4.0.
qplot(diary1$Aff1)## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
qplot(measure(diary1$Aff1), binwidth = .1)
sd(diary1$Aff1[6001:9000])## [1] NAdiary1 %>% group_by(id, neurot) %>%
summarise(Aff1 = mean(Aff1)) %>% ungroup() %>% summarise(cor(Aff1, neurot))## `summarise()` has grouped output by 'id'. You can override using the `.groups`
## argument.## # A tibble: 1 × 1
## `cor(Aff1, neurot)`
## <dbl>
## 1 0.842sd(diary1$neurot)## [1] 1.018331sd(diary1$Aff1)## [1] 0.7062864hist(diary1$Aff1)
qplot(diary1$neurot, diary1$Aff1)
diary2 <- people %>% full_join(
tibble(
id = rep(1:n, each = days_per_person),
), by = 'id') %>%
mutate(
Aff2 = -1.4 + 0 * neurot + 0.3* latentaffect +
rnorm(n_days, mean = 0, sd = exp(-1.1 +0.15 * neurot + 0.3*latentaffectsd))
)
sd(diary2$Aff2)## [1] 0.5024995qplot(diary2$Aff2)## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
qplot(diary2$neurot, diary2$Aff2, alpha = I(0.1)) + geom_hline(yintercept = -2, linetype = "dashed")
qplot(measure(diary2$Aff2), binwidth = .1)
diary3 <- people %>% full_join(
tibble(
id = rep(1:n, each = days_per_person),
), by = 'id') %>%
mutate(
Aff3 = -1.4 + 0.5 * neurot + 0.3*latentaffect +
rnorm(n_days, mean = 0, sd = exp(-1.1 + 0.15 * neurot + 0.3* latentaffectsd))
)
sd(diary3$neurot)## [1] 1.018331sd(diary3$Aff3)## [1] 0.7144778qplot(diary3$neurot, diary3$Aff3, alpha = I(0.1)) + geom_hline(yintercept = -2, linetype = "dashed")
qplot(measure(diary3$Aff3), binwidth = .1)
sd(measure(diary3$Aff3), na.rm=T)## [1] 0.6146095qplot(diary3$Aff3, binwidth = .1)
Add measured Affect to all three Simulations
diary1 <- diary1 %>%
mutate(
Affect_m = measure(Aff1)
)
sd(diary1$Affect_m)## [1] 0.5993947round(cor(diary1 %>% select(Aff1, Affect_m)),2) ## Aff1 Affect_m
## Aff1 1.00 0.96
## Affect_m 0.96 1.00qplot(diary1$Affect_m, binwidth=.1)
diary2 <- diary2 %>%
mutate(
Affect_m = measure(Aff2 )
)
sd(diary2$Affect_m)## [1] 0.4509627round(cor(diary2 %>% select(Aff2, Affect_m)),2) ## Aff2 Affect_m
## Aff2 1.00 0.97
## Affect_m 0.97 1.00qplot(diary2$Affect_m, binwidth=.1)
diary3 <- diary3 %>%
mutate(
Affect_m = measure(Aff3)
)
sd(diary3$Affect_m)## [1] 0.6146095#round(cor(diary3 %>% select(Aff3, Affect_m)),2)
qplot(diary3$Affect_m, binwidth=.1)
diary1$Acens <- case_when(diary1$Affect_m == 1 ~ "left",
diary1$Affect_m == 5 ~ "right",
TRUE ~ "none")
table(diary1$Acens)##
## left none
## 1246 4754diary2$Acens <- case_when(diary2$Affect_m == 1 ~ "left",
diary2$Affect_m == 5 ~ "right",
TRUE ~ "none")
table(diary2$Acens)##
## left none
## 735 5265diary3$Acens <- case_when(diary3$Affect_m == 1 ~ "left",
diary3$Affect_m == 5 ~ "right",
TRUE ~ "none")
table(diary3$Acens)##
## left none right
## 1361 4638 1Add measured neuroticism to all three Simulations
measure_n <- function(x) {
# expects that x is N(0,1)
x <- x
round(x,1)
}
diary1 <- diary1 %>%
mutate(
neurot_m = measure_n(neurot)
)
sd(diary1$neurot_m)## [1] 1.018262#round(cor(diary1 %>% select(neurot, neurot_m)),2)
qplot(diary1$neurot_m, binwidth=.1)
diary2 <- diary2 %>%
mutate(
neurot_m = measure_n(neurot)
)
sd(diary2$neurot_m)## [1] 1.018262#round(cor(diary2 %>% select(neurot, neurot_m)),2)
#qplot(diary2$neurot_m, binwidth=.1)
diary3 <- diary3 %>%
mutate(
neurot_m = measure_n(neurot)
)
sd(diary3$neurot_m)## [1] 1.018262#round(cor(diary3 %>% select(neurot, neurot_m)),2)
#qplot(diary3$neurot_m, binwidth=.1)2 Estimated models
model 1: naiv: only associations with the mean, normal distribution assumption model 2: associations with the mean, censored model 3: association with mean and variability, censored
2.1 Simulation 1 (effect on mean) BCLSM
w1model_neuro3 <- brm(bf(Affect_m | cens(Acens) ~ neurot_m + (1|id),
sigma ~ neurot_m + (1|id)), data = diary1,
iter = 6000, warmup = 2000, init = 0.1,
file = "w1model_neuro3")
print(w1model_neuro3)## Family: gaussian
## Links: mu = identity; sigma = log
## Formula: Affect_m | cens(Acens) ~ neurot_m + (1 | id)
## sigma ~ neurot_m + (1 | id)
## Data: diary1 (Number of observations: 6000)
## Draws: 4 chains, each with iter = 6000; warmup = 2000; thin = 1;
## total post-warmup draws = 16000
##
## Group-Level Effects:
## ~id (Number of levels: 200)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept) 0.33 0.02 0.29 0.36 1.00 2884 5365
## sd(sigma_Intercept) 0.32 0.02 0.28 0.36 1.00 5762 8747
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept 1.60 0.02 1.55 1.65 1.00 1813 3506
## sigma_Intercept -1.02 0.03 -1.07 -0.97 1.00 6030 8922
## neurot_m 0.50 0.02 0.46 0.55 1.00 2264 4147
## sigma_neurot_m 0.02 0.03 -0.03 0.08 1.00 7085 9626
##
## 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(w1model_neuro3)
2.2 Simulation 2 (effect on SD) BCLSM
w2model_neuro3 <- brm(bf(Affect_m | cens(Acens) ~ neurot_m + (1|id),
sigma ~ neurot_m + (1|id)), data = diary2,
control = list(adapt_delta = .99),chains = 8,
iter = 6000, warmup = 2000, init = 0.1,
file = "w2model_neuro3")
print(w2model_neuro3)## Family: gaussian
## Links: mu = identity; sigma = log
## Formula: Affect_m | cens(Acens) ~ neurot_m + (1 | id)
## sigma ~ neurot_m + (1 | id)
## Data: diary2 (Number of observations: 6000)
## Draws: 8 chains, each with iter = 6000; warmup = 2000; thin = 1;
## total post-warmup draws = 32000
##
## Group-Level Effects:
## ~id (Number of levels: 200)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept) 0.32 0.02 0.29 0.36 1.00 3288 7366
## sd(sigma_Intercept) 0.31 0.02 0.28 0.35 1.00 8970 14951
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept 1.60 0.02 1.56 1.65 1.01 1860 3868
## sigma_Intercept -1.03 0.02 -1.08 -0.98 1.00 6886 13048
## neurot_m 0.00 0.02 -0.04 0.05 1.00 2114 4438
## sigma_neurot_m 0.17 0.02 0.13 0.22 1.00 7359 13761
##
## 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(w2model_neuro3)
2.3 Simulation 3 (effects on both) BCLSM
w3model_neuro3 <- brm(bf(Affect_m| cens(Acens) ~ neurot_m + (1|id),
sigma ~ neurot_m + (1|id)), data = diary3,
iter = 7000, warmup = 2000,chains = 8,
control = list(adapt_delta = .99), inits = 0.1 , #options = list(adapt_delta = 0.99)
file = "w3model_neuro3")## Warning: Argument 'inits' is deprecated. Please use argument 'init' instead.print(w3model_neuro3)## Family: gaussian
## Links: mu = identity; sigma = log
## Formula: Affect_m | cens(Acens) ~ neurot_m + (1 | id)
## sigma ~ neurot_m + (1 | id)
## Data: diary3 (Number of observations: 6000)
## Draws: 8 chains, each with iter = 7000; warmup = 2000; thin = 1;
## total post-warmup draws = 40000
##
## Group-Level Effects:
## ~id (Number of levels: 200)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept) 0.33 0.02 0.30 0.37 1.00 4424 9617
## sd(sigma_Intercept) 0.32 0.02 0.28 0.36 1.00 11205 20632
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept 1.58 0.02 1.53 1.63 1.01 2669 5050
## sigma_Intercept -0.98 0.03 -1.04 -0.93 1.00 9782 17033
## neurot_m 0.52 0.03 0.47 0.57 1.00 3544 7172
## sigma_neurot_m 0.12 0.03 0.07 0.18 1.00 10156 17937
##
## 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(w3model_neuro3)
extract_param2 <- function(model, parameter) {
ci <- posterior_summary(model, variable = parameter)
est <- sprintf("%.4f %.4f %.4f", ci[,"Estimate"], ci[,"Q2.5"], ci[,"Q97.5"])
est
}
results_sim <- data.frame(matrix(nrow = 8,
ncol = 9+1))
names(results_sim) <- c("model", "w1_b_neuro", "w1_b_neuro_sigma", "w1_sigma",
"w2_b_neuro", "w2_b_neuro_sigma", "w2_sigma",
"w3_b_neuro", "w3_b_neuro_sigma", "w3_sigma"
)
results_sim$model <- c("model1", "model2", "BCLSM",
"RVI", "RVI_weight", "SD", "SD*", "BLSM")results_sim[3, "w1_b_neuro"] <- extract_param2(w1model_neuro3, "b_neurot_m")
results_sim[3, "w1_b_neuro_sigma"] <- extract_param2(w1model_neuro3, "b_sigma_neurot_m")
results_sim[3, "w1_sigma"] <- extract_param2(w1model_neuro3, "b_sigma_Intercept")results_sim[3, "w2_b_neuro"] <- extract_param2(w2model_neuro3, "b_neurot_m")
results_sim[3, "w2_b_neuro_sigma"] <- extract_param2(w2model_neuro3, "b_sigma_neurot_m")
results_sim[3, "w2_sigma"] <- extract_param2(w2model_neuro3, "b_sigma_Intercept")results_sim[3, "w3_b_neuro"] <- extract_param2(w3model_neuro3, "b_neurot_m")
results_sim[3, "w3_b_neuro_sigma"] <- extract_param2(w3model_neuro3, "b_sigma_neurot_m")
results_sim[3, "w3_sigma"] <- extract_param2(w3model_neuro3, "b_sigma_Intercept")3 RVI (Relative-Variability-Index)
3.1 Unweighted RVI for all three Simulations
# Neurot Measure
people <- people %>%
mutate(
neurot = measure_n(neurot)
)
id <- unique(diary1$id)
id <- as.data.frame(id)
people$RSD_d1 <- NA
for (i in 1:nrow(id)) {
if (id$id[i] %in% diary1$id) {
people$RSD_d1[i] <- relativeSD(diary1$Affect_m[diary1$id == id$id[i]],
1, 5)
}
} ## Warning in checkOutput(M, MIN, MAX): NaN returned. Data has a mean equal the
## minimum
## Warning in checkOutput(M, MIN, MAX): NaN returned. Data has a mean equal the
## minimum
## Warning in checkOutput(M, MIN, MAX): NaN returned. Data has a mean equal the
## minimumpeople$logrsd_d1 <- log(people$RSD_d1)
m_rsd_d1 <- brm(logrsd_d1 ~ neurot, data= people,
file="simneg1_logrsd_uw")## Warning: Rows containing NAs were excluded from the model.print(m_rsd_d1, digits=4)## Family: gaussian
## Links: mu = identity; sigma = identity
## Formula: logrsd_d1 ~ neurot
## Data: people (Number of observations: 197)
## 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.3997 0.0257 -1.4522 -1.3487 1.0007 4199 2860
## neurot -0.1933 0.0265 -0.2460 -0.1432 1.0008 4313 2920
##
## Family Specific Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma 0.3599 0.0181 0.3265 0.3982 1.0014 3535 2735
##
## 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_sim[4,3] <- extract_param2(m_rsd_d1, "b_neurot")
people$RSD_d2 <- NA
for (i in 1:nrow(id)) {
if (id$id[i] %in% diary2$id) {
people$RSD_d2[i] <- relativeSD(diary2$Affect_m[diary2$id == id$id[i]],
1, 5)
}
}
people$logrsd_d2 <- log(people$RSD_d2)
m_rsd_d2 <- brm( logrsd_d2~ neurot, data= people,
file="simneg2_logrsd_uw")
m_rsd_d2## Family: gaussian
## Links: mu = identity; sigma = identity
## Formula: logrsd_d2 ~ neurot
## Data: people (Number of observations: 200)
## 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.47 0.02 -1.51 -1.43 1.00 3876 2886
## neurot 0.14 0.02 0.09 0.18 1.00 4275 2773
##
## Family Specific Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma 0.32 0.02 0.29 0.35 1.00 4300 3055
##
## 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_sim[4,6] <- extract_param2(m_rsd_d2, "b_neurot")
people$RSD_d3 <- NA
for (i in 1:nrow(id)) {
if (id$id[i] %in% diary3$id) {
people$RSD_d3[i] <- relativeSD(diary3$Affect_m[diary3$id == id$id[i]],
1, 5)
}
} ## Warning in checkOutput(M, MIN, MAX): NaN returned. Data has a mean equal the
## minimum
## Warning in checkOutput(M, MIN, MAX): NaN returned. Data has a mean equal the
## minimum
## Warning in checkOutput(M, MIN, MAX): NaN returned. Data has a mean equal the
## minimum
## Warning in checkOutput(M, MIN, MAX): NaN returned. Data has a mean equal the
## minimum
## Warning in checkOutput(M, MIN, MAX): NaN returned. Data has a mean equal the
## minimum
## Warning in checkOutput(M, MIN, MAX): NaN returned. Data has a mean equal the
## minimumpeople$logrsd_d3 <- log(people$RSD_d3)
m_rsd_d3 <- brm(logrsd_d3 ~ neurot, data= people,
file="simneg3_logrsd_uw")## Warning: Rows containing NAs were excluded from the model.m_rsd_d3## Family: gaussian
## Links: mu = identity; sigma = identity
## Formula: logrsd_d3 ~ neurot
## Data: people (Number of observations: 194)
## 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.34 0.03 -1.40 -1.29 1.00 3447 2901
## neurot -0.14 0.03 -0.19 -0.08 1.00 3869 2940
##
## Family Specific Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma 0.37 0.02 0.34 0.41 1.00 3705 2702
##
## 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_sim[4,9] <- extract_param2(m_rsd_d3, "b_neurot")3.2 weighted RVI for all three Simulations
people$mean_Aff_d1 <- NA
for (i in 1:nrow(id)) {
if (id$id[i] %in% diary1$id) {
people$mean_Aff_d1[i] <- mean(diary1$Affect_m[diary1$id == id$id[i]],
na.rm = T)
}
}
range(people$mean_Aff_d1)## [1] 1.00 3.37people$mean_Aff_d2 <- NA
for (i in 1:nrow(id)) {
if (id$id[i] %in% diary2$id) {
people$mean_Aff_d2[i] <- mean(diary2$Affect_m[diary2$id == id$id[i]],
na.rm = T)
}
}
range(people$mean_Aff_d2)## [1] 1.01 2.48people$mean_Aff_d3 <- NA
for (i in 1:nrow(id)) {
if (id$id[i] %in% diary3$id) {
people$mean_Aff_d3[i] <- mean(diary3$Affect_m[diary3$id == id$id[i]],
na.rm = T)
}
}
range(people$mean_Aff_d3)## [1] 1.000000 3.056667people$weight_d1 <- NA
for (i in 1:nrow(id)) {
if (id$id[i] %in% diary1$id) {
people$weight_d1[i] <- maximumSD(people$mean_Aff_d1[i],
1, # Minimum
5, # Maximum
sum(!is.na(diary1$Affect_m[diary1$id == id$id[i]])) # Anzahl Beobachtungen in var eingeflossen/30
)
# W as reported in paper
people$weight_d1[i] <- people$weight_d1[i]^2
}
}
people$weight_d2 <- NA
for (i in 1:nrow(id)) {
if (id$id[i] %in% diary2$id) {
people$weight_d2[i] <- maximumSD(people$mean_Aff_d2[i],
1, # Minimum
5, # Maximum
sum(!is.na(diary2$Affect_m[diary2$id == id$id[i]])) # Anzahl Beobachtungen in var eingeflossen/30
)
# W as reported in paper
people$weight_d2[i] <- people$weight_d2[i]^2
}
}
people$weight_d3 <- NA
for (i in 1:nrow(id)) {
if (id$id[i] %in% diary3$id) {
people$weight_d3[i] <- maximumSD(people$mean_Aff_d3[i],
1, # Minimum
5, # Maximum
sum(!is.na(diary3$Affect_m[diary3$id == id$id[i]])) # Anzahl Beobachtungen in var eingeflossen/30
)
# W as reported in paper
people$weight_d3[i] <- people$weight_d3[i]^2
}
}m_rsd_d1_w <- brm(logrsd_d1| weights(weight_d1) ~ neurot, data= people,
file="simneg1_logrsd_w")## Warning: Rows containing NAs were excluded from the model.m_rsd_d1_w## Family: gaussian
## Links: mu = identity; sigma = identity
## Formula: logrsd_d1 | weights(weight_d1) ~ neurot
## Data: people (Number of observations: 197)
## 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.50 0.02 -1.54 -1.47 1.00 2784 2715
## neurot -0.08 0.02 -0.11 -0.04 1.00 2743 3030
##
## Family Specific Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma 0.33 0.01 0.31 0.35 1.00 3201 2845
##
## 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_sim[5,3] <- extract_param2(m_rsd_d1_w, "b_neurot")
m_rsd_d2_w <- brm(logrsd_d2| weights(weight_d2) ~ neurot, data= people,
file="simneg2_logrsd_w")
m_rsd_d2_w## Family: gaussian
## Links: mu = identity; sigma = identity
## Formula: logrsd_d2 | weights(weight_d2) ~ neurot
## Data: people (Number of observations: 200)
## 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.51 0.02 -1.55 -1.48 1.00 3911 2870
## neurot 0.15 0.02 0.12 0.18 1.00 4397 3178
##
## Family Specific Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma 0.31 0.01 0.29 0.33 1.00 3696 3010
##
## 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_sim[5,6] <- extract_param2(m_rsd_d2_w, "b_neurot")
m_rsd_d3_w <- brm(logrsd_d3| weights(weight_d3) ~ neurot, data= people,
file="simneg3_logrsd_w")## Warning: Rows containing NAs were excluded from the model.m_rsd_d3_w## Family: gaussian
## Links: mu = identity; sigma = identity
## Formula: logrsd_d3 | weights(weight_d3) ~ neurot
## Data: people (Number of observations: 194)
## 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.47 0.02 -1.51 -1.43 1.00 2825 2597
## neurot 0.02 0.02 -0.02 0.06 1.00 2837 2917
##
## Family Specific Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma 0.33 0.01 0.31 0.36 1.00 2869 2747
##
## 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_sim[5,9] <- extract_param2(m_rsd_d3_w, "b_neurot")4 SD
people$sd_Aff_d1 <- NA
for (i in 1:nrow(id)) {
if (id$id[i] %in% diary1$id) {
people$sd_Aff_d1[i] <- sd(diary1$Affect_m[diary1$id == id$id[i]],
na.rm = T)
}
}
people$sd_Aff_d2 <- NA
for (i in 1:nrow(id)) {
if (id$id[i] %in% diary2$id) {
people$sd_Aff_d2[i] <- sd(diary2$Affect_m[diary2$id == id$id[i]],
na.rm = T)
}
}
people$sd_Aff_d3 <- NA
for (i in 1:nrow(id)) {
if (id$id[i] %in% diary3$id) {
people$sd_Aff_d3[i] <- sd(diary3$Affect_m[diary3$id == id$id[i]],
na.rm = T)
}
}
people$sd_Aff_d1[people$sd_Aff_d1 == 0] <- NA
people$sd_Aff_d2[people$sd_Aff_d2 == 0] <- NA
people$sd_Aff_d3[people$sd_Aff_d3 == 0] <- NA
people$logsd_d1 <- log(people$sd_Aff_d1)
people$logsd_d2 <- log(people$sd_Aff_d2)
people$logsd_d3 <- log(people$sd_Aff_d3)
mean(people$sd_Aff_d1) ## [1] NAmean(people$sd_Aff_d2) ## [1] 0.3308355mean(people$sd_Aff_d3, na.rm = T) ## [1] 0.32808324.1 Regression with SD
m_sd_d1 <- brm(logsd_d1 ~ neurot, data= people,
file="simneg1_logsd")## Warning: Rows containing NAs were excluded from the model.m_sd_d1## Family: gaussian
## Links: mu = identity; sigma = identity
## Formula: logsd_d1 ~ neurot
## Data: people (Number of observations: 197)
## 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.03 -1.36 -1.24 1.00 4214 3013
## neurot 0.30 0.03 0.24 0.36 1.00 3792 2656
##
## Family Specific Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma 0.43 0.02 0.39 0.48 1.00 3049 2714
##
## 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_sim[6,3] <- extract_param2(m_sd_d1, "b_neurot")
m_sd_d2 <- brm(logsd_d2 ~ neurot, data= people,
file="simneg2_logsd")
m_sd_d2## Family: gaussian
## Links: mu = identity; sigma = identity
## Formula: logsd_d2 ~ neurot
## Data: people (Number of observations: 200)
## 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.18 0.03 -1.23 -1.13 1.00 4083 3019
## neurot 0.15 0.03 0.10 0.20 1.00 4058 2893
##
## Family Specific Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma 0.38 0.02 0.35 0.42 1.00 3982 3305
##
## 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_sim[6,6] <- extract_param2(m_sd_d2, "b_neurot")
m_sd_d3 <- brm(logsd_d3 ~ neurot, data= people,
file="simneg3_logsd")## Warning: Rows containing NAs were excluded from the model.m_sd_d3## Family: gaussian
## Links: mu = identity; sigma = identity
## Formula: logsd_d3 ~ neurot
## Data: people (Number of observations: 194)
## 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.28 0.03 -1.35 -1.22 1.00 3784 2824
## neurot 0.42 0.04 0.35 0.49 1.00 3803 2748
##
## Family Specific Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma 0.46 0.02 0.42 0.51 1.00 3773 2958
##
## 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_sim[6,9] <- extract_param2(m_sd_d3, "b_neurot")4.2 Regression with SD + controlling for mean values of negative Emotion
m_sd_d1c <- brm(logsd_d1 ~ neurot + mean_Aff_d1, data= people,
file="simneg1_logsd_m")## Warning: Rows containing NAs were excluded from the model.m_sd_d1c## Family: gaussian
## Links: mu = identity; sigma = identity
## Formula: logsd_d1 ~ neurot + mean_Aff_d1
## Data: people (Number of observations: 197)
## 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.87 0.18 -2.21 -1.52 1.00 2088 2373
## neurot 0.16 0.05 0.06 0.27 1.00 2060 2164
## mean_Aff_d1 0.34 0.10 0.14 0.54 1.00 2041 2227
##
## Family Specific Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma 0.42 0.02 0.38 0.47 1.00 3367 2460
##
## 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_sim[7,3] <- extract_param2(m_sd_d1c, "b_neurot")
m_sd_d2c <- brm(logsd_d2 ~ neurot + mean_Aff_d2, data= people,
file="simneg2_logsd_m")
m_sd_d2c## Family: gaussian
## Links: mu = identity; sigma = identity
## Formula: logsd_d2 ~ neurot + mean_Aff_d2
## Data: people (Number of observations: 200)
## 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 -2.02 0.14 -2.30 -1.74 1.00 4007 2995
## neurot 0.14 0.02 0.09 0.19 1.00 4595 3065
## mean_Aff_d2 0.51 0.08 0.35 0.68 1.00 3959 2923
##
## Family Specific Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma 0.35 0.02 0.32 0.39 1.00 4141 3076
##
## 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_sim[7,6] <- extract_param2(m_sd_d2c, "b_neurot")
m_sd_d3c <- brm(logsd_d3 ~ neurot + mean_Aff_d3, data= people,
file="simneg3_logsd_m")## Warning: Rows containing NAs were excluded from the model.m_sd_d3c## Family: gaussian
## Links: mu = identity; sigma = identity
## Formula: logsd_d3 ~ neurot + mean_Aff_d3
## Data: people (Number of observations: 194)
## 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.79 0.20 -2.20 -1.41 1.00 1865 2495
## neurot 0.29 0.06 0.18 0.41 1.00 1808 2168
## mean_Aff_d3 0.30 0.12 0.07 0.54 1.00 1814 2336
##
## Family Specific Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma 0.46 0.02 0.41 0.51 1.00 3313 2620
##
## 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_sim[7,9] <- extract_param2(m_sd_d3c, "b_neurot")5 BLSM without Censoring
5.1 Simulation 1 BLSM
people_1 <- people[order(people$mean_Aff_d1),]
people_1e <- people_1[-c(1:10),]
high_1 <- diary1[(diary1$id %in% people_1e$id),]
w1model_neuro4b <- brm(bf(Affect_m ~ neurot_m + (1|id),
sigma ~ neurot_m + (1|id)), data = high_1,
control = list(adapt_delta = .99999),chains = 4,
iter = 7000, warmup = 2000, init = 0.1,
file = "w1model_neuro4b")
print(w1model_neuro4b)## Family: gaussian
## Links: mu = identity; sigma = log
## Formula: Affect_m ~ neurot_m + (1 | id)
## sigma ~ neurot_m + (1 | id)
## Data: high_1 (Number of observations: 5700)
## Draws: 4 chains, each with iter = 7000; warmup = 2000; thin = 1;
## total post-warmup draws = 20000
##
## Group-Level Effects:
## ~id (Number of levels: 190)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept) 0.29 0.02 0.26 0.32 1.00 2177 4284
## sd(sigma_Intercept) 0.34 0.02 0.30 0.39 1.00 5256 9738
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept 1.68 0.02 1.64 1.72 1.00 1217 2785
## sigma_Intercept -1.23 0.03 -1.28 -1.18 1.00 3368 6121
## neurot_m 0.41 0.02 0.37 0.46 1.00 1197 2703
## sigma_neurot_m 0.25 0.03 0.19 0.30 1.00 3652 5806
##
## 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(w1model_neuro4b)
5.2 Simulation 2 BLSM
people_2 <- people[order(people$mean_Aff_d2),]
people_2e <- people_2[-c(1:10),]
high_2 <- diary2[(diary2$id %in% people_2e$id),]
w2model_neuro4b <- brm(bf(Affect_m ~ neurot_m + (1|id),
sigma ~ neurot_m + (1|id)), data = high_2,
control = list(adapt_delta = .9999),chains = 4,
iter = 6000, warmup = 2000, init = 0.1,
file = "w2model_neuro4b")
print(w2model_neuro4b)## Family: gaussian
## Links: mu = identity; sigma = log
## Formula: Affect_m ~ neurot_m + (1 | id)
## sigma ~ neurot_m + (1 | id)
## Data: high_2 (Number of observations: 5700)
## Draws: 4 chains, each with iter = 6000; warmup = 2000; thin = 1;
## total post-warmup draws = 16000
##
## Group-Level Effects:
## ~id (Number of levels: 190)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept) 0.27 0.01 0.24 0.29 1.00 1850 3665
## sd(sigma_Intercept) 0.28 0.02 0.25 0.32 1.00 4586 8137
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept 1.66 0.02 1.63 1.70 1.00 1040 2170
## sigma_Intercept -1.12 0.02 -1.16 -1.08 1.00 3518 6428
## neurot_m 0.02 0.02 -0.02 0.06 1.00 1113 1677
## sigma_neurot_m 0.15 0.02 0.11 0.19 1.00 3655 6710
##
## 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).5.3 Simulation 3 BLSM
people_3 <- people[order(people$mean_Aff_d3),]
people_3e <- people_3[-c(1:10),]
high_3 <- diary3[(diary3$id %in% people_3e$id),]
w3model_neuro4b <- brm(bf(Affect_m ~ neurot_m + (1|id),
sigma ~ neurot_m + (1|id)), data = high_3,
iter = 7000, warmup = 2000,chains = 4,
control = list(adapt_delta = .9999), inits = 0.1 , #options = list(adapt_delta = 0.99)
file = "w3model_neuro4b")## Warning: Argument 'inits' is deprecated. Please use argument 'init' instead.print(w3model_neuro4b)## Family: gaussian
## Links: mu = identity; sigma = log
## Formula: Affect_m ~ neurot_m + (1 | id)
## sigma ~ neurot_m + (1 | id)
## Data: high_3 (Number of observations: 5700)
## Draws: 4 chains, each with iter = 7000; warmup = 2000; thin = 1;
## total post-warmup draws = 20000
##
## Group-Level Effects:
## ~id (Number of levels: 190)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept) 0.28 0.02 0.25 0.31 1.00 2544 4975
## sd(sigma_Intercept) 0.36 0.02 0.32 0.41 1.00 4779 9496
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept 1.67 0.02 1.63 1.71 1.00 1133 2531
## sigma_Intercept -1.23 0.03 -1.28 -1.17 1.00 3711 6880
## neurot_m 0.42 0.02 0.38 0.46 1.01 1252 2757
## sigma_neurot_m 0.37 0.03 0.31 0.43 1.00 3722 6849
##
## 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(w3model_neuro4b)
results_sim[8, "w1_b_neuro"] <- extract_param2(w1model_neuro4b, "b_neurot_m")
results_sim[8, "w1_b_neuro_sigma"] <- extract_param2(w1model_neuro4b, "b_sigma_neurot_m")
results_sim[8, "w1_sigma"] <- extract_param2(w1model_neuro4b, "b_sigma_Intercept")
results_sim[8, "w2_b_neuro"] <- extract_param2(w2model_neuro4b, "b_neurot_m")
results_sim[8, "w2_b_neuro_sigma"] <- extract_param2(w2model_neuro4b, "b_sigma_neurot_m")
results_sim[8, "w2_sigma"] <- extract_param2(w2model_neuro4b, "b_sigma_Intercept")
results_sim[8, "w3_b_neuro"] <- extract_param2(w3model_neuro4b, "b_neurot_m")
results_sim[8, "w3_b_neuro_sigma"] <- extract_param2(w3model_neuro4b, "b_sigma_neurot_m")
results_sim[8, "w3_sigma"] <- extract_param2(w3model_neuro4b, "b_sigma_Intercept")