[R] quasi-Poisson regression in stats glm vs aod quasipois

From: John P. Burkett <burkett_at_uri.edu>
Date: Tue, 25 Mar 2008 15:16:57 -0400


For quasi-Poisson regression, I have tried to compare the output of the quasipois function in the aod package to the output of the glm function with family=quasi(link="log", variance="mu^2") in the stats package. The standard errors are about 32.5 times larger in the former than in the latter. If anyone can explain the difference, I would be most grateful.

I am running R version 2.6.1 and aod version 1.1-24 under Gentoo Linux.

In case they are of diagnostic value, I am attaching my R code (compare.R), my data set (visitsbymonth1990_2006.csv), and my output file (compare.Rout).

Best regards,
John

-- 
John P. Burkett
Department of Environmental and Natural Resource Economics
and Department of Economics
University of Rhode Island
Kingston, RI 02881-0808
USA

phone (401) 874-9195


# This is a file of R commands to compare to functions designed to estimate # quasi-Poisson regressions. # Load data on recreational visits to the NPS system, 1990-2006. monthlyvisits <- read.csv("visitsbymonth1990_2006.csv",header=TRUE,sep=",") attach(monthlyvisits) # Put monthlyvisits in R search path mv.ts <- ts(monthlyvisits, start = c(1990, 1), frequency = 1) # Form time series # so as to associate a date with each observation Jan.ts <- ts(Jan, start = c(1990, 1), frequency = 1) year.ts <- 1990:2006 jydf <- data.frame(Jan.ts,year.ts) # start regression analysis using data for 1990-2006 Jan.quasiP <- glm(Jan.ts ~ year.ts, family=quasipoisson()) summary(Jan.quasiP) Jan.quasi1 <- glm(Jan.ts ~ year.ts, family=quasi(link="log", variance="mu")) summary(Jan.quasi1) # Note that Jan.quasiP and Jan.quasi1 are identical. Jan.quasi2 <- glm(Jan.ts ~ year.ts, family=quasi(link="log", variance="mu^2")) summary(Jan.quasi2) library(aod) Jan.qp.aod <- quasipois(formula = Jan.ts ~ year.ts, data = jydf) Jan.qp.aod # I expected Jan.qp.aod to be the same as Jan.quasi2; but the # standard errors are about 32.5 times larger in the former than the latter. # Why is that?

R version 2.6.1 (2007-11-26) Copyright (C) 2007 The R Foundation for Statistical Computing ISBN 3-900051-07-0 R is free software and comes with ABSOLUTELY NO WARRANTY. You are welcome to redistribute it under certain conditions. Type 'license()' or 'licence()' for distribution details. Natural language support but running in an English locale R is a collaborative project with many contributors. Type 'contributors()' for more information and 'citation()' on how to cite R or R packages in publications. Type 'demo()' for some demos, 'help()' for on-line help, or 'help.start()' for an HTML browser interface to help. Type 'q()' to quit R. [Previously saved workspace restored]
> # This is a file of R commands to compare to functions designed to estimate
> # quasi-Poisson regressions.
> # Load data on recreational visits to the NPS system, 1990-2006.
> monthlyvisits <- read.csv("visitsbymonth1990_2006.csv",header=TRUE,sep=",")
> attach(monthlyvisits) # Put monthlyvisits in R search path
> mv.ts <- ts(monthlyvisits, start = c(1990, 1), frequency = 1) # Form time series
> # so as to associate a date with each observation
> Jan.ts <- ts(Jan, start = c(1990, 1), frequency = 1)
> year.ts <- 1990:2006
> jydf <- data.frame(Jan.ts,year.ts)
> # start regression analysis using data for 1990-2006
> Jan.quasiP <- glm(Jan.ts ~ year.ts, family=quasipoisson())
> summary(Jan.quasiP)
Call: glm(formula = Jan.ts ~ year.ts, family = quasipoisson()) Deviance Residuals: Min 1Q Median 3Q Max -596.92 -47.10 29.76 86.15 279.24 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) -2.335775 6.420994 -0.364 0.7211 year.ts 0.009274 0.003213 2.886 0.0113 * --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 (Dispersion parameter for quasipoisson family taken to be 45416.02) Null deviance: 1079945 on 16 degrees of freedom Residual deviance: 701426 on 15 degrees of freedom AIC: NA Number of Fisher Scoring iterations: 3
> Jan.quasi1 <- glm(Jan.ts ~ year.ts, family=quasi(link="log", variance="mu"))
> summary(Jan.quasi1) # Note that Jan.quasiP and Jan.quasi1 are identical.
Call: glm(formula = Jan.ts ~ year.ts, family = quasi(link = "log", variance = "mu")) Deviance Residuals: Min 1Q Median 3Q Max -596.92 -47.10 29.76 86.15 279.24 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) -2.335775 6.420994 -0.364 0.7211 year.ts 0.009274 0.003213 2.886 0.0113 * --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 (Dispersion parameter for quasi family taken to be 45416.02) Null deviance: 1079945 on 16 degrees of freedom Residual deviance: 701426 on 15 degrees of freedom AIC: NA Number of Fisher Scoring iterations: 3
> Jan.quasi2 <- glm(Jan.ts ~ year.ts, family=quasi(link="log", variance="mu^2"))
> summary(Jan.quasi2)
Call: glm(formula = Jan.ts ~ year.ts, family = quasi(link = "log", variance = "mu^2")) Deviance Residuals: Min 1Q Median 3Q Max -0.189726 -0.015445 0.009383 0.026420 0.085374 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) -2.168754 6.490702 -0.334 0.7429 year.ts 0.009190 0.003249 2.829 0.0127 * --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 (Dispersion parameter for quasi family taken to be 0.004305760) Null deviance: 0.103446 on 16 degrees of freedom Residual deviance: 0.068676 on 15 degrees of freedom AIC: NA Number of Fisher Scoring iterations: 4
> library(aod)
Package aod, version 1.1-24
> Jan.qp.aod <- quasipois(formula = Jan.ts ~ year.ts, data = jydf)
> Jan.qp.aod # I expected Jan.qp.aod to be the same as Jan.quasi2; but the
Quasi-likelihood generalized linear model ----------------------------------------- quasipois(formula = Jan.ts ~ year.ts, data = jydf) Fixed-effect coefficients: Estimate Std. Error z value Pr(>|z|) (Intercept) -2.1688 210.7997 -0.0103 0.9918 year.ts 0.0092 0.1055 0.0871 0.9306 Overdispersion parameter: phi 4.5416 Pearson's chi-squared goodness-of-fit statistic = 0.0142
> # standard errors are about 32.5 times larger in the former than the latter.
> # Why is that?
>
>
>
>
> proc.time()
user system elapsed 1.176 0.022 1.186

______________________________________________ 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.

Received on Tue 25 Mar 2008 - 20:10:31 GMT

Archive maintained by Robert King, hosted by the discipline of statistics at the University of Newcastle, Australia.
Archive generated by hypermail 2.2.0, at Tue 25 Mar 2008 - 20:30:24 GMT.

Mailing list information is available at https://stat.ethz.ch/mailman/listinfo/r-help. Please read the posting guide before posting to the list.

list of date sections of archive