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

[Viechtbauer Wolfgang (SP)]

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.?

Many thanks in advance
César