Dear all,
I am performing a meta-analysis in brms. My data is correlated, so I’m using a variance-covariance matrix. However, when I use the variance-covariance matrix and the conditional_effects
function, the following error appears: Error in get_dpar(object, dpar = dp) : !is.null(x) is not TRUE.
However, when I fit the model and I don’t use the variance-covariance matrix, I can plot the figures using the conditional_effects
function.
Hi and welcome! Would it be possible for you to post the model code you are running and either a snippet of the data or simulated data. Also your version of R and brms would be helpful.
thanks!
Hi, @Ara_Winter thanks for your answer.
- R version 3.6.1
- brms Version: 2.11.1
The dataset:
dataset.csv (660 Bytes)
The coding:
# covariance matrix of outcomes
cor_matrix <- function(x, r, v = rep(1, length(x)), na.rm = FALSE) {
mat <- diag(v)
se <- sqrt(v)
se[is.na(se)] <- 0
if (length(x) > 1L) {
for (i in 2:nrow(mat)) {
for (j in 1:(i-1)) {
if (x[i] == x[j]) {
mat[i, j] <- mat[j, i] <- r * se[i] * se[j]
}
}
}
}
dimnames(mat) <- list(1:nrow(mat), 1:ncol(mat))
if (na.rm) {
keep <- !is.na(diag(mat))
mat <- mat[keep, keep]
}
mat
}
effect_coding <- function(x, levels) {
# convenience function for effect (i.e. sum) coding
x <- factor(x, levels = levels)
contrasts(x) <- contr.sum(length(levels))
colnames(contrasts(x)) <- levels[-length(levels)]
x
}
prior <-c(prior(normal(0, 100), class = Intercept),
prior(normal(0, 10), class = b),
prior(cauchy(0, 10), class = sd))
v.m1 <- cor_matrix(dat$study_ID, r = 0.9, v = dat$vi)
cor_obs <- cor_matrix(dat$study_ID, r = 0.9)
m1 <- brm(yi ~ 1 + log_time + level + exer +
(1|study_ID),
family = gaussian(),
prior = prior, autocor = cor_fixed(v.m1),
cores = parallel::detectCores(),
control = list(adapt_delta = 0.99), data = dat)
conditions1 <- data.frame(level = 1)
plot(conditional_effects(m1, effects = "log_time:exer",
conditions = conditions1))
1 Like
Thanks! I am hoping someone more proficient in brms will be able to answer. I’ll see if I can tag some people in.
1 Like
I know the dev @paul.buerkner posts here but I wasn’t sure if there is a list of brms users who can also help answer this question.
Sorry, I have overlooked this post. Does it still fail in the latest github version of brms?
Note that, in the github version, I recommend specifying the autocor term directly in formula
via
... ~ ... + fcor(v.m1)
instead of in autocor
. Please also note that the residual SD sigma will be estimated by default in this case. To turns it off, specify
formula = bf(
yi ~ 1 + log_time + level + exer + (1|study_ID) + fcor(v.m1),
sigma = 1
)
1 Like
Dear @paul.buerkner, thanks for your response. Yes, it still fails in the latest GitHub version of brms. I followed your suggestion, using the formula function, as below:
formula = bf(
yi ~ 1 + log_time + level + exer + (1|study_ID) + fcor(v.m1), data = dat,
sigma = 1
)
and the following error message appeared:
Error: Expecting a single value when fixing parameter 'data.pred_ID'.
Ok. Can you please provide a minimal reproducible example. I don’t see definition of dat
in the provided code.
With the following code, it works for me:
# covariance matrix of outcomes
cor_matrix <- function(x, r, v = rep(1, length(x)), na.rm = FALSE) {
mat <- diag(v)
se <- sqrt(v)
se[is.na(se)] <- 0
if (length(x) > 1L) {
for (i in 2:nrow(mat)) {
for (j in 1:(i-1)) {
if (x[i] == x[j]) {
mat[i, j] <- mat[j, i] <- r * se[i] * se[j]
}
}
}
}
dimnames(mat) <- list(1:nrow(mat), 1:ncol(mat))
if (na.rm) {
keep <- !is.na(diag(mat))
mat <- mat[keep, keep]
}
mat
}
effect_coding <- function(x, levels) {
# convenience function for effect (i.e. sum) coding
x <- factor(x, levels = levels)
contrasts(x) <- contr.sum(length(levels))
colnames(contrasts(x)) <- levels[-length(levels)]
x
}
dat <- data.frame(
yi = rnorm(20),
level = factor(rep(0:1, each = 10)),
log_time = rnorm(20),
exer = rnorm(20),
study_ID = 1:20,
vi = 0.25
)
prior <-c(prior(normal(0, 100), class = Intercept),
prior(normal(0, 10), class = b),
prior(cauchy(0, 10), class = sd))
v.m1 <- cor_matrix(dat$study_ID, r = 0.9, v = dat$vi)
cor_obs <- cor_matrix(dat$study_ID, r = 0.9)
# brms 2.11.6
m1 <- brm(bf(
yi ~ 1 + log_time + level + exer + (1|study_ID) + fcor(v.m1),
sigma = 1
),
family = gaussian(),
prior = prior, chains = 1,
data2 = list(v.m1 = v.m1),
control = list(adapt_delta = 0.99), data = dat)
conditions1 <- data.frame(level = 1)
plot(conditional_effects(m1, effects = "log_time:exer",
conditions = conditions1))
2 Likes
Now I used the model you provided and replaced it with my data. The model worked perfectly and the figures were generated by conditional_effects
.
Thank you very much @paul.buerkner for your help.
1 Like