remlscor: REML for Heteroscedastic Regression

In statmod: Statistical Modeling

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 mu_i = E(y_i) and sigma_i^2 = var(y_i) for the expectation and variance of the i'th response. We assume the heteroscedastic regression model

mu_i = x_i^T beta

log(sigma_i^2 = z_i^T gamma ,

where x_i and z_i are vectors of covariates, and beta 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
iter	number of iterations used

Author(s)

Gordon Smyth

References

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

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)

statmod documentation built on Jan. 6, 2023, 5:14 p.m.