Load data
n <- 1000
set.seed(1)
vcs <- tibble(f0 = rnorm(n))
vcs$sex_c <- if_else(rnorm(n) > 0, 1, 0)
vcs$contra <- factor(if_else(vcs$sex_c == 1, "ep", "no"))
vcs$who_c <- if_else(rnorm(n) > 0, 1, 0)
vcs$who <- factor(if_else(vcs$who_c == 1, "extra", "in"))
vcs$dataset <- factor(sample(1:40, size = n, replace = T))
datasets <- tibble(dataset = factor(1:40), intercept = rnorm(40))
vcs <- vcs %>% left_join(datasets)
## Joining, by = "dataset"
xtabs(~ who + contra, vcs)
## contra
## who ep no
## extra 235 262
## in 246 257
warmup <- 2000
iter <- warmup + 2000
chains <- 4
control <- list(adapt_delta = 0.99)
priors <- c(
prior(normal(0, 1), class = b)
)
vcs <- vcs %>% mutate(
f0 = 0.3 * rnorm(n()) + -0.3 * sex_c + rnorm(n()),
dominance =
intercept +
0.6 * sex_c +
-0.6 * who_c +
-0.4 * f0 +
-0.2 * if_else(sex_c == 1, 1, 0) * (f0-1)^2 +
-0.2 * if_else(who_c == 1, 1, 0) * (f0)^2 +
0.1 * if_else(who_c + sex_c == 2, 1, 0) * (f0-1)^2 +
0.95 * rnorm(n())
)
library(gamm4)
## Loading required package: Matrix
##
## Attaching package: 'Matrix'
## The following object is masked from 'package:tidyr':
##
## expand
## Loading required package: lme4
##
## Attaching package: 'lme4'
## The following object is masked from 'package:brms':
##
## ngrps
## Loading required package: mgcv
## Loading required package: nlme
##
## Attaching package: 'nlme'
## The following object is masked from 'package:lme4':
##
## lmList
## The following object is masked from 'package:dplyr':
##
## collapse
## This is mgcv 1.8-27. For overview type 'help("mgcv-package")'.
##
## Attaching package: 'mgcv'
## The following objects are masked from 'package:brms':
##
## s, t2
## This is gamm4 0.2-5
gam_simple <- gamm4(dominance ~ who + contra + s(f0, by = interaction(who, contra), bs = "cc"), random = ~(1 | dataset), data = vcs)
par(mfrow=c(2,2))
plot(gam_simple$gam)
gam_brms <- brm(dominance ~ who + contra + s(f0, by = interaction(who, contra), bs = "cc") + (1 | dataset), data = vcs,
iter = iter, warmup = warmup, cores = chains, chains = chains,
file = "test_repro4",
# prior = priors,
control = control)
smos <- marginal_smooths(gam_brms)
smos$`mu: s(f0,by=interaction(who,contra),bs="cc")` %>%
ggplot(aes(f0, estimate__, colour = who, fill = who, ymin = lower__, ymax = upper__)) +
geom_smooth(stat = 'identity') +
facet_grid(who ~ contra, scales = "free_y") +
scale_color_viridis_d(guide = FALSE) +
scale_fill_viridis_d(guide = FALSE)
ggplot(vcs %>% drop_na(who), aes(f0, dominance)) +
# geom_jitter(aes(colour = who), alpha = 0.3) +
geom_smooth(aes(colour = who, fill = who), method = "gam",
formula = y ~ s(x, bs = "cc")) +
scale_color_viridis_d("Contra", breaks = c(-1,1), labels = c("female", "male")) +
scale_fill_viridis_d("Contra", breaks = c(-1,1), labels = c("female", "male")) +
facet_grid( who ~ contra)
summary(gam_brms)
## Family: gaussian
## Links: mu = identity; sigma = identity
## Formula: dominance ~ who + contra + s(f0, by = interaction(who, contra), bs = "cc") + (1 | dataset)
## Data: vcs (Number of observations: 1000)
## Samples: 4 chains, each with iter = 4000; warmup = 2000; thin = 1;
## total post-warmup samples = 8000
##
## Smooth Terms:
## Estimate Est.Error l-95% CI
## sds(sf0interactionwhocontraextra.ep_1) 0.17 0.09 0.06
## sds(sf0interactionwhocontrain.ep_1) 0.08 0.06 0.00
## sds(sf0interactionwhocontraextra.no_1) 0.17 0.09 0.07
## sds(sf0interactionwhocontrain.no_1) 0.30 0.17 0.09
## u-95% CI Eff.Sample Rhat
## sds(sf0interactionwhocontraextra.ep_1) 0.41 2478 1.00
## sds(sf0interactionwhocontrain.ep_1) 0.24 2817 1.00
## sds(sf0interactionwhocontraextra.no_1) 0.39 2884 1.00
## sds(sf0interactionwhocontrain.no_1) 0.73 1687 1.00
##
## Group-Level Effects:
## ~dataset (Number of levels: 40)
## Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
## sd(Intercept) 0.99 0.12 0.78 1.27 1351 1.00
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
## Intercept -0.29 0.17 -0.61 0.04 918 1.00
## whoin 0.68 0.07 0.55 0.82 12592 1.00
## contrano -0.22 0.07 -0.35 -0.09 14336 1.00
##
## Family Specific Parameters:
## Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
## sigma 1.00 0.02 0.96 1.05 10786 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).