Re: [R] Joint confidence intervals for GLS models?

From: Frank E Harrell Jr <f.harrell_at_vanderbilt.edu>
Date: Thu 10 Aug 2006 - 00:54:08 EST

Kennedy David wrote:
> Dear All,
>
> I would like to be able to estimate confidence intervals for a linear
> combination of coefficients for a GLS model. I am familiar with John
> Foxton's helpful paper on Time Series Regression and Generalised Least
> Squares (GLS) and have learnt a bit about the gls function.
>
> I have downloaded the gmodels package so I can use the estimable
> function. The estimable function is very useful because it allows me to
> calculate confidence intervals for a linear combination of coefficients,
> but only for OLS models. For example, the code below calculates the
> confidence interval for the sum of the coefficient of petrol_A and the
> coefficient of petrol_B:
>

>> results <- lm(all_rural_count_capita ~ petrol_A + petrol_B +

> gdp_capita)
>> estimable(results,cm=c(0,1,1,0),conf.int=0.95)

>
> However, the estimable function does not seem to work for GLS objects,
> as shown below. The estimable documentation confirm that the object
> must be one of the following: lm, glm, lme, lmer.
>
>> results.gls <- gls(all_rural_count_capita ~ petrol_A + petrol_B +

> gdp_capita, correlation=corARMA(p=1),method='ML')
>> estimable(results.gls,cm=c(0,1,1,0),conf.int=0.95)

> Error in estimable(results.gls, cm = c(0, 1, 1, 0), conf.int = 0.95) :
> obj must be of class 'lm', 'glm', 'aov', 'lme', 'lmer', 'gee',
> 'geese' or 'nlme'
>
>
> Therefore, I am looking for a solution to this problem. I think that
> the solution (if it exists) may be down one of the following paths:
>
> 1) An alternative command which allows me to generate joint confidence
> intervals for the objects generated by the gls function. p.s. I note
> that the intervals function only appears to produce confidence intervals
> for each coeffcient (not for a linear combination of coeffcients).
>
> 2) An alternative means of generating GLS estimates as lm, glm, lme or
> lmer objects, so they can be inputed into the estimable function.
>
> Regards,
> David

I'm glad you are using gls because I think it's underused. When you don't need random effects but want to handle serial correlation, things are much easier with gls. The Design package (which requires the Hmisc package) has a function glsD that allows easy anova( ) and contrast( ) usage. Contrasts are simple because they are done in terms of differences in predicted values for any user settings:

contrast(fit, list(sex='male', times=1:5), list(sex='female', times=1:5))

Confidence intervals from contrast( ) (actually contrast.Design) are not simultaneous in the sense of providing simultaneous confidence bands over the whole time axis. I would welcome code for that using a chi-square multiple d.f. approximation or other approach.

A case study using glsD will be in a 2nd edition of my book Regression Modeling Strategies which is still some time away.

Frank

-- 
Frank E Harrell Jr   Professor and Chair           School of Medicine
                      Department of Biostatistics   Vanderbilt University

______________________________________________
R-help@stat.math.ethz.ch 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.
Received on Thu Aug 10 01:06:26 2006

Archive maintained by Robert King, hosted by the discipline of statistics at the University of Newcastle, Australia.
Archive generated by hypermail 2.1.8, at Thu 10 Aug 2006 - 02:20:25 EST.

Mailing list information is available at https://stat.ethz.ch/mailman/listinfo/r-help. Please read the posting guide before posting to the list.