Using posterior estimates for separate calculations in brms

Modeling
1

Dear friends - using your help I have fitted a double exponential decay curve to data from 24 rats belonging to three groups (control, diabetes, hyperglycemic) and the combined model below fits random effects to all four nonlinear parameters, and fixed effects to two of the parameters (for good reasons that are not relevant, I believe). The model is fine and repoduces all 24 decay curves over time (TT), and yields a perfect pp_check and correlation > 0.99 between fitted and measured concentrations.

[Pleprior2 <- c(prior(normal(300,300),nlpar=“A”,lb=10,ub=1000),
prior(normal(300,300),nlpar=“B”,lb=10,ub=1000),
prior(normal(5,5),nlpar=“g1”,lb=.0001),
prior(normal(5,5),nlpar=“g2”,lb=.0001)
)

system.time(Thin_raw_H3<- brm(bf(MMN~Aexp(-g1TT/100)+Bexp(-g2TT/100),
A + B ~ 1+(1|RAT), g1 + g2 ~ 1+TRT+(1|RAT),sigma ~ TRT ,nl=TRUE),data=dats,cores=4,prior = prior2,chains=1,control = list(adapt_delta = 0.95,max_treedepth = 20),iter = 2000))

Now the reason for it all is to find the glomerular filtration rate and this is had from
GFR <- V1*(A+B) / (A/G2 + B/G1) where Gs are 1e-2*gs and V1 is plasma volume known (approximately) for each rat.

Now, is it legitimate to take the four parameters from the fitted model directly as the posteriors with g1 and g2 different between TRT (groups)? Then multiplying with the mean of plasma volumes corresponding to the groups produces a sample of pseudoposterior estimates of GFR according to the groups. If so, how should I then compare say diabetics to controls? Using the perhaps contrived method sketched diabetics have higher GFR which I had also thought and the values are reasonable compared to the result obtained when fitting animail per animal - but this gives much worse overall fit.

Thank you so much and happy New Year

1 Like
2

If I understand you correctly, it is fine as long as you are using all the posterior draws rather than some summary of them. This GFR thing is basically a generated quantity that you are calculating afterward in R, so it has a posterior distribution for each rat. I am not sure why you propose multiplying by the mean of plasma volumes, when you can multiply by the measured plasma volume for each rat. In other words, you should have like 4000 draws from the posterior distribution of GFR for each rat, and you can do comparisons between treated rats and control rats.

1 Like
3

Thanks a lot!! My understanding was that if I had fixed effects for the four parameters, two of them the same throughout and two of them differing between the three groups I could use these group specific fixed effects as parameters in the formula for calculating GFR in R as you say, and ignore the random effects on all four parameters, Therefore I have only three groups of parameters, and my guess was that multiplying that by the mean V1 for each group should be equivalent to picking each individual rat with its random effect and individual V1. I notice that the random effects are required to get the perfect match for measured and fitted FITC concentration in each rat. To track GFR in each rat it may be more convenient to edit the stan code produced from brms than picking the random effects for each rat - also a reason why I had hoped the other way was legitimate :-)

4

If we are talking about GFR for a new rat, then you need to draw realizations of the rat-specific parameter(s) from its assumed distribution conditioning on posterior draws for the hyperparameters of that distribution (typically a Gaussian).

5

Thanks - but really I was interested in the 24 rats I knew rather than forming predictive values for new similar rats. If I make the 24 posterior GFR estimates as you suggest from the estimated fixed and random effects and then gather them according to the three treatment groups I can make plots of their distribution combined in each group - I assume it would be legitimate to subtract say control from diabetes and make 95% density estimates? Am I right in that?

Thanks

6

Yes

7

Thanks a lot - I found iterations for each of the four parameters for the 24 rats using fixed and random effects. It was rather clumsy and as suggested produced a different result than what I first came up with. Since there is no interaction between the random effects and the treatment effects I really had thought the random effects should cancel but no so. But still the fit for each rat is perfect, so the model is probably as it should be. Below is the clumsy code to extract the four parameters - thanks a lot

DDD <- as.data.frame(Thin_raw_H3)
getparams <- function(i,GRP) #i is rotnr 1:24
{A <- DDD$b_A_Intercept + DDD$sd_RAT__A_Intercept +
DDD[c(paste(“r_RAT__A.”,i,".Intercept.",sep=""))]
B <- DDD$b_B_Intercept + DDD$sd_RAT__B_Intercept +
DDD[c(paste(“r_RAT__B.”,i,".Intercept.",sep=""))]
DDD$Zero <- rep(0,length(DDD[,1]))
GLV <- “Zero” #default if GRP == Ctrl
if(GRP[i]==“Diab”) GLV <- “b_g1_TRTDIAB” else if (GRP[i]==“Hyper”) GLV <- “b_g1_TRTHYPER"
G1 <- DDD$b_g1_Intercept + DDD[c(GLV)] + DDD$sd_RAT__g1_Intercept +
DDD[c(paste(“r_RAT__g1.”,i,”.Intercept.",sep=""))]

GLV1 <- “Zero” # if GRP == Ctrl
if(GRP[i]==“Diab”) GLV1 <- “b_g2_TRTDIAB” else if (GRP[i]==“Hyper”) GLV1 <- “b_g2_TRTHYPER"
G2 <- DDD$b_g2_Intercept + DDD[c(GLV1)] + DDD$sd_RAT__g2_Intercept +
DDD[c(paste(“r_RAT__g2.”,i,”.Intercept.",sep=""))]
names(A) <- "A"
names(B)<- "B"
names(G1) <- "G1"
names(G2) <- "G2"
c(A,B,G1,G2)
}