[R] question regarding variance function in gls

2010-03-15 Thread 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 

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.

[R] varComb in gls/lme

2010-03-08 Thread 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 

 

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

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.

[R] 2 different colors for 2 groups of lines using xyplot

2009-02-11 Thread Yu, Xuesong 
Hi All,
 
 
I am trying to use xyplot to plot several lines in one panel with red
color  for a group of lines, say 3 lines and blue color for another
group of lines, say 4. I know this can be easily  done using regular
plot function. But i could not figure out how to do this in xyplot. 
 
 
Any help would be greatly appreciated. 
 
Xuesong 
 
 
   


[[alternative HTML version deleted]]

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.