[R] bigglm() results different from glm()
Francisco J. Zagmutt
gerifalte28 at hotmail.com
Wed Mar 18 02:33:14 CET 2009
Dear Thomas and John,
Thanks for your prompt reply and for looking at your code to explain
these differences. I see you replied very late at night, so I am sorry
if my little problem kept you awake!
As you pointed out, the model indeed converges (as reported in
model$converged) once I specify a larger number of iterations.
A very minor comment: it seems that the reporting of the estimates in
summary.biglm() is not taking the parameters from options(digits). For
example, using the same data and models as before:
> require(biglm)
> options(digits=8)
> dat=data.frame(y =c(rpois(50000, 10),rpois(50000, 15)),
ttment=gl(2,50000))
> m1 <- glm(y~ttment, data=dat, family=poisson(link="log"))
> m1big <- bigglm(y~ttment , data=dat, family=poisson(link="log"),
maxit=20)
> summary(m1)
<Snipped>
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 2.3019509 0.0014147 1627.21 < 2.2e-16 ***
ttment2 0.4052002 0.0018264 221.86 < 2.2e-16 ***
> summary(m1big)
Coef (95% CI) SE p
(Intercept) 2.302 2.299 2.305 0.001 0
ttment2 0.405 0.402 0.409 0.002 0
To get more digits I can extract the point estimates using coef(m1big),
but after looking at str(m1big) the only way I could figure to extract
the p-values was:
> summary(m1big)$mat[,"p"]
(Intercept) ttment2
0 0
Is there a way I can get summary.biglm() to report more digits directly?
Thanks again,
Francisco
--
Francisco J. Zagmutt
Vose Consulting
2891 20th Street
Boulder, CO, 80304
USA
www.voseconsulting.com
Thomas Lumley wrote:
>
> Yes, the slow convergence is easier to get with the log link.
> Overshooting the correct coefficient vector has more dramatic effects on
> the fitted values and weights, and in this example the starting value of
> (0,0) is a long way from the truth. The same sort of thing happens in
> the Cox model, where there are real data sets that will cause numeric
> overflow in a careless implementation.
>
> It might be worth trying to guess better starting values: saving an
> iteration or two would be useful with large data sets.
>
> -thomas
>
>
> On Tue, 17 Mar 2009, John Fox wrote:
>
>> Dear Francisco,
>>
>> I was able to duplicate the problem that you reported, and in addition
>> discovered that the problem seems to be peculiar to the poisson family.
>> lm(y~ttment, data=dat) and biglm(y~ttment, data=dat) produce identical
>> results, as do glm(y~ttment, data=dat) and bigglm(y~ttment, data=dat).
>> Another example, with the binomial family:
>>
>>> set.seed(12345)
>>> dat=data.frame(y =c(rpois(50000, 10),rpois(50000, 15)),
>> ttment=gl(2,50000))
>>> dat$y0 <- ifelse(dat$y > 12, 1, 0)
>>> m1b <- glm(y0~ttment, data=dat, family=binomial)
>>> m1bigb <- bigglm(y0~ttment , data=dat, family=binomial())
>>> coef(m1b)
>> (Intercept) ttment2
>> -1.33508 2.34263
>>> coef(m1bigb)
>> (Intercept) ttment2
>> -1.33508 2.34263
>>> deviance(m1b)
>> [1] 109244
>>> deviance(m1bigb)
>> [1] 109244
>>
>> I'm copying this message to Thomas Lumley, who's the author and
>> maintainer
>> of the biglm package.
>>
>> Regards,
>> John
>>
>>
>>> -----Original Message-----
>>> From: r-help-bounces at r-project.org [mailto:r-help-bounces at r-project.org]
>> On
>>> Behalf Of Francisco J. Zagmutt
>>> Sent: March-16-09 10:26 PM
>>> To: r-help at stat.math.ethz.ch
>>> Subject: [R] bigglm() results different from glm()
>>>
>>> Dear all,
>>>
>>> I am using the bigglm package to fit a few GLM's to a large dataset (3
>>> million rows, 6 columns). While trying to fit a Poisson GLM I noticed
>>> that the coefficient estimates were very different from what I obtained
>>> when estimating the model on a smaller dataset using glm(), I wrote a
>>> very basic toy example to compare the results of bigglm() against a
>>> glm() call. Consider the following code:
>>>
>>>
>>> > require(biglm)
>>> > options(digits=6, scipen=3, contrasts = c("contr.treatment",
>>> "contr.poly"))
>>> > dat=data.frame(y =c(rpois(50000, 10),rpois(50000, 15)),
>>> ttment=gl(2,50000))
>>> > m1 <- glm(y~ttment, data=dat, family=poisson(link="log"))
>>> > m1big <- bigglm(y~ttment , data=dat, family=poisson(link="log"))
>>> > summary(m1)
>>>
>>> <snipped output for this email>
>>> Coefficients:
>>> Estimate Std. Error z value Pr(>|z|)
>>> (Intercept) 2.30305 0.00141 1629 <2e-16 ***
>>> ttment2 0.40429 0.00183 221 <2e-16 ***
>>>
>>> Null deviance: 151889 on 99999 degrees of freedom
>>> Residual deviance: 101848 on 99998 degrees of freedom
>>> AIC: 533152
>>>
>>> > summary(m1big)
>>> Large data regression model: bigglm(y ~ ttment, data = dat, family =
>>> poisson(link = "log"))
>>> Sample size = 100000
>>> Coef (95% CI) SE p
>>> (Intercept) 2.651 2.650 2.653 0.001 0
>>> ttment2 4.346 4.344 4.348 0.001 0
>>>
>>> > m1big$deviance
>>> [1] 287158986
>>>
>>>
>>> Notice that the coefficients and deviance are quite different in the
>>> model estimated using bigglm(). If I change the chunk to
>>> seq(1000,10000,1000) the estimates remain the same.
>>>
>>> Can someone help me understand what is causing these differences?
>>>
>>> Here is my version info:
>>>
>>> > version
>>> _
>>> platform i386-pc-mingw32
>>> arch i386
>>> os mingw32
>>> system i386, mingw32
>>> status
>>> major 2
>>> minor 8.1
>>> year 2008
>>> month 12
>>> day 22
>>> svn rev 47281
>>> language R
>>> version.string R version 2.8.1 (2008-12-22)
>>>
>>>
>>> Many thanks in advance for your help,
>>>
>>> Francisco
>>>
>>> --
>>> Francisco J. Zagmutt
>>> Vose Consulting
>>> 2891 20th Street
>>> Boulder, CO, 80304
>>> USA
>>> www.voseconsulting.com
>>>
>>> ______________________________________________
>>> R-help at r-project.org 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.
>>
>> ______________________________________________
>> R-help at r-project.org 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.
>>
>
> Thomas Lumley Assoc. Professor, Biostatistics
> tlumley at u.washington.edu University of Washington, Seattle
>
> ______________________________________________
> R-help at r-project.org 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.
>
More information about the R-help
mailing list