From: John P. Burkett <burkett_at_uri.edu>

Date: Tue, 25 Mar 2008 15:16:57 -0400

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-9195Received on Tue 25 Mar 2008 - 20:10:31 GMT# 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.

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