tfjaeger
commented
5 years ago As I should have done straight away (apologies), I have put together a reproducible example below. The session info is still the same as above. For me, this code results in a fatal error (R crashes without any specific error message and R studio offers to restart R session) though the specific moment when it crashes seems to differ between runs.
library(tidyverse)
library(magrittr)
library(brms)
library(bridgesampling)
# I've just grabbed this data generator from an old project. It generates a function that can generate
# samples of repeated measures data (repeated subject and repeated items) for a by-2 between-subject
# and within-item design. The outcome is a binomially distributed response. The function is inelegant
# in a number of ways, but for the present purpose all it needs to do is to generate example data.
data.generator <- function(
estimates = cbind("Intercept" = 0, "Condition" = 0),
n.subj.perCondition = 1,
n.subj = n.subj.perCondition * 2, # Total number of subjects (for between subject design)
n.obs = 1, # Number of items per subject
subject.ranef.covar = matrix( # By default sets a simple by-subject covar with unit variance for intercept
c(1, rep(0, length(estimates)**2-1)),
nrow = length(estimates),
dimnames = list(colnames(estimates), colnames(estimates))
),
item.ranef.covar = matrix( # By default sets a by-item covar with unit variance for intercept and condition
c(1, 0, 0, 1),
nrow = length(estimates),
dimnames = list(colnames(estimates), colnames(estimates))
)
) {
require(mvtnorm)
require(magrittr)
require(dplyr)
n = n.obs * n.subj
d = expand.grid(Subject = 1:n.subj,
Item = 1:n.obs,
Intercept = 1)
d = d[order(d$Subject),]
d$Cond = factor(
sort(rep(
c("Control", "Manipulation"),
n / 2)),
levels = c("Manipulation", "Control"))
contrasts(d$Cond) = cbind("Condition" = c(1, 0))
d %<>%
mutate(
Condition = case_when(
Cond == "Manipulation" ~ 1,
T ~ 0
)
)
generate.data <- function() {
if (all(dim(item.ranef.covar) == 1)) {
item.adjustment <- data.frame(sapply(1:n.obs, function(x) { rnorm(n = 1, mean = 0, eval(item.ranef.covar**.5)) } ))
} else {
item.adjustment <- as.data.frame(t(sapply(1:n.obs, function(x) t(rmvnorm(n = 1, sigma = eval(item.ranef.covar))))))
}
names(item.adjustment) <- colnames(eval(item.ranef.covar))
item.adjustment$Item <- 1:n.obs
if (all(dim(subject.ranef.covar) == 1)) {
subject.adjustment <- data.frame(sapply(1:n.subj, function(x) { rnorm(n = 1, mean = 0, subject.ranef.covar**.5) } ))
} else{
subject.adjustment <- as.data.frame(t(sapply(1:n.subj, function(x) t(rmvnorm(n = 1, sigma = subject.ranef.covar)))))
}
names(subject.adjustment) <- colnames(subject.ranef.covar)
subject.adjustment$Subject <- 1:n.subj
d$y <- rowSums(d[,colnames(estimates)] * (
matrix(
rep(as.numeric(estimates), n),
byrow = T,
nrow = n
) +
item.adjustment[d$Item, colnames(estimates)] +
subject.adjustment[d$Subject, colnames(estimates)]
))
d$Item = factor(d$Item)
d$Subject = factor(d$Subject)
d$Outcome = sapply(d$y, FUN = function(x) { rbinom(1,1,plogis(x))})
return(d)
}
return(generate.data)
}
# This is the function that conducts the replication analysis that will be called in the loop
# below.
fit.replication.models <- function(d.sim1, d.sim2) {
require(brms)
require(bridgesampling)
niter = 10000 # number of iteractions for each chain of brm model
n.bs.rep = 1 # number of repetions for each bridge sampler
my.priors = c(
prior(normal(0, 10), class = b),
prior(cauchy(0, 2), class = sd),
prior(lkj(1), class = cor)
)
b.orig <- brm(Outcome ~ 1 + Cond +
(1 | Subject) + (1 + Cond | Item),
data = d.sim1, family = bernoulli, iter = niter, warmup = 1000, chains = 4, cores = 4,
prior = my.priors)
# Here I would usually obtain priors for the replication analysis from the posterior of the analysis of
# the original data. To keep this code simpler, I am just using the same priors (so we don't expect
# support for replication but since the point is just that this code crashes that's ok).
my.newpriors = my.priors
my.newpriors.woCondition = my.newpriors[-1,]
b.rep <- brm(Outcome ~ 1 + Cond + (1 + Cond | Item) + (1 | Subject),
data = d.sim2, family = bernoulli, iter = niter, warmup = 1000, chains = 4, cores = 4,
prior = my.priors, save_all_pars = T
)
b.rep.woCondition <- brm(Outcome ~ 1 + (1 + Cond | Item) + (1 | Subject),
data = d.sim2, family = bernoulli, iter = niter, warmup = 1000, chains = 4, cores = 4,
prior = my.newpriors.woCondition, save_all_pars = T
)
# Run bridge sampler
bs.rep = bridge_sampler(samples = b.rep, cores = 4, method = "warp3", maxiter = 1000, repetitions = n.bs.rep)
bs.rep.woCondition = bridge_sampler(samples = b.rep.woCondition, cores = 4, method = "warp3", maxiter = 1000, repetitions = n.bs.rep)
# Store some aspects of the summary of each model, as well as the BF based on the marginal likelihoods obtained
# by the bridge sampler.
s.orig = summary(b.orig)$fixed
s.rep = summary(b.rep)$fixed
simulation <- data.frame(
Int.1.estimate = s.orig[1,1],
Int.1.lower = s.orig[1,3],
Condition.1.estimate = s.orig[2,1],
Condition.1.lower = s.orig[2,3],
Int.2.estimate = s.rep[1,1],
Int.2.lower = s.rep[1,3],
Condition.2.estimate = s.rep[2,1],
Condition.2.lower = s.rep[2,3],
BF = bf(bs.rep, bs.rep.woCondition)$bf
)
return(simulation)
}SIMULATION
nsubj = 10
nitem = 20
# create a function that generates repeated measures data.
my.data = data.generator(
estimates = cbind("Intercept" = 1, "Condition" = 1),
n.subj.perCondition = nsubj,
n.obs = nitem
)
# run 10 iterations of the replication simulation by comparing a original and a replication
# experiment drawn from the same ground truth.
num.sims = 10
d.sim.orig_rep = plyr::rdply(.n = num.sims,
fit.replication.models(my.data(), my.data()),
.progress = "text")
singmann
I've been fitting models with brm() from brms that I then get the marginal likelihood for via bridge_sampler(). When I run a sequence of this (brm fit, bridge_sampler), R eventually crashes with a fatal error and needs to be restarted. This often happens on the 2nd or 3rd application of bridge_sampler().
Unfortunately, I cannot post an easily reproducible example, as I have not yet figured out what causes the problem, but I wanted to state it here in case others have encountered similar issues.
‘bridgesampling’ version 0.6-0
platform x86_64-apple-darwin15.6.0
arch x86_64
os darwin15.6.0
system x86_64, darwin15.6.0
status
major 3
minor 5.2
year 2018
month 12
day 20
svn rev 75870
language R
version.string R version 3.5.2 (2018-12-20) nickname Eggshell Igloo
run through RStudio Version 1.1.442 – © 2009-2018 RStudio, Inc