Comparing models generated by nlrob to ones generated by nls, I've noticed that even though the models might be nearly identical, the log-likelihood of the models is sometimes significantly different, e.g. (reproducible code below):
coefficients(fit.nls) a b c 1.990388 -3.049477 1.019370 coefficients(fit.nlrob) a b c 1.990650 -3.049251 1.019687
logLik(fit.nls) 'log Lik.' 3.341953 (df=4) logLik(fit.nlrob) 'log Lik.' 5.114307 (df=4)
Looking into the Loglik.nls code, this seems to arise from using the following line for obtaining the model's residuals:
res <- object$m$resid()
Using the same code for calculating log-likelihood but with the residuals() method instead, the obtained log-likelihood is much closer to the one generated for the nls model:
logLik.nlrob.method  3.340076
the nlrob documentation does relate to the different residuals, however I'm not clear if this has to do with the above:
residuals(.), by default type = "response", returns the residuals e_i, defined above as e[i] = Y[i] - f(x[i], theta^). These differ from the standardized or weighted residuals which, e.g., are assumed to be normally distributed, and a version of which is returned in working.residuals component.
I'm not familiar with robust statistics, so I wonder:
Does the different log likelihood indeed reflects some deeper difference between the models? or is it just a technical difference arising from non-standard use of the resid() function?
Let's say I'd like to choose between nls and nlrob models based on the AIC criterion. for the sake of the question, I'll assume I don't know enough about the underlying model to say anything about the likelihood of outliers. Should I use the standard AIC/logLik functions for that? or should I use my own version, based on the residuals() function instead of the object\$m\$resid()?
code (intended for interactive session):
library(robustbase) library(ggplot2) #generate data set.seed(1357) f <- function(x) 2*x*x-3*x+1 x=rnorm(50,0,2) df <- data.frame(x=x,y=f(x)+rnorm(length(x),0,0.3)) #fit fit.nls <- nls(y~a*I(x*x)+b*x+c,start=c(a=2,b=-3,c=1),data=df) fit.nlrob <- nlrob(y~a*I(x*x)+b*x+c,start=c(a=2,b=-3,c=1),data=df) coefficients(fit.nls) coefficients(fit.nlrob) #plot s <- data.frame(x=seq(0,2,by=0.1)) ggplot(data=s,aes(x=x)) + geom_line(aes(y=f(x),color="f(x)")) + geom_line(aes(y=predict(fit.nls,s),color="nls")) + geom_line(aes(y=predict(fit.nlrob,s),color="nlrob")) #log likelihood logLik(fit.nls) logLik(fit.nlrob) #Let's compare the two versions of the residuals res.obj <- fit.nlrob$m$resid() res.method <- residuals(fit.nlrob) N <- length(res.obj) #all weights are 1 in both models, so the weight terms are cancelled logLik.nlrob.obj <- -N * (log(2 * pi) + 1 - log(N) + log(sum(res.obj^2)))/2 logLik.nlrob.method <- -N * (log(2 * pi) + 1 - log(N) + log(sum(res.method^2)))/2 #this one is identical to the original logLik.nls function logLik.nlrob.obj #and this one is very close to the fit.nls log likelihood logLik.nlrob.method