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