Additional Simulation with different Priors
1 Simulated 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. censored from 1-5
measure <- function(x) {
x[x < -2] <- -2
x[x > 2] <- 2
round(x,1) +3
}1.1 3 true models
True model 1: Associations with the mean (diary1) True model 2: Associations with SD (diary2) True 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: default priors vs. normally distributed priors
2.1 Simulation 1 (effect on mean)
prior1 <- c(prior("normal(0,1)", class = "b"),
prior("normal(0,1)", class = "b", dpar = "sigma"),
prior("normal(0,1)", class = "sd"))
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")
w1m3_prior <- brm(bf(Affect_m | cens(Acens) ~ neurot_m + (1|id),
sigma ~ neurot_m + (1|id)), data = diary1,
prior = prior1,init = 0.1,
iter = 6000, warmup = 2000,
file = "w1m3prior")
prior_summary(w1m3_prior)## prior class coef group resp dpar nlpar lb ub
## normal(0,1) b
## normal(0,1) b neurot_m
## normal(0,1) b sigma
## normal(0,1) b neurot_m sigma
## student_t(3, 1.6, 2.5) Intercept
## student_t(3, 0, 2.5) Intercept sigma
## normal(0,1) sd 0
## student_t(3, 0, 2.5) sd sigma 0
## normal(0,1) sd id 0
## normal(0,1) sd Intercept id 0
## student_t(3, 0, 2.5) sd id sigma 0
## student_t(3, 0, 2.5) sd Intercept id sigma 0
## source
## user
## (vectorized)
## user
## (vectorized)
## default
## default
## user
## default
## (vectorized)
## (vectorized)
## (vectorized)
## (vectorized)print(w1m3_prior)## 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 2900 5553
## sd(sigma_Intercept) 0.32 0.02 0.28 0.36 1.00 6411 9865
##
## 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 1866 3573
## sigma_Intercept -1.02 0.03 -1.07 -0.97 1.00 6301 9605
## neurot_m 0.50 0.02 0.46 0.55 1.00 2472 4665
## sigma_neurot_m 0.02 0.03 -0.03 0.08 1.00 7341 9851
##
## 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).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.37 1.00 2493 4999
## sd(sigma_Intercept) 0.32 0.02 0.28 0.36 1.00 5909 9422
##
## 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 1372 3511
## sigma_Intercept -1.02 0.03 -1.07 -0.97 1.00 5554 8118
## neurot_m 0.50 0.02 0.46 0.55 1.00 1790 3436
## sigma_neurot_m 0.02 0.03 -0.03 0.07 1.00 6097 8769
##
## 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)
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 3255 7040
## sd(sigma_Intercept) 0.31 0.02 0.28 0.35 1.00 9555 14943
##
## 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.00 1905 3821
## sigma_Intercept -1.03 0.02 -1.08 -0.98 1.00 6785 12692
## neurot_m 0.00 0.02 -0.04 0.05 1.00 2058 3863
## sigma_neurot_m 0.17 0.02 0.13 0.22 1.00 7577 13464
##
## 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)
w2m3_prior <- brm(bf(Affect_m | cens(Acens) ~ neurot_m + (1|id),
sigma ~ neurot_m + (1|id)), data = diary2,
prior = prior1,init = 0.1,
iter = 6000, warmup = 2000,
file = "w2m3prior")
prior_summary(w2m3_prior)## prior class coef group resp dpar nlpar lb ub
## normal(0,1) b
## normal(0,1) b neurot_m
## normal(0,1) b sigma
## normal(0,1) b neurot_m sigma
## student_t(3, 1.6, 2.5) Intercept
## student_t(3, 0, 2.5) Intercept sigma
## normal(0,1) sd 0
## student_t(3, 0, 2.5) sd sigma 0
## normal(0,1) sd id 0
## normal(0,1) sd Intercept id 0
## student_t(3, 0, 2.5) sd id sigma 0
## student_t(3, 0, 2.5) sd Intercept id sigma 0
## source
## user
## (vectorized)
## user
## (vectorized)
## default
## default
## user
## default
## (vectorized)
## (vectorized)
## (vectorized)
## (vectorized)print(w2m3_prior)## 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: 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.32 0.02 0.29 0.36 1.00 2400 4985
## sd(sigma_Intercept) 0.31 0.02 0.28 0.35 1.00 5823 9799
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept 1.61 0.02 1.56 1.65 1.00 1067 2120
## sigma_Intercept -1.03 0.02 -1.08 -0.98 1.00 4526 8431
## neurot_m 0.00 0.02 -0.04 0.05 1.00 1289 2677
## sigma_neurot_m 0.17 0.02 0.13 0.22 1.00 5052 8369
##
## 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).2.3 Simulation 3 (effects on both)
w3model_neuro3 <- brm(bf(Affect_m| cens(Acens) ~ neurot_m + (1|id),
sigma ~ neurot_m + (1|id)), data = diary3,
iter = 7000, warmup = 2000,chains = 4,
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: 4 chains, each with iter = 7000; warmup = 2000; thin = 1;
## total post-warmup draws = 20000
##
## 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.37 1.00 1794 4202
## sd(sigma_Intercept) 0.32 0.02 0.28 0.37 1.00 4401 9649
##
## 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.00 1039 2201
## sigma_Intercept -0.99 0.03 -1.04 -0.93 1.00 4180 7412
## neurot_m 0.52 0.03 0.47 0.57 1.00 1437 3186
## sigma_neurot_m 0.12 0.03 0.07 0.18 1.00 3704 8003
##
## 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).w3m3_prior <- brm(bf(Affect_m | cens(Acens) ~ neurot_m + (1|id),
sigma ~ neurot_m + (1|id)), data = diary3,
prior = prior1,init = 0.1,
iter = 6000, warmup = 2000,
file = "w3m3prior")
prior_summary(w3m3_prior)## prior class coef group resp dpar nlpar lb ub
## normal(0,1) b
## normal(0,1) b neurot_m
## normal(0,1) b sigma
## normal(0,1) b neurot_m sigma
## student_t(3, 1.6, 2.5) Intercept
## student_t(3, 0, 2.5) Intercept sigma
## normal(0,1) sd 0
## student_t(3, 0, 2.5) sd sigma 0
## normal(0,1) sd id 0
## normal(0,1) sd Intercept id 0
## student_t(3, 0, 2.5) sd id sigma 0
## student_t(3, 0, 2.5) sd Intercept id sigma 0
## source
## user
## (vectorized)
## user
## (vectorized)
## default
## default
## user
## default
## (vectorized)
## (vectorized)
## (vectorized)
## (vectorized)print(w3m3_prior)## 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: 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.37 1.00 3016 5683
## sd(sigma_Intercept) 0.32 0.02 0.28 0.37 1.00 5585 9713
##
## 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.00 1857 3621
## sigma_Intercept -0.98 0.03 -1.04 -0.93 1.00 6207 9273
## neurot_m 0.52 0.03 0.47 0.57 1.00 2657 5129
## sigma_neurot_m 0.12 0.03 0.07 0.18 1.00 6448 9891
##
## 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)
#qplot(diary3$neurot_m, diary3$Affect_m, alpha = I(0.1))