remlscore {statmod}R Documentation

REML for Heteroscedastic Regression

Description

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

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 μ_i=E(y_i) for the expectation of the $i$th response and $s_i=var(y_i)$. We assume the heteroscedastic regression model

μ_i=x_i^Tβ

log(σ^2_i)=z_i^Tgamma,

where $x_i$ and $z_i$ are vectors of covariates, and $β$ and $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:

beta vector of regression coefficients for predicting the mean
se.beta vector of standard errors for beta
gamma vector of regression coefficients for predicting the variance
se.gam vector of standard errors for gamma
mu estimated means
phi estimated variances
deviance minus twice the REML log-likelihood
h numeric vector of leverages
cov.beta estimated covariance matrix for beta
cov.gam estimated covarate matrix for gamma

Author(s)

Gordon Smyth

References

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

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)

[Package statmod version 1.2.4 Index]