I am attempting to use R for model selection based on the AIC statistic. When comparing linear models with or without weighting, my code in R informs me that weighting is preferable compared to no-weighting, and these results are confirmed in other software (GraphPad Prism). I have sample code using real data from a standard curve:
#Linear Curve Fitting
a <- c(0.137, 0.412, 1.23, 3.7, 11.1 ,33.3)
b <- c(0.00198, 0.00359, 0.00816, 0.0220, 0.0582, 0.184)
m1 <- lm(b ~ poly(a,1))
m2 <- lm(b ~ poly(a,1), weight=1/a)
n1 <- 6 #Number of observations
k1 <- 2 #Number of parameters
When I calculate AIC using either the internal function in R or via manual calculation in which:
AIC = n + n log 2π + n log(RSS/n) + 2(k + 1) with n observations and k parameters
I get equivalent AIC values for the non-weighted model. When I analyze the effect of weighting, the manual AIC value is lower, however the end result is that both the internal and manual AIC suggest that weighting is preferred.
> AIC(m1); n1+(n1*log(2*pi))+n1*(log(deviance(m1)/n1))+(2*(k1+1))
[1] -54.83171
[1] -54.83171
> AIC(m2); n1+(n1*log(2*pi))+n1*(log(deviance(m2)/n1))+(2*(k1+1))
[1] -64.57691
[1] -69.13025
When I try the same analysis using a nonlinear model, the difference in AIC between the internal function and manual calculation is more profound. Below is a code of examplar Michaelis-Menten kinetic data:
c <- c(0.5, 1, 5, 10, 30, 100, 300)
d <- c(3, 5, 20, 50, 75, 200, 250)
m3 <- nls(d ~ (V * c)/(K + c), start=list(V=10, K=1))
m4 <- nls(d ~ (V * c)/(K + c), start=list(V=10, K=1), weight=1/d^2)
n2 <- 7
k2 <- 2
The AIC are calculated as indicated for the first two models:
> AIC(m3); n2+(n2*log(2*pi))+n2*(log(deviance(m3)/n2))+(2*(k2+1))
[1] 58.48839
[1] 58.48839
> AIC(m4); n2+(n2*log(2*pi))+n2*(log(deviance(m4)/n2))+(2*(k2+1))
[1] 320.7105
[1] 0.1538546
Similar to the linear example, the internal AIC and manual AIC values are the same when data are not weighted (m3). The problem occurs with weighting (m4) as the manual AIC estimate is much lower. This situation is similar to what was asked in a related problem AIC with weighted nonlinear regression (nls).
I earlier mentioned GraphPad Prism, which for both the models and datasets given above showed lower AICs when weighting was used. My question then is why is there such a difference in the internal vs. manual AIC estimates in R when weighting the data (for which the outcome is different for nonlinear model compared to a linear one)? Ultimately, should I regard the internal AIC value or the manual value as being more correct, or am I using a wrong equation?
The discrepancy you are seeing is from using the unweighted log-likelihood formula in the manual calculations for a weighted model. For example, you can replicate the AIC results for m2 and m4 with the following adjustments:
In the case of m2, you simply need to subract sum(log(m2$weights)) from your calculation:
AIC(m2); n1+(n1*log(2*pi))+n1*(log(deviance(m2)/n1))+(2*(k1+1)) - sum(log(m2$weights))
[1] -64.57691
[1] -64.57691
In the case of m4, you would have to swap the deviance call with a weighted residuals calculation, and subtract n2 * sum(log(m4$weights)) from your results:
AIC(m4); n2+(n2*log(2*pi))+n2*(log(sum(m4$weights * m4$m$resid()^2)/n2))+(2*(k2+1)) - n2 * sum(log(m4$weights))
[1] 320.7105
[1] 320.7105
I believe the derivation for the formula used by logLikin m2 is pretty straight forward and correct, but I am not as sure about m4. From reading some other threads about logLik.nls() (example 1, example 2), it seems like there is some confusion about the correct approach for the nls estimate. To summarize, I believe AIC is correct for m2; I was not able to verify the math for the weighted nls model and would lean towards using the m2 formula again in that case (but replace deviance calculation with weighted residuals), or (maybe better) not use AIC for the nls model