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

[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 ‐‐‐‐‐‐‐
On Wednesday, September 11, 2019 7:12 PM, Viechtbauer, Wolfgang (SP) <wolfgang.viechtbauer using maastrichtuniversity.nl> wrote:

> 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 using 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 using 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 using maastrichtuniversity.nl>