Learn R Programming

onls (version 0.2)

logLik_o: Log-Likelihood for the orthogonal residuals

Description

Returns the log-likelihood as calculated from the orthogonal residuals obtained by residuals_o, including a Jacobian/normalizing-constant correction for the response and predictor precisions (weights, sigma_x, sigma_y) used by onls, so that the resulting value (and hence AIC/BIC) is comparable across onls fits made with different weighting schemes.

Usage

logLik_o(object)

Value

The log-likelihood as calculated from the orthogonal residuals, with attributes "df" (\(1 + q_{free}\), the number of free model parameters plus one for the residual scale) and "nobs"/"nall" (the number of observations).

Arguments

object

an object returned from onls.

Author

Andrej-Nikolai Spiess

Details

Let \(d_i\) be the weighted orthogonal distance of observation \(i\) at the fitted parameters (residuals_o, i.e. object$resid_o), and \(N\) the number of observations. A naive treatment of the \(d_i\) as i.i.d. \(N(0,\sigma^2)\) draws gives the usual concentrated Gaussian log-likelihood $$\ell_0 = -\frac{N}{2}\left(\log(2\pi) + 1 - \log(N) + \log\left(\sum_i d_i^2\right)\right),$$ which is what earlier versions of logLik_o returned. However, since onls already folds the response precision \(Qyy_i\) and predictor precision \(Qx_i\) into \(d_i\) itself (see 'Details' in onls), \(\ell_0\) alone omits the normalizing-constant (Jacobian) term that these precisions contribute to the underlying Gaussian density -- exactly as plain logLik/AIC on unweighted residuals from a weighted lm/nls fit would, if computed without R's own \(+\frac{1}{2}\sum_i \log(w_i)\) correction. logLik_o therefore adds the corresponding correction, $$\ell = \ell_0 \;+\; \frac{1}{2}\sum_{i=1}^{N} \log(Qyy_i) \;+\; \frac{1}{2}\sum_{i=1}^{N} \log\left|Qx_i\right|,$$ where \(|Qx_i|\) is the determinant of the (possibly per-observation) predictor precision matrix used by onls -- a scalar for single-predictor models, and the full \(p \times p\) determinant (correctly accounting for any predictor correlation) for multivariate models. For the default, fully unweighted onls fit, \(Qyy_i = 1\) and \(Qx_i = I_p\) for every \(i\), so both correction terms are exactly \(0\) and \(\ell = \ell_0\); the correction only changes the value for fits using weights, sigma_x, and/or sigma_y.

Note that \(d_i\) itself remains a simplification: it is the square root of a sum of \(p+1\) squared, precision-weighted Gaussian terms (one from the response, \(p\) from the predictors), not a single univariate normal draw, so treating \(\sum_i d_i^2 / N\) as the MLE of a common residual variance (as \(\ell_0\) does) is itself an approximation, inherited unchanged from the original (unweighted) formula. The correction above addresses only the missing weighting/precision Jacobian term, not this deeper simplification.

Examples

Run this code
DNase1 <- subset(DNase, Run == 1)
DNase1$density <- sapply(DNase1$density, function(x) rnorm(1, x, 0.1 * x))
mod1 <- onls(density ~ Asym/(1 + exp((xmid - log(conc))/scal)), 
             data = DNase1, start = list(Asym = 3, xmid = 0, scal = 1))
logLik_o(mod1)

## compare AIC of vertical versus orthogonal residuals.
AIC(mod1)
AIC(logLik_o(mod1))

## Comparing an unweighted and a weighted onls() fit on the same data:
## logLik_o() includes the precision Jacobian correction, so AIC/BIC
## from the two fits are validly comparable.
mod2 <- onls(density ~ Asym/(1 + exp((xmid - log(conc))/scal)), 
              data = DNase1, start = list(Asym = 3, xmid = 0, scal = 1),
              sigma_x = 0.05, sigma_y = 0.1)
AIC(logLik_o(mod1))
AIC(logLik_o(mod2))

Run the code above in your browser using DataLab