# [R] question regarding variance function in gls

Yu, Xuesong xyu at scharp.org
Mon Mar 15 19:00:55 CET 2010

```Dear R-help members,

I have a question regarding how to use varComb function to specify a
variance function for the "weights"  in the gls.  I  need to fit a
linear model with heteroscedasticity. The variance function is
exp(c0+nu0*W +nu1*W^2) where W is  a covariate.  Initially I want to use
varFunc to define my own variance function following the instruction in
the Pinheiro and Bates (2000), but I could not make it work. Then I used
varComb in gls with  weights=varComb(varExp(form=~W),
varExp(form=~I(W^2))))).  But the estimated variance parameters seems to
have a large discrepancy from the true values
(I used the simulated data).  This makes me wonder if  it is a right way
to model variance function  exp(c0+nu0*W +nu1*W^2) using "varComb".  The
codes and outputs are copied below.

Any suggestions and help are very apprecited

library(nlme)

simulate.pilot = function(m, mn, sigma, alpha0, alpha1, c0, nu0, nu1) {

pilot.dat=data.frame(W=rnorm(m, mean=mn, sd=sigma))
pilot.dat=transform(pilot.dat, Y=rnorm(m, mean=alpha0 + alpha1*W,
sd=sqrt(exp(c0+nu0*W+nu1*W^2))))

pilot.dat
}

mn=3.3
sigma=sqrt(0.5)

alpha0=0.1
alpha1=3

m=200
n=200

c0=-2.413;   nu0=-0.2; nu1=0.3

simu.dat=simulate.pilot(m, mn, sigma, alpha0, alpha1, c0, nu0, nu1)

fit1=try(gls(Y~W, data=simu.dat, weights=varComb(varExp(form=~W),
varExp(form=~I(W^2)))))
c0.hat= log(fit1\$sigma^2)
nu0.hat=2*fit1\$modelStruct\$varStruct\$A
nu1.hat=2*fit1\$modelStruct\$varStruct\$B

> c0.hat
 -1.570104
> nu0.hat
 -0.787264
> nu1.hat
 0.4057129
>

Thanks

Xuesong

```