REML for Heteroscedastic Regression
Fits a heteroscedastic regression model using residual maximum likelihood (REML).
remlscore(y, X, Z, trace=FALSE, tol=1e-5, maxit=40)
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 |
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).
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 |
Gordon Smyth
Smyth, G. K. (2002). An efficient algorithm for REML in heteroscedastic regression. Journal of Computational and Graphical Statistics 11, 836-847.
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)Please choose more modern alternatives, such as Google Chrome or Mozilla Firefox.