[R-meta] How to calculate predictions adjusted for the random effects of Experiment nested into Publication

Dear list members,

In the example dataset attached, the data have this hierarchical structure: 50 Treatments from 20 Experiments from 8 Publications. For each treatment, there is an observed (or reported) value and the corresponding predicted value by a model under evaluation. I can evaluate the model using the raw predictions and the predictions adjusted with the random effect of Experiment. However, I am not able to calculated the predictions adjusted with the random effect of Experiment nested into Publication.

I use the following steps:

?         The dataset is imported and renamed "d" and Residuals are calculated

#Calculate the residuals

d$Residuals <- d$Observed-d$Predicted

str(d)

'data.frame':      50 obs. of  8 variables:

 $ ID         : num  1 2 3 4 5 6 7 8 9 10 ...

 $ Publication: num  1 1 1 1 1 1 1 1 1 2 ...

 $ Experiment : num  1.1 1.1 1.1 1.2 1.2 1.2 1.3 1.3 1.3 2.1 ...

 $ Treatment  : num  1 2 3 4 5 6 7 8 9 10 ...

 $ SEM        : num  4 4 4 12 12 12 11 11 11 45 ...

 $ Observed   : num  141 65 97 178 210 185 147 141 174 217 ...

 $ Predicted  : num  136 71 78 158 200 255 141 112 162 259 ...

 $ Residuals  : num  5 -6 19 20 10 -70 6 29 12 -42 ...



?         Results from the RMSE function (see details of the function at the bottom):

#RMSE with raw predictions

RMSE(d$Observed, d$Predicted)

labels    output

1                  N  50.00000

2      Observed Mean 172.50000

3     Predicted Mean 186.24000

4               RMSE  42.14712

5       RMSE, % mean  24.43311 (% of observed mean)

6   Mean Bias, % MSE  10.62766

7  Slope Bias, % MSE   6.49443

8  Dispersion, % MSE  82.87792

9          Mean Bias -13.74000

10        Slope Bias  -0.09949

11       P-Mean Bias   0.01957

12      P-Slope Bias   0.05834

13               RSR   0.39921 (RMSE/sd Observed)

14               CCC   0.92198


?         Calculate predictions adjusted with the random effect of Experiment and weighting data with 1/SEM:

m1 <- rma.mv(Residuals  ~ 1, SEM, random = list(~1|Experiment, ~ 1|Treatment),

+              data=d, method = "REML", digits=5, sparse = TRUE)

#Source random effects, pull together with Experiments into a dataframe and name columns correctly

Study <- ranef.rma.mv(m1)$Experiment

Study <- subset (Study, select = c(intrcpt))

Experiment <- rownames(Study)

r <- cbind(Experiment, Study)

names(r) <- c("Experiment", "Study")

head(r)

Experiment  Study

1.1        1.1  16.01

1.2        1.2  -1.04

#Merge random effects back into data

d <- merge(d, r, by="Experiment")

#Calculate predictions with study adjustment

d$Predicted_Study <- d$Predicted+d$Study

str(d)

'data.frame':      50 obs. of  10 variables:

 $ Experiment     : num  1.1 1.1 1.1 1.2 1.2 1.2 1.3 1.3 1.3 2.1 ...

 $ ID             : num  1 2 3 4 5 6 7 8 9 10 ...

 $ Publication    : num  1 1 1 1 1 1 1 1 1 2 ...

 $ Treatment      : num  1 2 3 4 5 6 7 8 9 10 ...

 $ SEM            : num  4 4 4 12 12 12 11 11 11 10 ...

 $ Observed       : num  141 65 97 178 210 185 147 141 174 117 ...

 $ Predicted      : num  136 71 78 158 200 255 141 112 162 159 ...

 $ Residuals      : num  5 -6 19 20 10 -70 6 29 12 -42 ...

 $ Study          : num  16.01 16.01 16.01 -1.04 -1.04 ...

 $ Predicted_Study: num  152 87 94 157 199 ...

#RMSE with adjusted predictions

RMSE(d$Observed, d$Predicted_Study)

labels    output

1                  N  50.00000

2      Observed Mean 172.50000

3     Predicted Mean 184.66019

4               RMSE  20.45068

5       RMSE, % mean  11.85547

6   Mean Bias, % MSE  35.35618

7  Slope Bias, % MSE   0.48970

8  Dispersion, % MSE  64.15412

9          Mean Bias -12.16019

10        Slope Bias  -0.01367

11       P-Mean Bias   0.00010

12      P-Slope Bias   0.54783

13               RSR   0.19371

14               CCC   0.98101

The LRT test indicate that the hierarchical structure of data should be accounted for:

#  LRT test

Publication.Experiment.fit <- rma.mv(Residuals  ~ 1, SEM, random = list(~1|Publication, ~1|Experiment, ~ 1|Treatment),

+                                     data=d, method = "REML", digits=5, sparse = TRUE)

Publication.fit <- rma.mv(Residuals  ~ 1, SEM, random = list(~1|Publication, ~ 1|Treatment),

+                           data=d, method = "REML", digits=5, sparse = TRUE)

anova(Publication.Experiment.fit,Publication.fit)

df       AIC       BIC      AICc     logLik      LRT    pval         QE

Full     4 478.64551 486.21279 479.55460 -235.32275                  7451.72582

Reduced  3 502.23883 507.91429 502.77216 -248.11941 25.59332 <.00001 7451.72582

I would use the following model but I don't know how to proceed after that:

m2 <- rma.mv(Residuals  ~ 1, SEM, random = list(~1|Publication, ~1|Experiment, ~ 1|Treatment),

+              data=d, method = "REML", digits=5, sparse = TRUE)

#Source random effects, pull together with Experiments into a dataframe and name columns correctly

???

Thanks in advance,

Roger


#Specify RMSE function to be used ####
RMSE <- function(k,l) {
  if(is.list(l)==TRUE) {iter <- length(l)} else {iter=1}
  for(i in 1:iter) {
    if(is.list(k)==TRUE) { o <- k[[i]]} else {o <- k}
    if(is.list(l)==TRUE) { p <- l[[i]]} else {p <- l}
    d <- data.frame(o, p)
    d$res=d$o-d$p
    d <- subset(d, is.na(d$res)==FALSE)
    o <- d$o
    p <- d$p
    res <- d$res
    meano <- mean(o, na.rm=TRUE)
    meanp <- mean(p, na.rm=TRUE)
    PMeanBias <- t.test(res)$p.value
    PMeanBias <- ifelse(PMeanBias < 0.0001, 0.0001, PMeanBias)
    PSlope <- anova(lm(res~p))[1,5]
    PSlope <- ifelse(PSlope < 0.0001, 0.0001, PSlope)
    res2=res^2;
    rm=sqrt(mean(res2, na.rm=TRUE));
    uss=sum(res2, na.rm=TRUE);
    lo <- ifelse(is.na(o)==FALSE, 1, 0)
    n=sum(lo);
    meano=mean(o, na.rm=TRUE);
    mb=sum(res, na.rm=TRUE)/n;
    sse <- anova(lm(res~p))[2,2];
    msb <- mb^2;
    mspe <- rm^2;
    msre <- sse/n;
    msslope <- mspe-msre-msb;
    mean <- msb/mspe*100;
    slope <- msslope/mspe*100;
    residual <- msre/mspe*100;
    check <- mean+slope+residual
    rsr <- rm/sd(o, na.rm=TRUE)
    ccc <- epi.ccc(o,p)$rho.c[1]
    rmp = rm/meano*100
    mb <- mean(res, na.rm=TRUE)
    sb <- coef(lm(res~p))[2]
    output <- format(c(
      n, meano, meanp,
      rm, rmp,
      mean, slope, residual,
      mb, sb, PMeanBias, PSlope,
      rsr, ccc[,1]), digits=4, scientific=FALSE)
    labels <- c(
      "N", "Observed Mean", "Predicted Mean",
      "RMSE", "RMSE, % mean",
      "Mean Bias, % MSE", "Slope Bias, % MSE", "Dispersion, % MSE",
      "Mean Bias", "Slope Bias", "P-Mean Bias", "P-Slope Bias",
      "RSR", "CCC")
    if(i==1) {
      out <- data.frame(labels,output)
    }
    else {
      out <- cbind(out, output)
    }
  }
  return(out)
}



-------------- next part --------------
An HTML attachment was scrubbed...
URL: <https://stat.ethz.ch/pipermail/r-sig-meta-analysis/attachments/20220220/08c4a48c/attachment-0001.html>

-------------- next part --------------
A non-text attachment was scrubbed...
Name: Example.csv
Type: application/octet-stream
Size: 1193 bytes
Desc: Example.csv
URL: <https://stat.ethz.ch/pipermail/r-sig-meta-analysis/attachments/20220220/08c4a48c/attachment-0001.obj>