R help - Apr 2008 - function to generate weights for lm?

Hi,

I would like to use a weighted lm model to reduce heteroscendasticity. I am
wondering if the only way to generate the weights in R is through the
laborious process of trial and error by hand. Does anyone know if R has a
function that would automatically generate the weights need for lm?

Thanks,

-- 
Tom

	[[alternative HTML version deleted]]

On Tue, Apr 29, 2008 at 6:20 AM, tom soyer <tom.soyer at gmail.com>
wrote:
> Hi,
>
>  I would like to use a weighted lm model to reduce heteroscendasticity. I
am
>  wondering if the only way to generate the weights in R is through the
>  laborious process of trial and error by hand. Does anyone know if R has a
>  function that would automatically generate the weights need for lm?
Hi Tom,

The 'weights' argument to the 'gls' function in the nlme package
provides a great deal of flexibility in estimate weighting parameters
and model coefficients.  For example, if you want to model monotonic
heteroscedasticity by estimating the weights $E(Y)^{-2\alpha}$,
 you can use the varPower variance function class.  E.g., something like

f1 <- gls(y ~ x1 + x2, data = your.data, weights = varPower())

will estimate the regression coefficients and alpha parameter together
via maximum likelihood.  (note that the usual specification for varPower is
varPower(form = ~ your.formula), but by default the mean is used.  See
Ch 5 of the Pinheiro and Bates Mixed-effects Models book for details)

Kingsford Jones

Thanks Kingsford! I will cc r-help.

varPower() works for me. I want to use lm because predict.gls does not give
lwr and upr values, and also the bptest and ncv.test only work with lm
models...


On 4/29/08, Kingsford Jones <kingsfordjones@gmail.com>
wrote:
>
> hi Tom,
>
> Basically, a simple way to model non-constant variance when the
> variance increases (the common case) or decreases with the conditional
> mean of your response is to use "weights = varPower()", and gls
will
> estimate the parameter \alpha that defines the power relationship
> between the regression line and the error variances.  If for some
> reason you wanted to use lm rather then gls for your final model, then
> after estimating the gls model you could supply a weights argument to
> lm of the form:
>
> weights >
fitted(your.gls.model)^(-2*the.alpha.estimate.from.the.gls.model).
>
> and you should get the same (or very close) estimates as from the gls
> model but have a model of class 'lm'.
>
> There are many different ways to model non-constant variance and
> you'll need to choose one that is appropriate for your data.  If the
> model I described isn't appropriate then you should look at Ch 5 of
> P&B to learn about the other varFunc classes.
>
> good luck,
>
> Kingsford
>
> ps - would you mind forwarding to r-help in case this others have the
> same question.
>
>
> On Tue, Apr 29, 2008 at 3:24 PM, tom soyer <tom.soyer@gmail.com>
wrote:
> > Thanks Kingsford! What are the varClasses? I don't know how to use
> them....
> > If I choose varPower(), then does it mean the function would generate
> the
> > weights, w, so that w gives the most likely explanation of the
> relationship
> > between the variance and the independent variable?
> >
> >
> >
> >
> > On 4/29/08, Kingsford Jones <kingsfordjones@gmail.com> wrote:
> > > On Tue, Apr 29, 2008 at 6:20 AM, tom soyer
<tom.soyer@gmail.com>
> wrote:
> > > > Hi,
> > > >
> > > >  I would like to use a weighted lm model to reduce
> heteroscendasticity.
> > I am
> > > >  wondering if the only way to generate the weights in R is
through
> the
> > > >  laborious process of trial and error by hand. Does anyone
know if R
> has
> > a
> > > >  function that would automatically generate the weights need
for lm?
> > >
> > > Hi Tom,
> > >
> > > The 'weights' argument to the 'gls' function in
the nlme package
> > > provides a great deal of flexibility in estimate weighting
parameters
> > > and model coefficients.  For example, if you want to model
monotonic
> > > heteroscedasticity by estimating the weights $E(Y)^{-2\alpha}$,
> > > you can use the varPower variance function class.  E.g.,
something
> like
> > >
> > > f1 <- gls(y ~ x1 + x2, data = your.data, weights = varPower())
> > >
> > > will estimate the regression coefficients and alpha parameter
together
> > > via maximum likelihood.  (note that the usual specification for
> varPower
> > is
> > > varPower(form = ~ your.formula), but by default the mean is used.
See
> > > Ch 5 of the Pinheiro and Bates Mixed-effects Models book for
details)
> > >
> > > Kingsford Jones
> > >
> >
> >
> >
> > --
> > Tom
>


-- 
Tom

	[[alternative HTML version deleted]]