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?