[R] Joint confidence intervals for GLS models?

[Frank E Harrell Jr]

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