Dear Colleagues:
I am experiencing difficulties to extract standard deviations
of slope and interecept in "gls"
Consider these two functions:
The first one fits linear model, prints out summary and returns standard
errors of coefficients
___________________________________________________________________________
function(N = 500, alf = 0, bet = 1)
{
#now.lm
x <- rchisq(N, df = 3)
eps <- rexp(N)
y <- alf + bet * x + eps
LM <- lm(y ~ x, data.frame(cbind(x, y)))
suLM <- summary(LM)
cfLM <- suLM$coefficients
er.cfLM <- cfLM[, 2]
print(suLM)
list(er.coef = er.cfLM)
}
The second one fits generalized linear model, prints out summary and
returns standard errors of coefficients
________________________________________________________________________
function(N = 500, alf = 0, bet = 1)
{
#now.gls
x <- rchisq(N, df = 3)
eps <- rexp(N)
y <- alf + bet * x + eps
GL <- gls(y ~ x, data.frame(cbind(x, y)), weights = varPower())
suGL <- summary(GL)
print(suGL)
cfGL <- suGL$coefficients
er.cfGL <- cfGL[, 2]
list(er.coef = er.cfGL)
}
___________________________________________________________________
here is how they work:
> now.lm()
Call: lm(formula = y ~ x, data = data.frame(cbind(x, y)))
Residuals:
Min 1Q Median 3Q Max
-1.007 -0.6767 -0.2846 0.3867 5.435
Coefficients:
Value Std. Error t value Pr(>|t|)
(Intercept) 0.9075 0.0643 14.1220 0.0000
x 1.0140 0.0173 58.7814 0.0000
Residual standard error: 0.9113 on 498 degrees of freedom
Multiple R-Squared: 0.874
F-statistic: 3455 on 1 and 498 degrees of freedom, the p-value is 0
Correlation of Coefficients:
(Intercept)
x -0.7732
________________________________________________________________output
$er.coef:
(Intercept) x
0.0642591 0.01725025
_________________________________________________________
Now we do the same with gls
> now.gls()
Generalized least squares fit by REML
Model: y ~ x
Data: data.frame(cbind(x, y))
AIC BIC logLik
1449.718 1466.56 -720.859
Variance function:
Structure: Power of variance covariate
Formula: ~ fitted(.)
Parameter estimates:
power
0.1252992
Coefficients:
Value Std.Error t-value p-value
(Intercept) 0.897721 0.06998511 12.82731 <.0001
x 1.024476 0.02171861 47.17042 <.0001
Correlation:
(Intr)
x -0.764
Standardized residuals:
Min Q1 Med Q3 Max
-1.027471 -0.7095038 -0.3511753 0.3350945 5.44721
Residual standard error: 0.8791239
Degrees of freedom: 500 total; 498 residual
_________________________________________________________output
Error in [: No dim attribute for array subset: structure(....
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
Apparently, the structure of summary$coefficients is similar in both
cases
However, the command "summary$coefficients" or "coefficients(summary)"
returns 2X2 matrix in case of lm and 2X1 matrix in case of gls
Is it a bug, or I am missing something?
Is there a way around?
Thank you
Simon Rosenfeld
National Cancer Institute
Splus2000, WinNT 4.0
|