remlscore

0th

Percentile

REML for Heteroscedastic Regression

Fits a heteroscedastic regression model using residual maximum likelihood (REML).

Keywords
regression
Usage
remlscore(y, X, Z, trace=FALSE, tol=1e-5, maxit=40)
Arguments
y
numeric vector of responses
X
design matrix for predicting the mean
Z
design matrix for predicting the variance
trace
Logical variable. If true then output diagnostic information at each iteration.
tol
Convergence tolerance
maxit
Maximum number of iterations allowed
Details

Write $\mu_i=E(y_i)$ for the expectation of the $i$th response and $s_i=\var(y_i)$. We assume the heteroscedastic regression model $$\mu_i=\bold{x}_i^T\bold{\beta}$$ $$\log(\sigma^2_i)=\bold{z}_i^T\bold{\gamma},$$ where $\bold{x}_i$ and $\bold{z}_i$ are vectors of covariates, and $\bold{\beta}$ and $\bold{\gamma}$ are vectors of regression coefficients affecting the mean and variance respectively. Parameters are estimated by maximizing the REML likelihood using REML scoring as described in Smyth (2002).

Value

  • List with the following components:
  • betavector of regression coefficients for predicting the mean
  • se.betavector of standard errors for beta
  • gammavector of regression coefficients for predicting the variance
  • se.gamvector of standard errors for gamma
  • muestimated means
  • phiestimated variances
  • devianceminus twice the REML log-likelihood
  • hnumeric vector of leverages
  • cov.betaestimated covariance matrix for beta
  • cov.gamestimated covarate matrix for gamma
  • iternumber of iterations used

References

Smyth, G. K. (2002). An efficient algorithm for REML in heteroscedastic regression. Journal of Computational and Graphical Statistics 11, 836-847.

Aliases
  • remlscore
Examples
data(welding)
attach(welding)
y <- Strength
# Reproduce results from Table 1 of Smyth (2002)
X <- cbind(1,(Drying+1)/2,(Material+1)/2)
colnames(X) <- c("1","B","C")
Z <- cbind(1,(Material+1)/2,(Method+1)/2,(Preheating+1)/2)
colnames(Z) <- c("1","C","H","I")
out <- remlscore(y,X,Z)
cbind(Estimate=out$gamma,SE=out$se.gam)
Documentation reproduced from package statmod, version 1.4.9, License: LGPL (>= 2)

Community examples

Looks like there are no examples yet.