statmod (version 1.4.9)

remlscore: 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

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

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).

References

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

Examples

Run this code
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)

Run the code above in your browser using DataCamp Workspace