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[1]

 nu1.hat=2*fit1$modelStruct$varStruct$B[1]

 

> c0.hat

[1] -1.570104

> nu0.hat

[1] -0.787264

> nu1.hat

[1] 0.4057129

> 

 

Thanks

 

Xuesong 

 

CONFIDENTIALITY NOTICE: This e-mail message, including a...{{dropped:12}}

