Fit error-in-variables polynomial regression using mle2 (R)

Question

I need to fit a polynomial regression that accounts for measurement errors. I found out how to do it with a mcmc model (using RJags) and I would like to do it with a Maximum Likelihood Estimator (using mle2 function in R), since the model will be later more complex and mle2 will be faster than mcmc.

My model in RJags looks like this (I put some data to make the code reproducible):

modelFile = "model.txt"
modelString = "
model {
# Likelihood:
for (i in 1:N) {
y[i] ~ dnorm(y.hat[i], tauy[i])
y.hat[i] <- b[1] + b[2]*x.hat[i] + b[3]*z.hat[i] 

x.hat[i] ~ dnorm(x[i], taux[i])
z.hat[i] ~ dnorm(z[i], tauz[i])

taux[i] <- 1/pow(sdx[i],2)
tauy[i] <- 1/pow(sdy[i],2)
tauz[i] <- 1/pow(sdz[i],2)
}

for(j in 1:3) {b[j]~dunif(-2,2)}

}
"
writeLines(modelString,con=modelFile)


#Data

 ind <- data.frame(A = c(2.428, 2.601, 2.749, 2.553, 2.753, 2.421, 2.579, 2.415, 2.407, 2.509),
                  B = c(0.95, 0.99, 1.05, 1.00, 1.04, 0.96, 1.01, 0.95, 0.95, 1.01),
                  C = c(-0.04, -0.09, 0.01, 0.04, -0.15, 0.11, -0.17, -0.12, -0.13, 0.17),
                  eA=runif(10, 0, 0.5), eB=runif(10,0,0.5), eC=runif(10,0,0.2))

ml.data <- list(x=ind$A,
                    y=ind$B,
                z=ind$C,
                    sdy=ind$eA,
                sdx=ind$eB,
                    sdz=ind$eC,
                N=nrow(ind))

ml.par <- c("b")

ml.mod <- jags.model(modelFile,data=ml.data, n.chains=100, n.adapt=1000)

update(ml.mod, n.iter = 1000)

mcmc.out <- coda.samples(ml.mod, var=ml.par, n.iter=10000)

#summary of the posterior distributions of the parameters
summary(mcmc.out)

How can I translate this into an mle2, or how the function would be?

Answer 1

The model you are describing is a little bit strange... You have $i=1,...,N$ random variables $X_i$, $Y_i$, and $Z_i$'s, with each of them having their own mean $x_i$, $y_i$, $z_i$, and their own standard deviations $\tau_X$, $\tau_Y$, and $\tau_Z$ (in your case precision, but this is of less importance in here). There are three unknown parameters $\beta_0$, $\beta_1$, and $\beta_2$.

If you are familiar with R, than translating "likelihood" from your JAGS code to R would lead to

dnorm(y, mean = b[1] + b[2]*rnorm(n, x, sdx) + b[3]*rnorm(n, z, sdz), sd = sdy)

...are you sure you want $\hat x$ to be just some random number drawn from distribution with mean $x$ and standard deviation $\tau_X$..? (In most cases you rather estimate mean and don't have the data, rather than having mean and not having data.)

Notice also that since you know $\tau_{Yi}$ values, then you really do not have to care about $\tau_{Xi}$ and $\tau_{Zi}$ values since they are somehow already included in the overall variances (the $\tau_{Yi}$'s), a thing that you also know in advance!

However, since your "data" here is something that is random than each time you called the "likelihood" function it would be using different data (two rnorm's). This seems to be rather a simulation study and using a standard optimizer on data that is not static does not seem to be a sound idea in here - and won't be more efficient. Also this is not a maximum likelihood scenario because likelihood is probability of data given some parameter and you don't have data (you generate it at random).