I am doing a meta-regression with metafor package in R. The mixed-effect model for proportion is used to assess the linearity between study performed year and medication prevalence. Here below is my script in R:

model_A <- rma.glmm(xi=A, ni=Sample, measure="PLO", mods=~year)

And results I got from R are:

Mixed-Effects Model (k = 32; tau^2 estimator: ML)

tau^2 (estimated amount of residual heterogeneity):     1.6349
tau (square root of estimated tau^2 value):             1.2786
I^2 (residual heterogeneity / unaccounted variability): 99.40%
H^2 (unaccounted variability / sampling variability):   168.00

Tests for Residual Heterogeneity: 
Wld(df = 30) = 2221.4535, p-val < .0001
LRT(df = 30) = 3187.7073, p-val < .0001

Test of Moderators (coefficient(s) 2): 
QM(df = 1) = 22.7322, p-val < .0001

Model Results:

          estimate        se     zval    pval      ci.lb      ci.ub
intrcpt  -554.8145  116.4605  -4.7640  <.0001  -783.0728  -326.5561  ***
year        0.2767    0.0580   4.7678  <.0001     0.1630     0.3905  ***

Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 

Followed by this model, I would also like to perform a scatterplot in R. So my script is:

wi <- 0.5/sqrt(dat$vi)
preds <- predict(model_A, transf = transf.ilogit, addx=TRUE)
plot(year, transf.ilogit(dat$yi), cex=wi)
lines(year, preds$pred)

The plot I got is: enter image description here

Apparently, it doesn't seem right!. So my questions are:

  1. Did I use the right model with rma.glmm?

  2. How could I weight individual study (cex=wi?)? How to calculate standard error for individual study?

  3. How could I fit a right estimated line in scatterplot?

Many thanks.


Followed by Wolfgang's suggestions, I managed to rescale the bubble and get predicted line fitted (the model remains the same):

enter image description here

Obviously, the line wasn't straight! Should I change model into polynomial regression? Or is that normal with this graph?

I tried polynomial model like:

model1<-rma.glmm(xi=A, ni=Sample, measure="PLO", mods=~year+I(year^2))

The error came with "Error in print(model1) : error in evaluating the argument 'x' in selecting a method for function 'print': Error: object 'model1' not found"

And I tried another model:

model2: model2<-rma.glmm(xi=A, ni=Sample, measure="PLO", mods=~year+year^2)

I got exactly the same result as original model, which has only the year as covariate fitted. I am not sure where the problem is....

Many thanks!



1 Answer 1

  1. Assuming you are trying to model the relationship between year and the log odds of the outcome of interest using a logistic mixed-effects model, then yes, you used the right model.

  2. You may want to rescale wi a bit. Such as:

    wi <- 0.5 + 3.0 * (wi - min(wi))/(max(wi) - min(wi))
  3. Something like this should do:

    years <- 1998:2014
    preds <- predict(model_A, transf = transf.ilogit, newmods = years)
    plot(year, transf.ilogit(dat$yi), cex=wi)
    lines(years, preds$pred)
    lines(years, preds$ci.lb, lty="dashed")
    lines(years, preds$ci.ub, lty="dashed")

See also here for a related example:


  • $\begingroup$ You are right that I am investigating the relationship between year and the log-odds of medication prevalence. But my another question is: does this model also fit to non-linear meta regression? And furthermore, as you can see my results, I have a huge negative intercept by that model. I am thinking to change year into (year-1995) as my covariate. I found the same P value and coefficient. But I couldn't fit into the plot anymore... $\endgroup$
    – Min
    Feb 8, 2016 at 17:18
  • $\begingroup$ Furthermore, . I did check that example before I asked. But, I don't understand why I should weight study by: wi<-1/sqrt(vi) size<-0.5+3*(wi-min(wi))/(max(wi)-min(wi)) Thanks for your help! $\endgroup$
    – Min
    Feb 8, 2016 at 17:18
  • $\begingroup$ 1) You can always do polynomial regression (i.e., adding year squared, cubed, and so on to the model) if you want to model a non-linear relationship. If you want a truly non-linear logistic mixed-effects models, you will have to look elsewhere. 2) The huge negative intercept is because it is an extrapolation to the year 0. Indeed, by rescaling the year variable, you can make the intercept more interpretable. Not sure what "I couldn't fit into the plot anymore" means. 3) The rescaling is just a suggestion for the plot -- so that the points are not quite so huge. $\endgroup$
    – Wolfgang
    Feb 8, 2016 at 18:49
  • $\begingroup$ Thanks Wolfgang! I rescaled the year variable and now it is with more interpretable intercept. I managed to get predicted line fitted into this scatter plot and found the line wasn't a straight line (model<-rma.glmm(xi=A, ni=sample, measure="PLO", mods=~year)). Therefore, I am thinking to fit in polynomial regression. But I got some problem with fitting model. The model I used: model1<-rma.glmm(xi=A, ni=sample, measure="PLO", mods=~year+I(year^2)) and I got an error saying "Model matrix not of full rank. Cannot fit model." $\endgroup$
    – Min
    Feb 9, 2016 at 11:45
  • 1
    $\begingroup$ The line isn't straight because transf.ilogit is a non-linear transformation. On the log odds scale, you are modeling a linear relationship. If you do want to fit a polynomial model, then mods = ~ year + I(year^2) should work unless this yields a design matrix that is nearly singular -- which appears to be the case here. Centering year (e.g., at its minimum or mean) should help. If not, use mods = ~ poly(year, degree=2). $\endgroup$
    – Wolfgang
    Feb 9, 2016 at 18:27

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.