Measurement error for mixed models

This error-in-variables approach is not available in lme. I do have an R-based implementation of this for models with random intercepts. You can find this implementation at:

https://shiny.airast.org/METRICS/

And a complete tutorial is under the Help tab.

-----Original Message-----
From: R-sig-mixed-models <r-sig-mixed-models-bounces at r-project.org> On Behalf Of Krzysztof Bartoszek via R-sig-mixed-models
Sent: Wednesday, September 11, 2019 10:02 AM
To: r-sig-mixed-models at r-project.org
Subject: [R-sig-ME] Measurement error for mixed models

Dear all,
As far as I managed to see the weights parameter in nlme::lme(), mgcv::gamm(), gamm4::gamm4(), can be used to pass some specific residual variance structure based on nlme's varFunc class. I was wondering if the following variance structure is possible to be obtained from the already implemented instances in varClasses, or I will need to code it myself.

I want the variance of the response for observation i to be of the form v_i^2 = s^2 + s_i^2, where s^2 is a common for all observations unknown variability and s_i^2 is known, individual specific measurement error variance (can be 0).

Thank you
Best wishes
Krzysztof Bartoszek

Sent with ProtonMail Secure Email.

_______________________________________________
R-sig-mixed-models at r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models

Hi Krzysztof,

This can be done with lme(). Let si2 denote the known variance of observation i. With varFixed(), you can specify that the error variances are known up to a proportionality constant, sigma^2. With lmeControl(sigma=1), you can fix that proportionality constant to 1, so the error variances are equal to s2i. Then add random effects for each row of the dataset. That will be your s^2. So, putting this all together:

dat$id <- 1:nrow(dat)
lme(yi ~ x1 + ..., random = ~ 1 | id, 
    weights = varFixed(~ si2), 
    control = lmeControl(sigma = 1), 
    data = dat)

Interestingly, this is identical to what is done in meta-analysis in the 'random-effects model', where we have known sampling variances for the outcomes and we add a random effect for each row (i.e., study) to the model. You might find this of interest:

http://www.metafor-project.org/doku.php/tips:rma_vs_lm_lme_lmer

But also read the note at the end -- there is (was?) some issue with lme() when fixing sigma to a constant. Not sure if this has been fixed in the meantime.

Best,
Wolfgang

-----Original Message-----
From: R-sig-mixed-models [mailto:r-sig-mixed-models-bounces at r-project.org] On Behalf Of Krzysztof Bartoszek via R-sig-mixed-models
Sent: Wednesday, 11 September, 2019 16:02
To: r-sig-mixed-models at r-project.org
Subject: [R-sig-ME] Measurement error for mixed models

Dear all,
As far as I managed to see the weights parameter in nlme::lme(), mgcv::gamm(), gamm4::gamm4(), can be used to pass some specific residual variance structure based on nlme's varFunc class. I was wondering if the following variance structure is possible to be obtained from the already implemented instances in varClasses, or I will need to code it myself.

I want the variance of the response for observation i to be of the form v_i^2 = s^2 + s_i^2, where s^2 is a common for all observations unknown variability and s_i^2 is known, individual specific measurement error variance (can be 0).

Thank you
Best wishes
Krzysztof Bartoszek

Hello!
Thank you Wolfgang, it was really helpful.
It seems that the issue with sigma fixed in nlme::lme() still persists. With nlme_3.1-141 I reran the code by Maciej Beresewicz from https://bugs.r-project.org/bugzilla/show_bug.cgi?id=16975 and obtain the same output of nlme::lme() as he did. Hence, I am using method="ML" as suggested on https://www.jepusto.com/bug-in-nlme-with-fixed-sigma/

However, interestingly one cannot have 0 measurement error for some observations.

library(nlme)
library(mgcv)
library(gamm4)
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")
nlme::lme(y~x,random=~1|id, weights=varFixed(~merror), control=lmeControl(sigma=1), data=df_Dummy, method="ML") ## works fine
mgcv::gamm(y~x,random= list(id=~1), weights=varFixed(~merror), control=lmeControl(sigma=1), data=df_Dummy, family=gaussian, method="ML") ## works fine

df_Dummy$merror[1:4]<-0
nlme::lme(y~x,random=~1|id, weights=varFixed(~merror), control=lmeControl(sigma=1), data=df_Dummy, method="ML")
## produces:
## Error in MEestimate(lmeSt, grps) :
##   NA/NaN/Inf in foreign function call (arg 1)
## In addition: Warning messages:
## 1: In nlminb(c(oldPars), function(lmePars) -logLik(lmeSt, lmePars),  :
##   NA/NaN function evaluation
## 2: In nlminb(c(oldPars), function(lmePars) -logLik(lmeSt, lmePars),  :
##   NA/NaN function evaluation
## 3: In nlminb(c(oldPars), function(lmePars) -logLik(lmeSt, lmePars),  :
##   NA/NaN function evaluation
## 4: In nlminb(c(oldPars), function(lmePars) -logLik(lmeSt, lmePars),  :
##   NA/NaN function evaluation

mgcv::gamm(y~x,random= list(id=~1), weights=varFixed(~merror), control=lmeControl(sigma=1), data=df_Dummy, family=gaussian, method="ML")
## produces:
## Error in MEestimate(lmeSt, grps) :
##  NA/NaN/Inf in foreign function call (arg 1)
## In addition: Warning messages:
## 1: In nlminb(c(oldPars), function(lmePars) -logLik(lmeSt, lmePars),  :
##  NA/NaN function evaluation
## 2: In nlminb(c(oldPars), function(lmePars) -logLik(lmeSt, lmePars),  :
##  NA/NaN function evaluation
## 3: In nlminb(c(oldPars), function(lmePars) -logLik(lmeSt, lmePars),  :
##  NA/NaN function evaluation
## 4: In nlminb(c(oldPars), function(lmePars) -logLik(lmeSt, lmePars),  :
##  NA/NaN function evaluation

df_Dummy$merror<-df_Dummy$merror+0.000000001
nlme::lme(y~x,random=~1|id, weights=varFixed(~merror), control=lmeControl(sigma=1), data=df_Dummy, method="ML") ## works fine again
mgcv::gamm(y~x,random= list(id=~1), weights=varFixed(~merror), control=lmeControl(sigma=1), data=df_Dummy, family=gaussian, method="ML") ## works fine again

Also, unfortunately the same mechanism cannot be used in gamm4::gamm4(), if I am not missing something.

gamm4::gamm4(y~x,random= ~(1|id), weights=varFixed(~merror), control=lmeControl(sigma=1), data=df_Dummy, family=gaussian,methof="ML")
## produces:
## Error in model.frame.default(formula = y ~ x, data = df_Dummy, weights = varFixed(~merror),  :
##  variable lengths differ (found for '(weights)')

Best wishes
Krzysztof

??????? Original Message ???????

> Hi Krzysztof,
>
> This can be done with lme(). Let si2 denote the known variance of observation i. With varFixed(), you can specify that the error variances are known up to a proportionality constant, sigma^2. With lmeControl(sigma=1), you can fix that proportionality constant to 1, so the error variances are equal to s2i. Then add random effects for each row of the dataset. That will be your s^2. So, putting this all together:
>
> dat$id <- 1:nrow(dat)
> lme(yi ~ x1 + ..., random = ~ 1 | id,
> weights = varFixed(~ si2),
> control = lmeControl(sigma = 1),
> data = dat)
>
> Interestingly, this is identical to what is done in meta-analysis in the 'random-effects model', where we have known sampling variances for the outcomes and we add a random effect for each row (i.e., study) to the model. You might find this of interest:
>
> http://www.metafor-project.org/doku.php/tips:rma_vs_lm_lme_lmer
>
> But also read the note at the end -- there is (was?) some issue with lme() when fixing sigma to a constant. Not sure if this has been fixed in the meantime.
>
> Best,
> Wolfgang
>
> -----Original Message-----
> From: R-sig-mixed-models [mailto:r-sig-mixed-models-bounces at r-project.org] On Behalf Of Krzysztof Bartoszek via R-sig-mixed-models
> Sent: Wednesday, 11 September, 2019 16:02
> To: r-sig-mixed-models at r-project.org
> Subject: [R-sig-ME] Measurement error for mixed models
>
> Dear all,
> As far as I managed to see the weights parameter in nlme::lme(), mgcv::gamm(), gamm4::gamm4(), can be used to pass some specific residual variance structure based on nlme's varFunc class. I was wondering if the following variance structure is possible to be obtained from the already implemented instances in varClasses, or I will need to code it myself.
>
> I want the variance of the response for observation i to be of the form v_i^2 = s^2 + s_i^2, where s^2 is a common for all observations unknown variability and s_i^2 is known, individual specific measurement error variance (can be 0).
>
> Thank you
> Best wishes
> Krzysztof Bartoszek

</wolfgang.viechtbauer at maastrichtuniversity.nl>

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

-----Original Message-----
From: Krzysztof Bartoszek [mailto:krzbar at protonmail.ch] 
Sent: Thursday, 12 September, 2019 22:46
To: Viechtbauer, Wolfgang (SP)
Cc: r-sig-mixed-models at r-project.org
Subject: RE: [R-sig-ME] Measurement error for mixed models

Hello!
Thank you Wolfgang, it was really helpful.
It seems that the issue with sigma fixed in nlme::lme() still persists. With nlme_3.1-141 I reran the code by Maciej Beresewicz from https://bugs.r-project.org/bugzilla/show_bug.cgi?id=16975 and obtain the same output of nlme::lme() as he did. Hence, I am using method="ML" as suggested on https://www.jepusto.com/bug-in-nlme-with-fixed-sigma/

However, interestingly one cannot have 0 measurement error for some observations.

library(nlme)
library(mgcv)
library(gamm4)
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")
nlme::lme(y~x,random=~1|id, weights=varFixed(~merror), control=lmeControl(sigma=1), data=df_Dummy, method="ML") ## works fine
mgcv::gamm(y~x,random= list(id=~1), weights=varFixed(~merror), control=lmeControl(sigma=1), data=df_Dummy, family=gaussian, method="ML") ## works fine

df_Dummy$merror[1:4]<-0
nlme::lme(y~x,random=~1|id, weights=varFixed(~merror), control=lmeControl(sigma=1), data=df_Dummy, method="ML")
## produces:
## Error in MEestimate(lmeSt, grps) :
##   NA/NaN/Inf in foreign function call (arg 1)
## In addition: Warning messages:
## 1: In nlminb(c(oldPars), function(lmePars) -logLik(lmeSt, lmePars),  :
##   NA/NaN function evaluation
## 2: In nlminb(c(oldPars), function(lmePars) -logLik(lmeSt, lmePars),  :
##   NA/NaN function evaluation
## 3: In nlminb(c(oldPars), function(lmePars) -logLik(lmeSt, lmePars),  :
##   NA/NaN function evaluation
## 4: In nlminb(c(oldPars), function(lmePars) -logLik(lmeSt, lmePars),  :
##   NA/NaN function evaluation

mgcv::gamm(y~x,random= list(id=~1), weights=varFixed(~merror), control=lmeControl(sigma=1), data=df_Dummy, family=gaussian, method="ML")
## produces:
## Error in MEestimate(lmeSt, grps) :
##  NA/NaN/Inf in foreign function call (arg 1)
## In addition: Warning messages:
## 1: In nlminb(c(oldPars), function(lmePars) -logLik(lmeSt, lmePars),  :
##  NA/NaN function evaluation
## 2: In nlminb(c(oldPars), function(lmePars) -logLik(lmeSt, lmePars),  :
##  NA/NaN function evaluation
## 3: In nlminb(c(oldPars), function(lmePars) -logLik(lmeSt, lmePars),  :
##  NA/NaN function evaluation
## 4: In nlminb(c(oldPars), function(lmePars) -logLik(lmeSt, lmePars),  :
##  NA/NaN function evaluation

df_Dummy$merror<-df_Dummy$merror+0.000000001
nlme::lme(y~x,random=~1|id, weights=varFixed(~merror), control=lmeControl(sigma=1), data=df_Dummy, method="ML") ## works fine again
mgcv::gamm(y~x,random= list(id=~1), weights=varFixed(~merror), control=lmeControl(sigma=1), data=df_Dummy, family=gaussian, method="ML") ## works fine again

Also, unfortunately the same mechanism cannot be used in gamm4::gamm4(), if I am not missing something.

gamm4::gamm4(y~x,random= ~(1|id), weights=varFixed(~merror), control=lmeControl(sigma=1), data=df_Dummy, family=gaussian,methof="ML")
## produces:
## Error in model.frame.default(formula = y ~ x, data = df_Dummy, weights = varFixed(~merror),  :
##  variable lengths differ (found for '(weights)')

Best wishes
Krzysztof

??????? Original Message ???????

> Hi Krzysztof,
>
> This can be done with lme(). Let si2 denote the known variance of observation i. With varFixed(), you can specify that the error variances are known up to a proportionality constant, sigma^2. With lmeControl(sigma=1), you can fix that proportionality constant to 1, so the error variances are equal to s2i. Then add random effects for each row of the dataset. That will be your s^2. So, putting this all together:
>
> dat$id <- 1:nrow(dat)
> lme(yi ~ x1 + ..., random = ~ 1 | id,
> weights = varFixed(~ si2),
> control = lmeControl(sigma = 1),
> data = dat)
>
> Interestingly, this is identical to what is done in meta-analysis in the 'random-effects model', where we have known sampling variances for the outcomes and we add a random effect for each row (i.e., study) to the model. You might find this of interest:
>
> http://www.metafor-project.org/doku.php/tips:rma_vs_lm_lme_lmer
>
> But also read the note at the end -- there is (was?) some issue with lme() when fixing sigma to a constant. Not sure if this has been fixed in the meantime.
>
> Best,
> Wolfgang
>
> -----Original Message-----
> From: R-sig-mixed-models [mailto:r-sig-mixed-models-bounces at r-project.org] On Behalf Of Krzysztof Bartoszek via R-sig-mixed-models
> Sent: Wednesday, 11 September, 2019 16:02
> To: r-sig-mixed-models at r-project.org
> Subject: [R-sig-ME] Measurement error for mixed models
>
> Dear all,
> As far as I managed to see the weights parameter in nlme::lme(), mgcv::gamm(), gamm4::gamm4(), can be used to pass some specific residual variance structure based on nlme's varFunc class. I was wondering if the following variance structure is possible to be obtained from the already implemented instances in varClasses, or I will need to code it myself.
>
> I want the variance of the response for observation i to be of the form v_i^2 = s^2 + s_i^2, where s^2 is a common for all observations unknown variability and s_i^2 is known, individual specific measurement error variance (can be 0).
>
> Thank you
> Best wishes
> Krzysztof Bartoszek

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

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

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

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 at 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

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 at protonmail.ch] 
Sent: Friday, 13 September, 2019 11:01
To: Viechtbauer, Wolfgang (SP)
Cc: r-sig-mixed-models at 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 ???????

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 at 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