/ Home

DESCRIPTION

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

USAGE

remlscore(y, X, Z, trace=F, tol=1e-5, maxit=40)

 

REQUIRED ARGUMENTS

y		numeric vector of responses.
X		design matrix for predicting the mean.
Z		design matrix for predicting the variance.

 

OPTIONAL ARGUMENTS

trace		Logical variable. If true then output diagnostic information at each iteration.
tol		Convergence tolerance.
maxit		Maximum number of iterations allowed.

 

VALUE

beta		Vector of regression coefficients for predicting the mean.
se.beta		Standard errors for beta.
gamma		Vector of regression coefficients for predicting the variance.
se.gam		Standard errors for gamma.
mu		Estimated means.
phi		Estimated variances.
dev		Minus twice the REML log-likelihood.
h		Leverages.

 

DETAILS

Write m_i = E( y_i) for the expectation of the ith response and  s_i = Var( y_i). We assume the heteroscedastic regression model

m_i = x_i^Tbeta,    log(s_i) = z_i^Tgam,

where x_i and z_i are vectors of covariates, and beta and gam 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 (2000).

 

REFERENCES

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

EXAMPLES

> X
   Intercept B C 
 1         1 0 0
 2         1 0 1
 3         1 1 1
 4         1 1 0
 5         1 1 1
 6         1 1 0
 7         1 0 0
 8         1 0 1
 9         1 1 1
10         1 1 0
11         1 0 0
12         1 0 1
13         1 0 0
14         1 0 1
15         1 1 1
16         1 1 0
> Z
   Intercept C H I 
 1         1 0 0 0
 2         1 1 0 1
 3         1 1 1 0
 4         1 0 1 1
 5         1 1 0 1
 6         1 0 0 0
 7         1 0 1 1
 8         1 1 1 0
 9         1 1 0 1
10         1 0 0 0
11         1 0 1 1
12         1 1 1 0
13         1 0 0 0
14         1 1 0 1
15         1 1 1 0
16         1 0 1 1
> y
[1] 43.7 40.2 42.4 44.7 42.4 45.9 42.2 40.6 42.4 45.5 43.6 40.6 44.0 40.2
[15] 42.5 46.5
> library(Matrix)
> out <- remlscore(y,X,Z)
> cbind(Estimate=out$gamma,SE=out$se.gam)
             Estimate        SE 
Intercept -3.15886017 0.8313270
        C -2.73542576 0.8224828
        H -0.08588713 0.8351308
        I  3.33238821 0.8250499

Gordon Smyth. Copyright © 1996-2016. Last modified: 10 February 2004