remlscore {statmod} | R Documentation |

## 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 $i$th response and $s_i=var(y_i)$.
We assume the heteroscedastic regression model

*μ_i=***x**_i^T**β**

*log(σ^2_i)=***z**_i^T**gamma**,

where $**x**_i$ and $**z**_i$ are vectors of covariates, and $**β**$ 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 |

### Author(s)

Gordon Smyth

### References

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

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

[Package

*statmod* version 1.2.4

Index]