5
$\begingroup$

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 

Fits plot

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:

  1. 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?

  2. 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
$\endgroup$
  • 2
    $\begingroup$ You should only expect consistency within the same model/function even with identical fits, as likelihood is only defined up to a multiplicative constant (but it must be the same one on each item being compared). e.g. see Press, S. James (2003) Subjective and Objective Bayesian Statistics, 2E (see p35) and Box GEP and Tiao GC (1973), Bayesian Inference in Statistical Analysis (see p11, at least in the classics edition reprinted 1992, but I assume that's the same in the original) - different implementations of even the same model may yield different likelihoods. ... ctd $\endgroup$ – Glen_b Aug 16 '15 at 3:23
  • $\begingroup$ ctd ... so a big difference between likelihoods for almost identical model fits can happen even when everything is correct, if they don't keep/drop/use the same constant $\endgroup$ – Glen_b Aug 16 '15 at 3:26
  • 1
    $\begingroup$ This is a good (and basic) point, thanks. Basically it implies it's not possible to compare e.g. nls to nlrob models based on this criterion - unless the code is modified to use the same calculation/constant. $\endgroup$ – etov Aug 16 '15 at 15:09
2
$\begingroup$

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:

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

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

$\endgroup$
1
$\begingroup$

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

$\endgroup$
  • $\begingroup$ Thanks for your answer. However, you're not directly refering to the the robust log-likelihood calculation. So you seem to imply that the nlrob calculation is correct, while the nls calculation is faulty? $\endgroup$ – etov Nov 26 '14 at 8:02
  • $\begingroup$ Yes the loglikelihood calculation for nls is faulty in R. The following is the correct log-likelihood.. loglik3 = -N * (log(2 * pi) + 1 - log(N) - (1/N)*sum(log(w + zw)) + log(sum(w * res^2)))/2 $\endgroup$ – Arkajyoti Bhattacharya Dec 6 '14 at 20:19
0
$\begingroup$

I agree with Arkajyoti Bhattacharya.

See related discussion in:

Issue about 'logLik' command of R.

Perhaps the two threads should be merged.

$\endgroup$

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.