0
$\begingroup$

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?

$\endgroup$
0
$\begingroup$

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).

| cite | improve this answer | |
$\endgroup$
  • $\begingroup$ Ok, I see the point. Actually my data correspond to physical measurements, not random. But anyway from your comments I see that in that case applying a mle would be more complicated and probably less efficient than directly working with the mcmc. I thought it would be the contrary, but now I see that not here. I need to learn more about all these methods...! but thanks a lot for your reply, it helps me. $\endgroup$ – LRD Nov 28 '15 at 14:00
  • $\begingroup$ Maybe... But what I am suggesting is that probably your model is ill-specified in JAGS. Are you sure you did not meant x[i] ~ dnorm(xhat[i], taux[i]) (i.e. "x comes from normal distribution with sd = taux and mean = xhat") ..? If you are conducting a simulation study this may be ok, but from your description is seems that you want to estimate something using some data, so it seems the model definition is wrong - maybe if it is corrected it could be estimated in mle2, nlme or lme4... $\endgroup$ – Tim Nov 28 '15 at 15:13

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.