[R] extend summary.lm for hccm?

Achim Voß achim.voss at uni-muenster.de
Tue Dec 30 14:15:37 CET 2008


Hi!


I am trying to estimate Engel curves using a big sample (>42,000) using lm and
taking heteroskedasticity into account by using the summaryHCCM posted here by
John Fox (Mon Dec 25 16:01:59 CET 2006).

Having used the SIC (with MASS stepAIC) to determine how many powers to use I
estimate the model:

> # =========================================
> summary.lm(fit.lm.5)

Call:
lm(formula = energyshare ~ 1 + I(x.log) + I(x.log^2) + I(x.log^3) +
    I(x.log^4) + I(x.log^5), data = ev)

Residuals:
      Min        1Q    Median        3Q       Max
-0.098819 -0.023682 -0.007043  0.013924  0.486615

Coefficients:
              Estimate Std. Error t value Pr(>|t|)
(Intercept) -5.452e+01  1.234e+01  -4.418 9.97e-06 ***
I(x.log)     3.177e+01  6.966e+00   4.560 5.13e-06 ***
I(x.log^2)  -7.330e+00  1.567e+00  -4.677 2.93e-06 ***
I(x.log^3)   8.395e-01  1.757e-01   4.778 1.78e-06 ***
I(x.log^4)  -4.775e-02  9.814e-03  -4.865 1.15e-06 ***
I(x.log^5)   1.079e-03  2.185e-04   4.939 7.90e-07 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.03748 on 42738 degrees of freedom
Multiple R-squared: 0.09236,    Adjusted R-squared: 0.09226
F-statistic: 869.8 on 5 and 42738 DF,  p-value: < 2.2e-16
> # =========================================




Now I use summaryHCCM(fit.lm.5):

> # =========================================
> summaryHCCM(fit.lm.5)
Fehler in solve.default(L %*% V %*% t(L)) :
  System ist für den Rechner singulär: reziproke Konditionszahl = 6.98689e-19
> # =========================================

("Error in solve.default(L %*% V %*% t(L)) :
  System is singulary for the computer: reciprocal number of conditions =
  6.98689e-19")



This does not happen if I omit I(x.log^5). I do not know what it means and I'd
be grateful if anyone could help. And I'd like to add a (more or less) related
question: Can I still use AIC, SIC etc. if I know there's a heteroskedasticity
problem?


Thanks in advance,
Achim



More information about the R-help mailing list