# remlscore

From statmod v1.3.0
by Gordon Smyth

##### REML for Heteroscedastic Regression

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

- Keywords
- regression

##### 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)$ 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 $**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

##### 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)
```

*Documentation reproduced from package statmod, version 1.3.0, License: LGPL version 2 or newer*

### Community examples

Looks like there are no examples yet.