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}}
