[R-meta] How to calculate predictions adjusted for the random effects of Experiment nested into Publication
Martineau, Roger
roger@m@rt|ne@u @end|ng |rom AGR@GC@CA
Sun Feb 20 21:44:15 CET 2022
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>
More information about the R-sig-meta-analysis
mailing list