Become an expert in R — Interactive courses, Cheat Sheets, certificates and more!
Get Started for Free

remlscor

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 ith response and s_i=\var(y_i). We assume the heteroscedastic regression model

μ_i=\bold{x}_i^T\bold{β}

\log(σ^2_i)=\bold{z}_i^T\bold{γ},

where \bold{x}_i and \bold{z}_i are vectors of covariates, and \bold{β} and \bold{γ} 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.

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

Statistical Modeling

v1.4.36
GPL-2 | GPL-3
Authors
Gordon Smyth [cre, aut], Yifang Hu [ctb], Peter Dunn [ctb], Belinda Phipson [ctb], Yunshun Chen [ctb]
Initial release
2021-05-10

We don't support your browser anymore

Please choose more modern alternatives, such as Google Chrome or Mozilla Firefox.