# 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?



