[R] question regarding variance function in gls

[Yu, Xuesong]

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