R/lmInt.R

In tsieger/tsiMisc: Miscelaneous utility functions, mostly related to exploratory statistical analysis

lmInt<-structure(
function # Confidence and prediction intervals for regression.
##description<<
## \code{\link{lmInt}} contructs confidence and prediction intervals for
## linear regression. It computes confidence intervals around the regression
## line (i.e. the point-wise confidence bands of \eqn{E(Y | X = x)}
## for each individual \eqn{x}), confidence intervals for the
## regression line (i.e. the simultaneous confidence bands of
## \eqn{E(Y)} for all \eqn{x}), and prediction intervals (i.e.
## point-wise confidence bands for new observations \eqn{Y} for each
## individual \eqn{x}).
##
##details<<
## The confidence intervals around the regression line and the
## prediction intervals are computed using the
## \code{\link[stats]{predict.lm}} function.
## The confidence intervals for the regression line are computed
## according to eq. (4.15) in Zvara2008.
##
##references<< Karel Zv\'{a}ra: Regrese, Matfyzpress Praha 2008
##
##seealso<< \code{\link{plot.lmInt}}, \code{\link[stats]{predict.lm}},
## \code{\link[stats]{lm}}
##
(m, ##<< a fitted linear model
newdata = NULL, ##<< a data frame to look for predictors. If omitted,
## fitted values will be used. If a single numeric and if the model
## specifies a single predictor only, a new data frame will get
## generated by taking \code{newdata} values to uniformly span the range
## of the predictor.
level = .95 ##<< confidence level
) {
  if (!inherits(m,'lm')) {
    stop('need an object of class \'lm\'')
  }
  if (is.null(m$model)) {
    stop('the first argument of type \'lm\' has no \'model\' set')
  }
  #tt<-terms(m)
  #hasintercept<-attr(tt,"intercept")>0L

  preds<-delete.response(terms(m))
  # the number of predictors (excluding offset terms)
  k<-length(as.list(attr(preds,'variables')))-1-length(attr(preds,'offset'))
  if (k!=ncol(m$model)-1-length(attr(preds,'offset'))) {
    .pn(k)
    .pn(ncol(m$model))
    .pn(m$model)
    stop('FIXME: unrecognized model')
  }

  # construct newdata if requested (and possible)
  if (!missing(newdata)) {
    if (!inherits(newdata,'data.frame')) {
      if (!is.numeric(newdata) || length(newdata)!=1) {
        stop('invalid \'newdata\' argument')
      }
      if (k>1) {
        stop('don\'t know how to generate \'newdata\' for several predictors')
      }
      newdata<-data.frame(x=seq(min(m$model[,2]),max(m$model[,2]),len=newdata))
      colnames(newdata)[1]<-colnames(m$model)[2]
    }
  }

  # compute predictions intervals
  if (!missing(newdata)) {
    confintForLm<-predict(m,newdata,interval="confidence",level=level)
    pred<-predict(m,newdata,se.fit=T,level=level)
    predint<-predict(m,newdata,interval="predict",level=level)
  } else {
    confintForLm<-predict(m,interval="confidence",level=level)
    pred<-predict(m,se.fit=T,level=level)
    predint<-predict(m,interval="predict",level=level)
  }

  li<-data.frame(
    fit=pred$fit,
    ciaLwr=confintForLm[,'lwr'],
    ciaUpr=confintForLm[,'upr'],
    cifLwr=pred$fit-pred$se.fit*sqrt((k+1)*qf(level,k+1,pred$df)),
    cifUpr=pred$fit+pred$se.fit*sqrt((k+1)*qf(level,k+1,pred$df)),
    piLwr=predint[,'lwr'],
    piUpr=predint[,'upr'])

  if (!missing(newdata)) {
    if ('(Intercept)'%in%colnames(newdata)) {
      newdata<-newdata[,!colnames(newdata)%in%'(Intercept)']
    }
    li<-cbind(newdata,li)
  } else {
    mm<-model.matrix(m)
    mm<-mm[,!colnames(mm)%in%'(Intercept)',drop=FALSE]
    li<-cbind(mm,li)
  }

  class(li)<-c('lmInt',class(li))

  return(li)
  ### An object of class \code{lmInt} - a data frame consisting of the model
  ### matrix followed by columns of the:
  ### \code{fit} holding the mean value fitted by the regression model,
  ### \code{ciaLwr} and \code{ciaUpr} holding confidence intervals
  ### \emph{a}round the regression line,
  ### \code{cifLwr} and \code{cifUpr} holding confidence intervals
  ### \emph{f}or the regression line, and
  ### \code{piLwr} and \code{piUpr} holding predictions intervals.
},ex=function() {
  iris.setosa<-iris[iris$Species=='setosa',]
  attach(iris.setosa)

  # simple line fitting (one predictor only)
  m <- lm(Sepal.Width ~ Sepal.Length)
  lmi <- lmInt(m)
  plot(Sepal.Length, Sepal.Width)
  plot(lmi, fit = TRUE, lty = 1, col='red')
  plot(lmi, cia = TRUE, lty = 1)
  plot(lmi, cif = TRUE, lty = 2)
  plot(lmi, pi = TRUE, lty = 3)
  legend('bottomright', bg='white',
    c('fitted regression line',
    'confidence interval around the regression line',
    'confidence interval for the regression line',
    'prediction int.'),
    col = c('red', 'black', 'black', 'black'), lty = c(1, 1, 2, 3))

  # using more predictors
  m2 <- lm(Sepal.Width ~ Sepal.Length + Petal.Length + offset(Petal.Width))
  lmi2 <- lmInt(m2, newdata = data.frame(
    Sepal.Length = Sepal.Length,
    Petal.Length = mean(Petal.Length),
    Petal.Width = mean(Petal.Width)))
  plot(Sepal.Length, Sepal.Width)
  plot(lmi2, fit = TRUE, lty = 1, col='red')
  plot(lmi2, cia = TRUE, band = TRUE, border = NA, col = scales::alpha('red',.3))
  plot(lmi2, cif = TRUE, band = TRUE, border = NA, col = scales::alpha('blue',.2))
  plot(lmi2, pi = TRUE, band = TRUE, border = NA, col = scales::alpha('black',.2))
  legend('bottomright', bg='white',
    c('fitted regression line',
    'confidence interval around the regression line',
    'confidence interval for the regression line',
    'prediction int.'),
    col = c('red', scales::alpha('red',.3), scales::alpha('blue',.2), scales::alpha('black',.1)),
    lty = 1, lwd = c(1, 15, 15, 15))

  detach(iris.setosa)
})