lognlm {logNormReg}

Multiple linear regression with log Normal errors

Description

The function fits simple multiple linear regression models with log Normals erros. Two objectives as well as two optimizing functions can be used.

Usage

lognlm(formula, data, subset, weights, na.action, y = TRUE, start, model = TRUE, 
     lik = FALSE, opt = c("nlminb", "optim"), contrasts=NULL, ...)

Arguments

formula

a standard R formula with response and explanatory variables (and possible offset) specifying the regression model being fitted.

data

an optional data frame, list or environment containing some or all the variables in the model.

subset

an optional vector specifying a subset of observations to be used in the fitting process.

weights

an optional vector of (positive) weights to be used in the fitting process. Currently implemented only if lik=FALSE.

na.action

a function indicating what should happen when the data contain NAs. The default is set by the na.action setting of options.

y

logical. If TRUE the response vector is returned as y in the fit object.

start

(optional) starting values of the parameter to be estimated. If start is missing they are computed via ordinary least squares with the intercept \hat\beta_0 replaced by max(\hat\beta_0,median\{y_i\}).
If lik=TRUE (i.e. a log Normal model is fitted), start refers to the regression parameters and the error standard deviation; if lik=FALSE, start does not include the starting guess for the standard deviation.

model

logical. If TRUE the model frame is returned as model in the fit object.

lik

If TRUE the log Normal log likelihood is optimized, otherwise the sum of squared residuals based on the logs (see Details).

opt

the optimization function to be used. nlminb appears to be more efficient, probably because it uses (unlike optim) also the hessian matrix (supplied in the code). However results are often indistinguishable.

contrasts

an optional list. See the contrasts.arg of model.matrix.default.

...

optional arguments passed on to the optimizing functions (nlminb or optim), (and therefore should be consistent with that).

Details

lognlm fits simple linear regression models with log Normal errors and identity link. Actually two objectives could be used.

If lik=TRUE the usual log Normal likelihood is optimized, otherwise estimation is based on minimization of the following loss function

\sum_i(\log y_i-\log \mu_i )^2

where \mu_i=x_i^T\beta is the mean function equal to the linear predictor (as an identity link is exploited).

Value

A list with components

coefficients

the regression parameters estimate.

loglik

The objective function value, namely the log Normal log likelihood or the sum of the squared ‘log residuals’ (depending on lik option).

s2

the error variance estimate.

fitted.values

the fitted values.

residuals

the raw residuals on the original scale, i.e. 'observed - fitted'.

grad

the gradient at solution.

hessian

the hessian matrix at solution.

Ehessian

the expected hessian matrix at solution (only if lik=FALSE).

convergence

the convergence code coming from the fitter function. 0 means succefull convergence.

call

the matched call.

y

the response vector (provided that y=TRUE has been set).

opt

the employed optimizer.

lik

logical, indicating if the fit comes from a log Normal likelihood approach.

xlevels

(only where relevant) a record of the levels of the factors used in fitting.

terms

the terms object used.

contrasts

(only where relevant) the contrasts used.

model

if requested, i.e. model=TRUE has been set (the default), the model frame used.

offset

the (possible) offset used.

Author(s)

Vito M.R. Muggeo

See Also

See also print.lognlm and summary.lognlm to display results.

Examples

n=300
s=.4
set.seed(123)      #just to get reproducible results..

x<-seq(.1,10,l=n) #covariate
mu<- 10+2*x  #linear regression function
y<-rlnorm(n, log(mu)-s^2/2, s) #data..

o0<-lm(log(y)~x) #the usual but WRONG model
o<- lognlm(y~x, lik=TRUE) #fit the 'right' model by ML

plot(x,y)
lines(x, mu, lwd=2)
points(x, exp(fitted(o0)), col=2, type="l", lwd=2)
points(x, fitted(o), col=3, type="l", lwd=2)
legend("topleft", legend=c("true", "lm(log(y)~x)", "lognlm(y~x)"), 
    col=c(1,2,3), lwd=2)

#Sometimes people would estimate parameters by minimizing a least square objective 
# (i.e. by setting 'lik=FALSE', see Details), wherein data would come from 
# Y = mu * exp(eps) where eps~N(0,s)..
y1<-mu*exp(rnorm(n,0,1)) #data..
o1<-lognlm(y1~x, lik=FALSE) #set 'lik=FALSE', see Details