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