Modeling treatment effects with different variances for repeated measures within each treatment

Modeling
1

I have reason to believe that different treatments have influence on the variance of replicate measurements within a treatment. However I can’t figure out how my brm model formula should be in order to take this difference in variances into account.

To make things maybe a bit clearer I have simulated data which can be seen in below figure. This data is normally distributed with muA = 1, muB = 1.3 and muC = 0.8; the sd for all groups was 0.1. The sd for duplicate samples varies by group and were respectively sd_dupA = 0.05, sd_dupB = 0.01 and sd_dupC = 0.09.

synth_data
synth_data900×1050 34.5 KB

My data has three treatment groups (A,B and C), each treatment group has 6 samples and each sample is measured twice (Dupl). I have tried the following formula which do not seem to give me estimates on the different sd’s for the duplicates.

brmsformula(Value ~ Treatment + (1|Sample))
brmsformula(Value ~ Treatment + (1|gr(Sample, by = Treatment)))

Can someone please help me out on how to formulate a model in brms that does take into account the different variances for duplicates.

2

Have you tried bf(Value ~ Treatment + (1|Sample/Dupl))?

3

I don’t think this would work since this formula would make separate coefficients for each measurement i.e. Sample_Dupl.

4

Edit: Correction to previous post - you’re right, if you want within Sample variation (i.e. between duplicates) then you don’t want (1|Sample/Dupl), but I think you do want bf(Value ~ Treatment + (1|Sample), sigma ~ Treatment). You should be able to obtain different sigma for each level of Treatment, where sigma is the within Sample, i.e. between Dupl, variability.

Here is a reprex:

library(brms)
library(ggplot2)

n <- 3000
Treatment <- c(rep(rep(1:3), each=n/3))
Sample <- c(rep(rep(1:(n/2)), each=2))
Dupl <- c(rep(rep(1:2), times=n/2))
b <- c(0, -1, 1.5)
z <- rnorm(n/2, 0, 1)
mu <- b[Treatment] + z[Sample]
s_b <- c(0, -1, 1)
sigma <- exp(s_b[Treatment])
Value <- rnorm(n, mu, sigma)

d <- cbind.data.frame(Value,Treatment,Sample,Dupl)
d$Treatment <- factor(d$Treatment)
d$Sample <- factor(d$Sample)
d$Dupl <- factor(d$Dupl)
str(d)
d[1:10,]

ggplot(d, aes(x=Dupl, y=Value, group=Sample)) + geom_point() + geom_line() + facet_wrap(~Treatment, nrow=1) +  theme(legend.position="none")

m <- brm(bf(Value ~ Treatment + (1|Sample), sigma ~ Treatment), data=d, family=gaussian, cores=4)
m
1 Like
5

Hey thanks for your help, I think this is working! However I have a hard time interpreting the results. It seems that I have to exponentiate the sigma_TreatmentN parameters in order to get them on the right scale but I can’t get my head around of why that is. Could you clarify this?

Family: gaussian 
  Links: mu = identity; sigma = log 
Formula: Value ~ Treatment + (1 | Sample) 
         sigma ~ Treatment
   Data: d (Number of observations: 3000) 
  Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
         total post-warmup draws = 4000

Group-Level Effects: 
~Sample (Number of levels: 1500) 
              Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept)     0.99      0.03     0.94     1.05 1.01      927     1681

Population-Level Effects: 
                 Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept           -0.04      0.06    -0.15     0.07 1.00     2351     3139
sigma_Intercept      0.05      0.03    -0.01     0.11 1.00     2564     2718
Treatment2          -1.03      0.07    -1.17    -0.88 1.00      831     1810
Treatment3           1.62      0.11     1.40     1.84 1.00     4006     3139
sigma_Treatment2    -1.05      0.04    -1.13    -0.96 1.00     2356     3045
sigma_Treatment3     0.98      0.04     0.90     1.05 1.00     3380     3307

Draws were sampled using sampling(NUTS). 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).
6

Maybe I can jump in to address this quickly:
In order to keep \sigma positive all the computation is done on the unconstrained \log-scale. This means that also the coefficients will be on the \log-scale (\log(\sigma) = \sigma_{intercept}+\sigma_{Treatment_{k}}) . For your interpretation it means that you have to think about sigma_Intercept, sigma_Treatment2 and sigma_Treatment3 as multiplicative effects, because summing them up would be the equivalence of \sigma = e^{\sigma_{intercept}} \cdot e^{\sigma_{Treatment_{k}}}. Let’s say for Treatment 2 it would mean you have a \sigma of e^{0.05}\cdot e^{-1.05} = 0.368.

2 Likes
7

Thank you so much! this is exactly what I needed to know. Doing this transform recovers the sd’s used for simulation very nicely.