4
$\begingroup$

In my thesis I use both 'lm' and 'nls'. I learned that the 'logLik' command of R produces different results even if the 'lm' and 'nls' objects represent the same model (see numerical example below).

If I understand correctly:

  • For 'lm' objects, the log-likelihood is calculated as:

enter image description here

  • For 'nls' objects, the log-likelihood is calculated as:

enter image description here

The differences between the two are the multiplication by n in the first term and the use of squared weights in the last term (both are marked in bold in the second equation).

The 'logLik' documentation that I found contains a reference for the log-likelihood of 'lm' objects. The derivation of this equation is clear to me.

Regarding the log-likelihood equation for 'nls' objects, I did not find a reference in the documentation. Could someone please provide guidance where can I find explanations about this equation, its derivation, and the reason for the difference from the 'lm' case?

Numerical Example

Data Generation:

> set.seed(1)
> E <- numeric()
> Y <- numeric()
> W <- numeric()
> X <- c(101:150)
> for(i in 1:50){
+   E <- c(E, rnorm(1, 0, sqrt(1+X[i])))
+ }
> Y <- 1+X+E
> for(i in 1:50){
+   W <- c(W, 1/(1+X[i]))
+ }

Model Fitting:

> L <- lm(Y ~ X, weights = W)
> summary(L)
Call:
lm(formula = Y ~ X, weights = W)

Weighted Residuals:
     Min       1Q   Median       3Q      Max
-2.32374 -0.46873  0.01821  0.62034  1.47708

Coefficients:
            Estimate Std. Error t value Pr(>|t|)
(Intercept)  2.64566   11.48112    0.23    0.819
X            0.99582    0.09209   10.81 1.84e-14 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.8399 on 48 degrees of freedom
Multiple R-squared:  0.709, Adjusted R-squared:  0.7029
F-statistic: 116.9 on 1 and 48 DF,  p-value: 1.84e-14


> NL <- nls(Y ~ A+B*X, start = list(A = 1, B = 1), weights = W)
> summary(NL)
Formula: Y ~ A + B * X

Parameters:
  Estimate Std. Error t value Pr(>|t|)
A  2.64566   11.48113    0.23    0.819
B  0.99582    0.09209   10.81 1.84e-14 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.8399 on 48 degrees of freedom

Number of iterations to convergence: 1
Achieved convergence tolerance: 1.068e-07

Log-Likelihood Value:

> L_Log_Lik1 <- logLik(L)
> L_Log_Lik2 <- 0.5*(sum(log(W))-50*log(2*pi)-50+50*log(50)
+                    -50*log(sum(W*(Y-predict(L))^2)))
> L_Log_Lik1
'log Lik.' -182.0457 (df=3)
> L_Log_Lik2
[1] -182.0457


> NL_Log_Lik1 <- logLik(NL)
> NL_Log_Lik2 <- 0.5*(50*sum(log(W))-50*log(2*pi)-50+50*log(50)
+                     -50*log(sum(W^2*(Y-predict(NL))^2)))
> NL_Log_Lik1
'log Lik.' -5983.21 (df=3)
> NL_Log_Lik2
[1] -5983.21

Best regards,

Pini

$\endgroup$

migrated from stackoverflow.com Aug 3 '16 at 23:32

This question came from our site for professional and enthusiast programmers.

  • $\begingroup$ A quick look at the code for stats:::logLik.nls suggests its formula differs in several ways from the one you quote and that in fact it is the same as for lm. $\endgroup$ – whuber Aug 7 '16 at 20:42
  • $\begingroup$ I mentioned in my previous post that there are differences between the two formulas. I looked at the code of 'logLik' command for 'lm' and 'nls' objects – the formulas above are identical to those that appear in the code. I formed the formulas that appear in the code because I wanted to show that they are almost identical (the differences are minor). $\endgroup$ – Pini Aug 16 '16 at 6:32
  • $\begingroup$ Furthermore, I ran simple example in R: I fitted simple linear model to data that I generated by using 'lm' and 'nls' commands. I got the same coefficients in the two models. Then, for each object, I calculated the log-likelihood twice – first time using the 'logLik' command and second time using the formula above directly. I got the same results. $\endgroup$ – Pini Aug 16 '16 at 6:35
  • $\begingroup$ But the fact remains that the two formulas you give are algebraically inequivalent. We have to conclude that the code you describe does not implement the formulas you have posted. However, the formulas can be made equivalent up to a constant difference (depending only on $n$ and the weights but not the data) provided you interpret "$w_i^2$" in the second formula as being "$w_i$" in the first. $\endgroup$ – whuber Aug 16 '16 at 13:07
  • 1
    $\begingroup$ I usually don't respond to efforts to carry on a public conversation in private, because that circumvents the structure and operation of this site and makes implicit demands on my time that I cannot meet. If you would like readers here to help, then you must edit your post to include the information they need to understand it correctly. $\endgroup$ – whuber Aug 18 '16 at 13:43
3
$\begingroup$

When I run Pini’s code, L$weights==NL$weights. So the weights and their interpretation appear to be the same. Predictions are also practically the same:

max(abs(NL$m$fitted()-L$fitted.values)) < 2e-13.

The key difference is in the residuals. resid() in nls objects appear to be weighted, while residuals in lm objects are not weighted (by itself this difference may or may not be a problem). Indeed,

max(abs(NL$m$resid()-L$residuals*L$weights^0.5)) < 2e-14.

In stats:::nlsModel.plinear, line 69 reads: dev <- sum(resid^2). This makes sense if resid are weighted. But stats:::logLik.nls, line 12 contains the term log(sum(w*res^2)). The same term appears in the exact same way in stats:::logLik.lm. As a result, the code in Pini’s example, using w^2 in the case of nls, does reproduce the values of logLik, for both lm and nls.

In addition, the difference between the logLik values of lm and nls does depend on the data. For example with set.seed(1) I get L_Log_Lik1-NL_Log_Lik1 = 5801.164; while with set.seed(2) I get L_Log_Lik1-NL_Log_Lik1 = 5800.496.

This is a confusing situation, and quite possibly we are missing something. Clarifications will be very highly appreciated.

$\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.