# [R-meta] Nonlinear meta-regression with factorial moderator

Cesar Terrer Moreno cesar.terrer at me.com
Sun Apr 1 12:23:15 CEST 2018

```Hi Wolfgang,

I have followed your code and adapted for a factor with 5 levels (data attached)

df <- read.csv("~/OneDrive/OneDrive - Universitat Autònoma de Barcelona/FACEreview/data/AGB_effects.csv",na.strings=c("",NA))
df <- df[complete.cases(df\$id),]
df2 <- filter(df, Site.Name != "Nevada FACE",!is.na(CNr)) %>%
rename(es=yi, var=vi)
df2\$var <- ifelse(is.na(df2\$var),max(df2\$var,na.rm=T), df2\$var)

am.df <- filter(df2,Myc=="AM") %>% mutate(Biome = relevel(Biome, ref= "Grassland")) %>% droplevels()

# Biome has 5 levels. Then create 4 dummy variables (so level 1 is the reference level). Let's call these dummies B2, B3, B4, B5
library(dummies)
am.df <- dummy.data.frame("Biome",data=am.df, sep = "")

nlfun <- function(x, B2, B3, B4, B5, p1, p2, p3, p4, p5, p6)
(p1 + B2*p2 + B3*p3 + B4*p4 + B5*p5) * exp(-p6*x)

# optimization function
llfun <- function(par, yi, vi, x, B2, B3, B4, B5, random=TRUE) {
p1 <- par[1]
p2 <- par[2]
p3 <- par[3]
p4 <- par[4]
p5 <- par[5]
p6 <- par[6]
if (random) {
tau2 <- exp(par[7])
} else {
tau2 <- 0
}
mu <- nlfun(x, B2, B3, B4, B5, p1, p2, p3, p4, p5, p6)
-sum(dnorm(yi, mean=mu, sd=sqrt(vi + tau2), log=TRUE))
}

# optimize
res <- optim(par=c(8,0,0,0,0,0.4,log(.01)), llfun, yi=am.df\$es, vi=am.df\$var, x=am.df\$CNr,
B2=am.df\$BiomeCropland, B3=am.df\$BiomeShrubland,
B4=am.df\$BiomeTemperate_Forest, B5=am.df\$BiomeTropical_Forest,
hessian=TRUE)
# back-transform log(tau2) to tau2
res\$par[7] <- exp(res\$par[7])
tau2_ME <- res\$par[7]

# fit model with tau2=0
res0 <- optim(par=c(8,0,0,0,0,0.4,log(.01)), llfun, yi=am.df\$es, vi=am.df\$var, x=am.df\$CNr,
B2=am.df\$BiomeCropland, B3=am.df\$BiomeShrubland,
B4=am.df\$BiomeTemperate_Forest, B5=am.df\$BiomeTropical_Forest,
random=FALSE, hessian=TRUE)

# LRT of H0: tau2=0
x2 <- -2 * (res0\$value - res\$value)
x2 ### test statistic
pchisq(x2, df=1, lower.tail=FALSE) ### p-value

# compute standard errors
vb <- chol2inv(chol(res\$hessian))

However, in this last step, I get the following error:
Error in chol.default(res\$hessian) :
the leading minor of order 6 is not positive definite

Google says this may be due to overfitting. What can I do?
Thanks
César

> On 24 Mar 2018, at 16:39, Viechtbauer Wolfgang (SP) <wolfgang.viechtbauer at maastrichtuniversity.nl> wrote:
>
> Let's say F has three levels. Then create two dummy variables for levels 2 and 3 (so level 1 is the reference level). Let's call these dummies F2 and F2. Then I think this should be it (I haven't tested this):
>
> nlfun <- function(x, F2, F3, p1, p2, p3, p4)
>  (p1 + F2*p2 + F3*p3) * exp(-p4*x)
>
> # optimization function
> llfun <- function(par, yi, vi, x, F2, F3, random=TRUE) {
>  p1 <- par[1]
>  p2 <- par[2]
>  p3 <- par[3]
>  p4 <- par[4]
>  if (random) {
>    tau2 <- exp(par[5])
>  } else {
>    tau2 <- 0
>  }
>  mu <- nlfun(x, F2, F3, p1, p2, p3, p4)
>  -sum(dnorm(yi, mean=mu, sd=sqrt(vi + tau2), log=TRUE))
> }
>
> # optimize
> res <- optim(par=c(8,0,0,0.4,log(.01)), llfun, yi=am.df\$es, vi=am.df\$var, x=am.df\$CNr, F2=am.df\$F2, F3=am.df\$F3, hessian=TRUE)
> # back-transform log(tau2) to tau2
> res\$par[5] <- exp(res\$par[5])
> tau2_ME <- res\$par[5]
>
> # fit model with tau2=0
> res0 <- optim(par=c(8,0,0,0.4), llfun, yi=am.df\$es, vi=am.df\$var, x=am.df\$CNr, F2=am.df\$F2, F3=am.df\$F3, random=FALSE, hessian=TRUE)
>
> # LRT of H0: tau2=0
> x2 <- -2 * (res0\$value - res\$value)
> x2 ### test statistic
> pchisq(x2, df=1, lower.tail=FALSE) ### p-value
>
> # compute standard errors
> vb <- chol2inv(chol(res\$hessian))
> se <- sqrt(diag(vb))
>
> # ignore SE for tau2
> se[5] <- NA
>
> # estimates, standard errors, z-values, and p-values
> data.frame(estimate=round(res\$par,4), se=round(se,4), zval=round(res\$par/se,3), pval=round(2*pnorm(abs(res\$par/se), lower.tail=FALSE),5), row.names=c("p1","p2","p3","p4","tau2"))
>
> # Function to predict the nonlinear effect size (ES)
> ESpred <- function(x, F2, F3) {
>  p1 <- res\$par[1]
>  p2 <- res\$par[2]
>  p3 <- res\$par[3]
>  p4 <- res\$par[4]
>  (p1 + F2*p2 + F3*p3)*exp(-p4*x)
> }
>
> # Function to predict the nonlinear standard error (SE)
> SEpred <- function(x, F2, F3) {
>  p1 <- res\$par[1]
>  p2 <- res\$par[2]
>  p3 <- res\$par[3]
>  p4 <- res\$par[4]
>  g <- matrix(c((1 + F2*p2 + F3*p3) * exp(-p4*x), (p1 + F2 + F3*p3) * exp(-p4*x), (p1 + F2*p2 + F3) * exp(-p4*x), -(p1 + F2*p2 + F3*p3)*exp(-p4*x)*x), ncol=1, nrow=2)
>  vb <- vb[1:4,1:4]
>  c(sqrt(t(g) %*% vb %*% g))
> }
>
> Best,
> Wolfgang
>
> -----Original Message-----
> From: R-sig-meta-analysis [mailto:r-sig-meta-analysis-bounces at r-project.org] On Behalf Of Cesar Terrer Moreno
> Sent: Monday, 19 March, 2018 15:57
> To: r-sig-meta-analysis at r-project.org
> Subject: [R-meta] Nonlinear meta-regression with factorial moderator
>
> Dear all,
>
> I have an effect size that is driven by two moderators: a continuous “C” and a factorial “F”. Let’s say, in this case, that C is temperature, while F is ecosystem type (e.g. grassland, forest, cropland):
>
> rma(ES, VAR, data= dat, mods= ~ 1 + C + F)
>
> The relationship between ES and C is nonlinear, with high ES when C is low, but ES quickly approaching 0 when C becomes higher towards infinite, but never becoming negative.
>
> Previously, you guys helped me define a nonlinear meta-regression ES ~ C like: ES ~ p1 × exp(-p2*C), including a function to predict the nonlinear ES and SE.
>
> nlfun <- function(x, p1, p2)
>  p1 * exp(-p2*x)
>
> # optimization function
> llfun <- function(par, yi, vi, x, random=TRUE) {
>  p1 <- par[1]
>  p2 <- par[2]
>  if (random) {
>    tau2 <- exp(par[3])
>  } else {
>    tau2 <- 0
>  }
>  mu <- nlfun(x, p1, p2)
>  -sum(dnorm(yi, mean=mu, sd=sqrt(vi + tau2), log=TRUE))
> }
>
> # optimize
> res <- optim(par=c(8,0.4,log(.01)), llfun, yi=am.df\$es, vi=am.df\$var, x=am.df\$CNr, hessian=TRUE)
> # back-transform log(tau2) to tau2
> res\$par[3] <- exp(res\$par[3])
> tau2_ME <- res\$par[3]
>
> # fit model with tau2=0
> res0 <- optim(par=c(8,0.4), llfun, yi=am.df\$es, vi=am.df\$var, x=am.df\$CNr, random=FALSE, hessian=TRUE)
>
> # LRT of H0: tau2=0
> x2 <- -2 * (res0\$value - res\$value)
> x2 ### test statistic
> pchisq(x2, df=1, lower.tail=FALSE) ### p-value
>
> # compute standard errors
> vb <- chol2inv(chol(res\$hessian))
> se <- sqrt(diag(vb))
>
> # ignore SE for tau2
> se[3] <- NA
>
> # estimates, standard errors, z-values, and p-values
> data.frame(estimate=round(res\$par,4), se=round(se,4), zval=round(res\$par/se,3), pval=round(2*pnorm(abs(res\$par/se), lower.tail=FALSE),5), row.names=c("p1","p2","tau2"))
>
> # Function to predict the nonlinear effect size (ES)
> ESpred <- function(x) {
>  p1 <- res\$par[1]
>  p2 <- res\$par[2]
>  p1*exp(-p2*x)
> }
>
> # Function to predict the nonlinear standard error (SE)
> SEpred <- function(x) {
>  p1 <- res\$par[1]
>  p2 <- res\$par[2]
>  g <- matrix(c(exp(-p2*x), -p1*exp(-p2*x)*x), ncol=1, nrow=2)
>  vb <- vb[1:2,1:2]
>  c(sqrt(t(g) %*% vb %*% g))
> }
>
> I wouldn’t be able to figure out how to model the nonlinear relationship ES ~ C in my entire life, so I really appreciate the help.
>
> My question is: how can I incorporate a factorial moderator F in these functions to model ES ~ C + F.?
>