| ## some tests of inverse-gaussian GLMs based on a file supplied by |
| ## David Firth, Feb 2009. |
| |
| options(digits=5) |
| have_MASS <- requireNamespace('MASS', quietly = TRUE) |
| |
| ## Data from Whitmore, G A (1986), Inverse Gaussian Ratio Estimation. |
| ## Applied Statistics 35(1), 8-15. |
| ## |
| ## A "real, but disguised" set of data (Whitmore, 1986, p8). |
| ## |
| ## For each of 20 products, x is projected sales and y is actual sales. |
| |
| x <- c(5959, 3534, 2641, 1965, 1738, 1182, 667, 613, 610, 549, |
| 527, 353, 331, 290, 253, 193, 156, 133, 122, 114) |
| y <- c(5673, 3659, 2565, 2182, 1839, 1236, 918, 902, 756, 500, |
| 487, 463, 225, 257, 311, 212, 166, 123, 198, 99) |
| |
| ## Whitmore's model (2.4) is |
| |
| fit <- glm(y ~ x - 1, weights = x^2, |
| family = inverse.gaussian(link = "identity"), |
| epsilon = 1e-12) |
| fit |
| coef(summary(fit)) |
| ## Alternatively, use the explicit formula that's available for the MLE |
| ## in this example. It's just a ratio estimate: |
| (beta.exact <- sum(y)/sum(x)) |
| stopifnot(all.equal(beta.exact, as.vector(coef(fit)))) |
| ## and for a confidence interval via confint |
| if(have_MASS) { |
| ci <- confint(fit, 1, level = 0.95) |
| print(ci) |
| } |
| ## and via asymptotic normality |
| sterr <- coef(summary(fit))[, "Std. Error"] |
| coef(fit) + (1.96 * sterr * c(-1, 1)) |
| |
| ## David suggested the use of an inverse link |
| |
| fit2 <- glm(y ~ I(1/x) - 1, weights = x^2, |
| family = inverse.gaussian(link = "inverse"), |
| epsilon = 1e-12) |
| coef(summary(fit2)) |
| ## which gives the same CIs both ways |
| if(have_MASS) { |
| ci1 <- rev(1/(confint(fit2, 1, level = 0.95))) |
| print(ci1) |
| sterr <- (summary(fit2)$coefficients)[, "Std. Error"] |
| ci2 <- 1/(coef(fit2) - (1.96 * sterr * c(-1, 1))) |
| print(ci2) |
| stopifnot(all.equal(as.vector(ci), as.vector(ci1), tolerance = 1e-5), |
| all.equal(as.vector(ci), ci2, tolerance = 1e-3)) |
| } |
| ## because the log likelihood for 1/beta is exactly quadratic. |
| |
| ## The approximate intervals above differ slightly from the exact |
| ## confidence interval given in Whitmore (1986) -- as is to be |
| ## expected (they are based on asymptotic approximations, not the |
| ## exact pivot). |
| |
| |
| ## Now simulate from this model |
| if(requireNamespace("SuppDists")) { |
| print( ys <- simulate(fit, nsim = 3, seed = 1) ) |
| for(i in seq_len(3)) |
| print(coef(summary(update(fit, ys[, i] ~ .)))) |
| } |