Log-likelihood (and AIC) of robust nlrob model differs from standard nls model 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 


However:
 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
[1] 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

 A: There are some errors in the source code of logLik.nls.
(1/N) should be multiplied with sum(log(weights)) for the calculation of val in the source code. On the other hand, "res" is taken as weighted residuals and even after that weights are multiplied by residuals for calculating "val".
That calculation mistake is causing the difference in AIC values for nls model with weights.
I would expect that the R developers would recover this problem. I have written a sample code for pointing out this problem.
You can see this by the following code:
methods(logLik)   
stats:::logLik.nls 

set.seed(12345)
n=2000
npar=3 # alpha, beta & sigma
yt<-rnorm(n, mean=1, sd=0.25)
xt<-rnorm(n, mean=0, sd=0.16)
wt = runif(n)

data<-cbind(yt,xt)
modt=list()
modt[[1]]<-nls(yt~alpha+beta*xt, start=list(alpha=1, beta=1), weights = wt)

i=1
sse<-deviance(modt[[i]])

######################
#residuals used in the source code
######################
res1 = modt[[1]]$m$resid()
######################
#Our calculated Residuals without weights
######################
coef = coef(modt[[1]])

yhat = coef[[1]] +coef[[2]]*xt

res2 = yt - yhat
######################
#Residuals obtained by "residuals" command
######################

res3 = residuals(modt[[1]])

######################
#Matching between the residuals
#If the result is 2000 then all are matching
#Otherwise all are not same
######################

sum(res2==res3)

sum(res1==res2)

#######################
#matching the calculated residuals mulitplied by sqrt(wt) with the residuals obtained by the source code formula
#######################

sum((sqrt(wt)*res2)==res1)


#######################
#Calculating Loglikelihood
#######################

sigma = sqrt(sse/n)

########################
#Our likelihood
########################
loglik1<-(-(n/2)*log(2*pi)-n*log(sigma) + sum(log(wt))/2-sse*(1/(2*sigma^2)))

########################
#Source code likelihood
########################

#Using Source code residuals

res <- res1
N <- length(res)
w = wt
zw <- w == 0

loglik2 = -N * (log(2 * pi) + 1 - log(N) - sum(log(w + zw)) + 
        log(sum(w * res^2)))/2


########################
#Corrected Likelihood
########################

#Using Our Calculated Residuals
res <- res2
N <- length(res)
w = wt
zw <- w == 0

loglik3 = -N * (log(2 * pi) + 1 - log(N) - (1/N)*sum(log(w + zw)) + 
                  log(sum(w * res^2)))/2

(1/N) should be multiplied with sum(log(w+zw)). It is correctly written for logLik.lm. you can check its source code in the same manner as mentioned above
print(loglik1)#Our Calculation

print(loglik2)#source code

print(loglik3)#corrected

logLik(modt[[i]])#Using the command "logLik"

aic1<-(-2*loglik1+2*npar)
aic2<-(-2*loglik2+2*npar)
aic3<-(-2*loglik3+2*npar)

print(aic1)#Our Calculation

print(aic2)#source code

print(aic3)#corrected

AIC(modt[[i]])#Using the command AIC





#########################################################################
#Results
#########################################################################

######################
> #Matching between the residuals
  > #If the result is 2000 then all are matching
  > #Otherwise all are not same
  > ######################
> 
  > sum(res2==res3)
[1] 2000

Our Calculated residuals match with the outcome that "residuals" command produce
> 
  > sum(res1==res2)
[1] 0

Our Calculated residuals do not match with the residuals used for logLik.nls source code
> 
  > #######################
> #matching the calculated residuals mulitplied by sqrt(wt) with the residuals obtained by the source code formula
  > #######################
> 
  > sum((sqrt(wt)*res2)==res1)
[1] 2000

Our residuals multiplied with sqrt(wt) match with the residuals used for logLik.nls source code

   #Our Calculation

print(loglik1)
      [1] -372.4776
print(loglik2) #source code
      [1] -2010632
print(loglik3) #corrected
      [1] -372.4776
logLik(modt[[i]]) #Using the command "logLik"
      'log Lik.' -2010632 (df=3)
aic1<-(-2*loglik1+2*npar)
  aic2<-(-2*loglik2+2*npar)
  aic3<-(-2*loglik3+2*npar)
print(aic1) #Our Calculation
      [1] 750.9552
print(aic2) #source code
      [1] 4021270
print(aic3) #corrected
      [1] 750.9552
AIC(modt[[i]]) #Using the command AIC
      [1] 4021270  

A: This is just a partial answer, but since the question didn't get any attention until now, I figured it's better than nothing.
As noted in the question, the log-likelihood difference stems from the difference in the calculated residuals.
It turns out that calling residuals(fit.nlrob) ends up running a different code than fit.nlrob\$m$resid():
The class of fit.nlrob\$m is "nlsModel". 
nlsModel$resid() uses weights for calculating the residulas:
(from nlsModel):
 .swts <- if (!missing(wts) && length(wts)) 
        sqrt(wts)
    #...
    resid <- .swts * (lhs - rhs)
    #...
    resid = function() resid

However, residuals(fit.nlrob) calls residuals.nls (can be verified by using trace()), which ignores weights:
(from stats:::residuals.nls):
val <- as.vector(object$m$lhs() - object$m$fitted())

So, in answer to the questions above:


*

*Apparently the difference in log-likelihood is just a technical side-effect. It seems to me that it can be calculated both using weights or without weights; but for the log-likelihoods to be comparable, they should probably be calculated on the same basis.

*Consequently, you should probably use your own version of logLik (or any other method) to make sure the calculation is done on the same basis.
A: I agree with Arkajyoti Bhattacharya.
See related discussion in:
Issue about 'logLik' command of R.
Perhaps the two threads should be merged.
