[R-meta] testing for and visualizing correlation between dependent traits in bivariate meta-analysis

[Sigurd Einum]

Dear Wolfgang and James,
Thanks for your answers! I will start with the comments from James:
The data come from experiments that consists of a transfer of aquatic organisms from one salinity to a new one. Following the transfer, physiological traits (for these data this is NaK ATPase activity) are measured at subsequent time points. I then fit an asymptotic function to the the data for each experiment, and extract the estimates of two parameters (lambda which is the rate of change in the trait, and lFC.inf which is the asymptotic log-fold change in the trait) together with their variances and covariances. These then go into the meta-analyses. 

You are correct that there is skewness in the data, however the distribution of the BLUPs (which to my understanding is more crucial, https://stats.stackexchange.com/questions/451306/log-transforming-a-mean-what-to-do-with-sd) is more symmetric. Also, I am not sure how log-transforming the estimates would help, since I cannot obtain equivalent transformed measures of variances and covariances to use in the meta-analysis? Or is there a way to do this that I am not aware of?

And then, for Wolfgang’s comments:
Thank you for the pointer on the differences between the estimates of the correlation in the model and the correlation between the random effects. 

Yes, two of the species seem to have a very strong influence on the model. I don’t see any obvious issues with the data from these species. However, after you pointed this out, I went on to fit a model without a species random effect, and this one is actually much better in terms of AIC than one with a species random effect. So, given this apparent lack of a species effect my initial question about a correlation does not seem to make much sense anyways. 

mod1 <- rma.mv(yi = yi, V = V, mods = ~dep.trait,  data = new_data, random = list(~ dep.trait | species, ~1|id), struct ="DIAG",  method = "ML")
 
mod2 <- rma.mv(yi = yi, V = V, mods = ~dep.trait,  data = new_data, random = list(~dep.trait|id), struct ="DIAG",  method = "ML")

 AICc(mod1,mod2)
           df        AICc
mod1  5 11587.41340
mod2  4    11.92288

>-----Original Message-----
>From: R-sig-meta-analysis <r-sig-meta-analysis-bounces using r-project.org> On
>Behalf Of James Pustejovsky via R-sig-meta-analysis
>Sent: Thursday, October 31, 2024 5:46 PM
>To: R Special Interest Group for Meta-Analysis <r-sig-meta-analysis using r-
>project.org>
>Cc: James Pustejovsky <jepusto using gmail.com>
>Subject: Re: [R-meta] testing for and visualizing correlation between
>dependent traits in bivariate meta-analysis
>
>Hi Sigurd,
>To add to Wolfgang's note on outliers, could you share more about how these
>traits are measured and what the units of the effect sizes are? Just based on
>descriptive plots of the effect estimates, they appear to be strongly
>right-skewed:
>
>new_data_transf <-
>  new_data %>%
>  mutate(
>    vi = if_else(dep.trait == "lambda", v1i, v2i),
>    sei = sqrt(vi),
>    lnyi = log(yi),
>    sei_ln = sei / yi
>  )
>
>ggplot( new_data_transf, aes(yi, fill = dep.trait)) + geom_density(alpha =
>0.4) + facet_wrap(~ dep.trait)
>ggplot(new_data_transf, aes(yi, sei)) + geom_point() + facet_wrap(~
>dep.trait)
>
>Log transformation leads to a distribution much closer to symmetric:
>
>ggplot(new_data_transf, aes(yi, fill = dep.trait)) + geom_density(alpha =
>0.4) + facet_wrap(~ dep.trait) + scale_x_continuous(transform = "log")
>ggplot(new_data_transf, aes(lnyi, sei_ln)) + geom_point() + facet_wrap(~
>dep.trait) + scale_x_continuous()
>
>Best,
>James
>
>On Thu, Oct 31, 2024 at 7:46 AM Viechtbauer, Wolfgang (NP) via R-sig-meta-
>analysis <r-sig-meta-analysis using r-project.org> wrote:
>
>> Hi Sigurd,
>>
>> Thanks for the fully reproducible code/data.
>>
>> As for the difference between the parameter estimates from the model
>> and computations based on the random effects, you will find some
>> discussions of this in the archives. A quick search yielded:
>>
>> https://stat.ethz.ch/pipermail/r-sig-meta-analysis/2023-May/004627.htm
>> l
>>
>> but I vaguely remember other threads where this has come up (not in
>> connection with the correlation, but the variance).
>>
>> For you data, I examined the species-level standardized residuals:
>>
>> rstandard(mod1, cluster=new_data$species)$cluster
>>
>> and the Cook's distances:
>>
>> cds <- cooks.distance(mod1, cluster=new_data$species)
>> par(mar=c(4,12,2,2))
>> barplot(cds, horiz=TRUE, las=1, xlab="Cook's Distance")
>>
>> I would start by checking what is going with the data for these species.
>>
>> Best,
>> Wolfgang
>>
>> > -----Original Message-----
>> > From: R-sig-meta-analysis
>> > <r-sig-meta-analysis-bounces using r-project.org>
>> On Behalf
>> > Of Sigurd Einum via R-sig-meta-analysis
>> > Sent: Thursday, October 31, 2024 09:56
>> > To: R-sig-meta-analysis using r-project.org
>> > Cc: Sigurd Einum <sigurd.einum using ntnu.no>
>> > Subject: [R-meta] testing for and visualizing correlation between
>> dependent
>> > traits in bivariate meta-analysis
>> >
>> > I have a data set where two traits (lambda and lFC.inf) are measured
>> > simultaneously in different species, and the same species may be
>> observed in
>> > more than one experiment. I apply a bivariate model to these data,
>> > where
>> I want
>> > to check whether the two traits are correlated across species (e.g.
>> > does
>> species
>> > that have a large value for one trait have a small value for the other).
>> For
>> > this, I compare two models with struct ="UN" or struct = "DIAG",
>> > using
>> AICc.
>> > This gives strong support for the model with struct = "UN", with an
>> estimated
>> > negative correlation of -0.71.
>> >
>> > I want to visualize this pattern, and naively thought that I could
>> > just
>> obtain
>> > the species-level random effects for lambda and lFC.inf from this
>> > model (obtained using ranef) and make a plot between them. However,
>> > this plot
>> shows a
>> > very (unrealistic) tight relationship, and the negative correlation
>> between
>> > these is considerably stronger (-0.94) than what the model estimated.
>> > Interestingly, when I do the same for the model with struct =
>> > "DIAG", the correlation becomes more similar to the one estimated in
>> > the model
>> (-0.67).
>> >
>> > So my questions are (1) is my interpretation of the model comparison
>> correct?
>> > and (2) how do I best visualize the result in a figure?
>> >
>> > See below for script.
>> > Any input on this would be greatly appreciated.
>> >
>> > Best,
>> > Sigurd Einum
>>
>> _______________________________________________
>> R-sig-meta-analysis mailing list @ R-sig-meta-analysis using r-project.org
>> To manage your subscription to this mailing list, go to:
>> https://stat.ethz.ch/mailman/listinfo/r-sig-meta-analysis
>>
>
>	[[alternative HTML version deleted]]
>
>_______________________________________________
>R-sig-meta-analysis mailing list @ R-sig-meta-analysis using r-project.org To
>manage your subscription to this mailing list, go to:
>https://stat.ethz.ch/mailman/listinfo/r-sig-meta-analysis