R version 3.1.1 (2014-07-10) -- "Sock it to Me"
Copyright (C) 2014 The R Foundation for Statistical Computing
Platform: x86_64-unknown-linux-gnu (64-bit)

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.

> 
>  pdf(file = "glm2-plots.pdf")
> 
>  set.seed(42)
> 
>  x <- 1:51
>  z <- (x - mean(x)) / sd(x)
>  z <- z / max(z) * 2.0
>  range(z)
[1] -2  2
>  fred <- function(x) ifelse(x > 0, 1 / (1 + exp(- x)), exp(x) / (1 + exp(x)))
>  moo <- fred(z)
> 
>  plot(x, moo, ylim = c(0, 1))
>  title(main = "Simulation Truth Mean Values")
> 
>  y <- as.numeric(runif(length(x)) < moo)
>  plot(x, y)
>  title(main = "Simulated Data")
> 
>  # all links
>  links <- c("logit", "probit", "cloglog", "cauchit", "log", "identity")
>  cc <- c("black", "red", "darkgreen", "blue", "magenta", "cyan")
>  plot(x, y)
>  for (i in seq(along = links)) {
+      if (links[i] != "identity") {
+          out <- glm(y ~ x, family = binomial(link = links[i]),
+              start = c(-1, 0))
+      } else {
+          out <- glm(y ~ x, family = quasi(variance = "mu(1-mu)"),
+              start = c(0.5, 0))
+      }
+      curve(predict(out, type = "response", newdata = data.frame(x = x)),
+          add = TRUE, col = cc[i])
+  }
>  legend(0, 0.8, legend = links, col = cc, lwd = 2)
>  title(main = "All Link Functions R Function GLM Knows About")
> 
>  ########## Poisson Regression ##########
> 
>  x <- 1:51
>  mu <- exp(x / max(x)) * 10
>  y <- rpois(length(mu), mu)
> 
>  out <- glm(y ~ x, family = poisson)
>  summary(out)

Call:
glm(formula = y ~ x, family = poisson)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-3.4425  -0.9087   0.0390   0.6036   2.5460  

Coefficients:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept) 2.321295   0.077668  29.888   <2e-16 ***
x           0.019325   0.002327   8.304   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for poisson family taken to be 1)

    Null deviance: 130.123  on 50  degrees of freedom
Residual deviance:  59.498  on 49  degrees of freedom
AIC: 299.51

Number of Fisher Scoring iterations: 4

>  plot(x, y)
>  curve(predict(out, type = "response", newdata = data.frame(x = x)),
+      add = TRUE)
> 
>  ########## Overdispersed Poisson Regression ##########
> 
>  y <- rnbinom(length(mu), size = pi, mu = mu)
>  out <- glm(y ~ x, family = quasipoisson)
>  summary(out)

Call:
glm(formula = y ~ x, family = quasipoisson)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-3.8828  -1.8111  -0.3239   1.2142   5.4433  

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 2.143135   0.201617  10.630 2.55e-14 ***
x           0.022577   0.005943   3.799 0.000402 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for quasipoisson family taken to be 5.92317)

    Null deviance: 367.05  on 50  degrees of freedom
Residual deviance: 278.74  on 49  degrees of freedom
AIC: NA

Number of Fisher Scoring iterations: 5

>  plot(x, y)
>  curve(predict(out, type = "response", newdata = data.frame(x = x)),
+      add = TRUE)
>  out <- glm(y ~ x, family = poisson)
>  curve(predict(out, type = "response", newdata = data.frame(x = x)),
+      add = TRUE, col = "red")
> 
> 
> 
> proc.time()
   user  system elapsed 
  0.278   0.027   0.396