
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 
