[R-sig-ME] Measurement error for mixed models

[Viechtbauer, Wolfgang (SP)]

I am referring to the variance of the random effect specified via 'random = ~ 1 | id' (in rma(), this gets added by default). You are fitting the model:

y_i = beta0 + beta1 * x_i + u_i + e_i,

where u_i ~ N(0, sigma^2) (or tau^2 in the notation of rma()) and e_i ~ N(0, merror_i), where merror_i is the known measurement variance of the ith observation. If the estimate of sigma^2 is larger than 0, then the marginal variance of y_i is positive, even if merror_i = 0 for some observations. But if sigma^2 is estimated to be 0, then the marginal variance of y_i is 0 for observations where merror_i = 0 and then we are in trouble, because then we get division by zero when computing the MLEs of beta0 and beta1.

Note that the code below *does* work, but whether you can fit the model or not depends on whether the estimate of sigma^2 is larger than 0. If you run it multiple times, sometimes this will be the case and sometimes not.

If you add some small value, say epsilon, to merror_i, then you have merror_i = epsilon for those observations where merror_i = 0 to begin with and merror_i = merror_i + epsilon otherwise. That can be a problem, because the weights given to the observations (when computing the MLEs) can then be very different, especially if the estimate of sigma^2 = 0. In that case, the ratio of max((merror_i + epsilon) / epsilon) may be so extreme that rma() issues a warning (that is what the "Ratio of largest to smallest sampling variance extremely large" warning means). That should be warning, not an error (it used to be an error in earlier version of metafor, but should now only be warning -- so make sure you use the latest version of metafor).

Best,
Wolfgang

-----Original Message-----
From: Krzysztof Bartoszek [mailto:krzbar using protonmail.ch] 
Sent: Friday, 13 September, 2019 11:01
To: Viechtbauer, Wolfgang (SP)
Cc: r-sig-mixed-models using r-project.org
Subject: RE: [R-sig-ME] Measurement error for mixed models

PS
The below code does not work directly
df_Dummy<-as.data.frame(matrix(rnorm(20),ncol=2,nrow=10))
df_Dummy$id<-as.factor(1:nrow(df_Dummy))
df_Dummy$merror<-rexp(10)
colnames(df_Dummy)[1:2]<-c("x","y")
df_Dummy$merror[1:4]<-0
library(metafor)
rma(y ~ x, merror, data=df_Dummy, method="ML")

## produces
## Error in rma(y ~ x, merror, data = df_Dummy, method = "ML") :
##   Division by zero when computing the inverse variance weights.
## In addition: Warning message:
## In rma(y ~ x, merror, data = df_Dummy, method = "ML") :
##   There are outcomes with non-positive sampling variances.
## Taking as I need for nlme::lme() ang mgcv::gamm()

df_Dummy$merror<-df_Dummy$merror+0.0000000001
rma(y ~ x, merror, data=df_Dummy, method="ML")

## produces
## Error in rma(y ~ x, merror, data = df_Dummy, method = "ML") :
##  Ratio of largest to smallest sampling variance extremely large.
## Cannot obtain stable results.
## and

df_Dummy$merror[1:4]<-0;df_Dummy$merror<-df_Dummy$merror+0.000001
rma(y ~ x, merror, data=df_Dummy, method="ML")

## works, hopefully the bounding away from 0 variance is not too
## large to have an effect on estimation

Krzysztof

‐‐‐‐‐‐‐ Original Message ‐‐‐‐‐‐‐
On Friday, September 13, 2019 8:30 AM, Krzysztof Bartoszek <krzbar using protonmail.ch> wrote:

> Thanks for the code and suggestions. Trying to run the analysis in metafor was on the TODO list. Thank you for saving me the time on figuring out the syntax.
> If I may ask for clarification, what do you mean by "the variance component is larger than 0" ?
>
> Best wishes
> Krzysztof
>
> ‐‐‐‐‐‐‐ Original Message ‐‐‐‐‐‐‐
> On Thursday, September 12, 2019 9:04 PM, Viechtbauer, Wolfgang (SP) wolfgang.viechtbauer using maastrichtuniversity.nl wrote:
>
> > Thanks for the feedback.
> > For what it's worth, you can fit the model with those 0 variances with metafor:
> > df_Dummy<-as.data.frame(matrix(rnorm(20),ncol=2,nrow=10))
> > df_Dummy$id<-as.factor(1:nrow(df_Dummy))
> > df_Dummy$merror<-rexp(10)
> > colnames(df_Dummy)[1:2]<-c("x","y")
> > df_Dummy$merror[1:4]<-0
> > library(metafor)
> > rma(y ~ x, merror, data=df_Dummy, method="ML")
> > This will work as long as the variance component is larger than 0.
> > Best,
> > Wolfgang