``````## Paul E. Johnson
## 2014-09-18

## Here's a puzzle I presented to our group this week.

## Generate data like so

set.seed(234234)
dat <- data.frame(x1 = rnorm(100, m = 50, sd = 40), x2 = gl(4, 25,
labels = c("l", "m", "h", "s")), x3 = rgamma(100, 1, 2))
dat\$y <- rpois(100, lambda = exp(0.005 + 0.014 * dat\$x1 - 0.1 * dat\$x3))
m1 <- glm(y ~ x1 + x2 + x3, data = dat, family = "poisson")
m2 <- glm(y ~ x1, data = dat, family = "poisson")

## Exercises

## 1. Make a plot demonstrating the relationship between x1 and y,
## including the predicted values from both m1 and m2.

## 2. Suppose you want to plot the predicted value of y, for each
## categorical value of x2.  What would you do?

## 2A. What is the difference between

predict(m1)
``````
``````##        1        2        3        4        5        6        7        8
##  0.78761  0.57589  0.93523  0.57266  0.76530  1.28713  0.93793  0.98241
##        9       10       11       12       13       14       15       16
##  1.31750  0.68068  1.43293  0.73529  1.10346  0.72862  0.56668  0.94110
##       17       18       19       20       21       22       23       24
##  1.94896  0.04967  0.96010  0.19668  1.41698  1.35944  0.04705  0.42981
##       25       26       27       28       29       30       31       32
## -0.22183  1.34179  0.33170 -0.01971  0.39195  0.99532  0.54241  0.48098
##       33       34       35       36       37       38       39       40
##  0.73244  0.50376  0.30670  0.64904  1.21870 -0.16999  1.12084  1.01777
##       41       42       43       44       45       46       47       48
##  0.90110  1.19787  0.81861  0.21943  1.06838  0.07575  0.76713  2.13945
##       49       50       51       52       53       54       55       56
##  0.15473  0.74920 -0.52417 -0.31611  0.68070 -0.02657  1.10092  0.02224
##       57       58       59       60       61       62       63       64
##  1.01728  0.32843 -0.33360  1.40265  0.01630  0.40897  0.22931  0.51323
##       65       66       67       68       69       70       71       72
##  0.32363  0.78681  0.29791  0.85760  0.64160 -0.13928 -0.12438  0.25258
##       73       74       75       76       77       78       79       80
##  0.59972  0.85563 -0.08751  1.02280  0.46102  0.07333  0.34797 -0.01053
##       81       82       83       84       85       86       87       88
## -0.01263  1.20155 -0.50061  0.57326  1.03216  0.26116 -0.45316  0.63864
##       89       90       91       92       93       94       95       96
##  0.56602  0.72965  0.38354  0.75460  0.42533  0.39742  0.97632  0.52491
##       97       98       99      100
##  0.50616  0.78810  0.83690  0.62237
``````
``````## and

predict(m1, type = "response")
``````
``````##      1      2      3      4      5      6      7      8      9     10
## 2.1981 1.7787 2.5478 1.7730 2.1496 3.6224 2.5547 2.6709 3.7341 1.9752
##     11     12     13     14     15     16     17     18     19     20
## 4.1909 2.0861 3.0146 2.0722 1.7624 2.5628 7.0214 1.0509 2.6120 1.2174
##     21     22     23     24     25     26     27     28     29     30
## 4.1246 3.8940 1.0482 1.5370 0.8010 3.8259 1.3933 0.9805 1.4799 2.7056
##     31     32     33     34     35     36     37     38     39     40
## 1.7202 1.6177 2.0802 1.6549 1.3589 1.9137 3.3828 0.8437 3.0674 2.7670
##     41     42     43     44     45     46     47     48     49     50
## 2.4623 3.3131 2.2673 1.2454 2.9107 1.0787 2.1536 8.4948 1.1673 2.1153
##     51     52     53     54     55     56     57     58     59     60
## 0.5920 0.7290 1.9753 0.9738 3.0069 1.0225 2.7657 1.3888 0.7163 4.0660
##     61     62     63     64     65     66     67     68     69     70
## 1.0164 1.5053 1.2577 1.6707 1.3821 2.1964 1.3470 2.3575 1.8995 0.8700
##     71     72     73     74     75     76     77     78     79     80
## 0.8830 1.2873 1.8216 2.3529 0.9162 2.7810 1.5857 1.0761 1.4162 0.9895
##     81     82     83     84     85     86     87     88     89     90
## 0.9874 3.3253 0.6062 1.7740 2.8071 1.2984 0.6356 1.8939 1.7612 2.0744
##     91     92     93     94     95     96     97     98     99    100
## 1.4675 2.1268 1.5301 1.4880 2.6547 1.6903 1.6589 2.1992 2.3092 1.8633
``````
``````## 2B. What is the difference between

predict(m1, type = "response")
``````
``````##      1      2      3      4      5      6      7      8      9     10
## 2.1981 1.7787 2.5478 1.7730 2.1496 3.6224 2.5547 2.6709 3.7341 1.9752
##     11     12     13     14     15     16     17     18     19     20
## 4.1909 2.0861 3.0146 2.0722 1.7624 2.5628 7.0214 1.0509 2.6120 1.2174
##     21     22     23     24     25     26     27     28     29     30
## 4.1246 3.8940 1.0482 1.5370 0.8010 3.8259 1.3933 0.9805 1.4799 2.7056
##     31     32     33     34     35     36     37     38     39     40
## 1.7202 1.6177 2.0802 1.6549 1.3589 1.9137 3.3828 0.8437 3.0674 2.7670
##     41     42     43     44     45     46     47     48     49     50
## 2.4623 3.3131 2.2673 1.2454 2.9107 1.0787 2.1536 8.4948 1.1673 2.1153
##     51     52     53     54     55     56     57     58     59     60
## 0.5920 0.7290 1.9753 0.9738 3.0069 1.0225 2.7657 1.3888 0.7163 4.0660
##     61     62     63     64     65     66     67     68     69     70
## 1.0164 1.5053 1.2577 1.6707 1.3821 2.1964 1.3470 2.3575 1.8995 0.8700
##     71     72     73     74     75     76     77     78     79     80
## 0.8830 1.2873 1.8216 2.3529 0.9162 2.7810 1.5857 1.0761 1.4162 0.9895
##     81     82     83     84     85     86     87     88     89     90
## 0.9874 3.3253 0.6062 1.7740 2.8071 1.2984 0.6356 1.8939 1.7612 2.0744
##     91     92     93     94     95     96     97     98     99    100
## 1.4675 2.1268 1.5301 1.4880 2.6547 1.6903 1.6589 2.1992 2.3092 1.8633
``````
``````## and

predict(m2, type = "response")
``````
``````##      1      2      3      4      5      6      7      8      9     10
## 1.7928 1.2972 2.1088 1.4514 1.7367 2.7586 2.0980 2.0982 3.1104 1.3746
##     11     12     13     14     15     16     17     18     19     20
## 3.4362 1.6375 2.3996 1.5612 1.3870 2.0594 5.8368 0.8409 2.1651 0.9724
##     21     22     23     24     25     26     27     28     29     30
## 3.5021 3.2569 0.8194 1.1892 0.6045 3.8685 1.3651 0.9339 1.1128 2.6954
##     31     32     33     34     35     36     37     38     39     40
## 1.5926 1.6230 2.0055 1.6702 1.3396 1.9145 3.5025 0.8329 3.1139 2.6333
##     41     42     43     44     45     46     47     48     49     50
## 2.4356 3.4123 2.1988 1.2007 2.8922 1.0642 2.0026 7.1653 1.1177 1.9934
##     51     52     53     54     55     56     57     58     59     60
## 0.7139 0.8753 2.5741 1.2623 3.5495 1.1535 3.6667 1.8268 0.8701 5.5313
##     61     62     63     64     65     66     67     68     69     70
## 1.3202 1.8811 1.6449 2.0411 1.6702 2.6356 1.6910 3.1345 2.4700 1.1113
##     71     72     73     74     75     76     77     78     79     80
## 1.1154 1.5746 2.3968 2.8858 1.1378 2.8355 1.5522 1.1793 1.5106 1.0968
##     81     82     83     84     85     86     87     88     89     90
## 1.0317 3.8087 0.6507 2.0004 2.9785 1.3665 0.6506 1.9637 1.9971 2.3624
##     91     92     93     94     95     96     97     98     99    100
## 1.6549 2.4183 1.6263 1.6550 2.9992 1.8874 1.8589 2.4455 2.6275 1.9267
``````
``````## 3. What's wrong with this approach?

pred1 <- predict(m1)
plot(pred1 ~ dat\$x1, type = "l")
`````` ``````## 4. What's wrong with this approach?

dat <- dat[order(dat\$x1), ]
m1 <- glm(y ~ x1 + x2 + x3, data = dat, family = "poisson")
pred2 <- predict(m1)
plot(y ~ x1, data = dat)
lines(pred2 ~ x1, data = dat)

## 5. What's wrong with this approach?

library(rockchalk)
`````` ``````nd <- newdata(m1, predVals = c("x1"), n = 20)
nd\$pred <- predict(m1, nd)
plot(y ~ x1, data = dat)
lines(pred ~ x1, data = nd)
`````` ``````## 6. Can you make a pleasant looking regression table for m1 and m2?
## If not, try R packages that can like memisc, texreg, stargazer,
## apsrtable, rockchalk, or ...

## 7. More advanced: what's the problem with getting confidence
## intervals on that predicted value line?

## Can you explain why I feel like a cheater when I wrote this function?

plotSlopes(m1, plotx = x1, interval = "conf")
``````
``````## rockchalk:::predCI: model's predict method does not return an interval. We will improvize with a Wald type approximation to the confidence interval
`````` ``````## If you did not ever try to work with predicted value plots, I
## eagerly suggest you try it.

## 1. Make a plot demonstrating the relationship between x1 and y,
## including the predicted values from both m1 and m2.

## m2 is easy. Do that first
m2nd <- newdata(m2, predVals = list(x1 = "quantile"), n = 30)
m2nd\$fit <- predict(m2, newdata = m2nd, type = "response")

plot(y ~ x1, data = dat, col = "gray70")
lines(fit ~ x1, data = m2nd)
`````` ``````## Perhaps it is more clear to look at it like this
plot(y ~ x1, data = dat, col = "gray70")
lines(fit ~ x1, data = m2nd, type = "b")
title("30 quantile-spaced values of x1")
`````` ``````## Obvious problem: What values of x1 should be used? 30 quantiles?

m2nd2 <- newdata(m2, predVals = list(x1 = "seq"), n = 30)
m2nd2\$fit <- predict(m2, newdata = m2nd2, type = "response")
plot(y ~ x1, data = dat, col = "gray70")
lines(fit ~ x1, data = m2nd2, type = "b")
title("Evenly spaced values of x1")

## Obvious problem. What values are to be supposed for x2 and x3?
m1nd <- newdata(m1, predVals = list(x1 = "quantile", x2 = levels(dat\$x2)), n = 20)
m1nd
``````
``````##         x1 x2     x3
## 1  -28.737  l 0.4613
## 2   -7.123  l 0.4613
## 3    3.132  l 0.4613
## 4   12.508  l 0.4613
## 5   16.146  l 0.4613
## 6   23.781  l 0.4613
## 7   27.231  l 0.4613
## 8   35.789  l 0.4613
## 9   38.822  l 0.4613
## 10  40.438  l 0.4613
## 11  47.617  l 0.4613
## 12  51.380  l 0.4613
## 13  52.823  l 0.4613
## 14  56.367  l 0.4613
## 15  64.634  l 0.4613
## 16  69.572  l 0.4613
## 17  74.276  l 0.4613
## 18  79.927  l 0.4613
## 19  88.331  l 0.4613
## 20  93.273  l 0.4613
## 21 138.443  l 0.4613
## 22 -28.737  m 0.4613
## 23  -7.123  m 0.4613
## 24   3.132  m 0.4613
## 25  12.508  m 0.4613
## 26  16.146  m 0.4613
## 27  23.781  m 0.4613
## 28  27.231  m 0.4613
## 29  35.789  m 0.4613
## 30  38.822  m 0.4613
## 31  40.438  m 0.4613
## 32  47.617  m 0.4613
## 33  51.380  m 0.4613
## 34  52.823  m 0.4613
## 35  56.367  m 0.4613
## 36  64.634  m 0.4613
## 37  69.572  m 0.4613
## 38  74.276  m 0.4613
## 39  79.927  m 0.4613
## 40  88.331  m 0.4613
## 41  93.273  m 0.4613
## 42 138.443  m 0.4613
## 43 -28.737  h 0.4613
## 44  -7.123  h 0.4613
## 45   3.132  h 0.4613
## 46  12.508  h 0.4613
## 47  16.146  h 0.4613
## 48  23.781  h 0.4613
## 49  27.231  h 0.4613
## 50  35.789  h 0.4613
## 51  38.822  h 0.4613
## 52  40.438  h 0.4613
## 53  47.617  h 0.4613
## 54  51.380  h 0.4613
## 55  52.823  h 0.4613
## 56  56.367  h 0.4613
## 57  64.634  h 0.4613
## 58  69.572  h 0.4613
## 59  74.276  h 0.4613
## 60  79.927  h 0.4613
## 61  88.331  h 0.4613
## 62  93.273  h 0.4613
## 63 138.443  h 0.4613
## 64 -28.737  s 0.4613
## 65  -7.123  s 0.4613
## 66   3.132  s 0.4613
## 67  12.508  s 0.4613
## 68  16.146  s 0.4613
## 69  23.781  s 0.4613
## 70  27.231  s 0.4613
## 71  35.789  s 0.4613
## 72  38.822  s 0.4613
## 73  40.438  s 0.4613
## 74  47.617  s 0.4613
## 75  51.380  s 0.4613
## 76  52.823  s 0.4613
## 77  56.367  s 0.4613
## 78  64.634  s 0.4613
## 79  69.572  s 0.4613
## 80  74.276  s 0.4613
## 81  79.927  s 0.4613
## 82  88.331  s 0.4613
## 83  93.273  s 0.4613
## 84 138.443  s 0.4613
``````
``````## See the problem?:
m1nd\$fit <- predict(m1, newdata = m1nd, type = "response")
text(m1nd\$x1, m1nd\$fit, labels = as.character(m1nd\$x2))
`````` ``````plot(y ~ x1, data = dat, col = "gray70")
mycol <- gray.colors(4, start = 0.3, end = 0.6)
x2levels <- levels(m1nd\$x2)
for(i in seq_along(x2levels)){
lines(fit ~ x1, data = m1nd[m1nd\$x2 %in% x2levels[i], ], lty = i + 1, col = mycol[i], lwd = 2)
}
lines(fit ~ x1, data = m2nd)
legend("topleft", legend = c("Model 2", paste("Model 1: x2 =", x2levels)), col = c(1, mycol), lty = 1:5, lwd = 2)
`````` ``````## The rockchalk package has a short-cut function for most of this work, but it won't integrate
## m1 and m2 for us.

## ## 4. What's wrong with this approach?

## dat <- dat[order(dat\$x1), ]
## m1 <- glm(y ~ x1 + x2 + x3, data = dat, family = "poisson")
## pred2 <- predict(m1)
## plot(y ~ x1, data = dat)
## lines(pred2 ~ x1, data = dat)

## Predicted values are calculated for each row in the data. The
## plot looks/is wrong because we have not controlled for x2 and x3
## in the calculation of the predicted values for plotting.

## ## 5. What's wrong with this approach?

## library(rockchalk)
## nd <- newdata(m1, predVals = c("x1"), n = 20)
## nd\$pred <- predict(m1, nd)
## plot(y ~ x1, data = dat)
## lines(pred ~ x1, data = nd)

## That's almost right. Inspect nd, see it has adjusted x2 and x3.
## However, the predict statement needs type = "response"

library(rockchalk)
nd <- newdata(m1, predVals = c("x1"), n = 20)
nd\$fit <- predict(m1, nd, type = "response")
plot(y ~ x1, data = dat)
lines(fit ~ x1, data = nd)
`````` ``````## do you want various predictive lines?

plotCurves(m1, plotx = "x1", modx = "x2")
`````` ``````##
##
## Now, I want to take you on 2 detours.

## Detour 1. Plot predictions by quantile.

## Lets just ask for 5 values

nd2 <- newdata(m1, predVals = list(x1 = "quantile", x2 = "l"), n = 5)
nd2\$fit <- predict(m1, nd2, type = "response")

par(mfcol = c(2,1))

plot(y ~ x1, data = dat, col = "gray70")
lines(fit ~ x1, data = nd2, type = "b")

nd2\$q1 <- c("min", "25%", "median", "75%", "max")
nd2\$q2 <- 1:5

plot(fit ~ q2, data = nd2, type = "b", xlab = "x1 quantiles", axes = FALSE)
axis(2)
axis(1, at = nd2\$q2, labels = nd2\$q1)
box(bty = "L")
`````` ``````par(par.orig)

## If x1's distribution is less symmetric,the problem is exaggerated

set.seed(234234)
dat2 <- data.frame(x1 = 1.5 *  rgamma(100, shape = 2, scale = 18.3),
x2 = gl(4, 25, labels = c("l", "m", "h", "s")),
x3 = rgamma(100, 1, 2))
dat2\$y <- rpois(100, lambda = exp(0.005 + 0.014 * dat2\$x1 - 0.1 * dat2\$x3))

summary(dat)
``````
``````##        x1        x2           x3               y
##  Min.   :-28.7   l:25   Min.   :0.0046   Min.   :0.00
##  1st Qu.: 23.8   m:25   1st Qu.:0.1284   1st Qu.:1.00
##  Median : 47.6   h:25   Median :0.3293   Median :2.00
##  Mean   : 46.0   s:25   Mean   :0.4613   Mean   :2.06
##  3rd Qu.: 69.6          3rd Qu.:0.6377   3rd Qu.:3.00
##  Max.   :138.4          Max.   :2.8342   Max.   :9.00
``````
``````summary(dat2)
``````
``````##        x1         x2           x3              y
##  Min.   :  2.89   l:25   Min.   :0.002   Min.   : 0.00
##  1st Qu.: 27.75   m:25   1st Qu.:0.149   1st Qu.: 1.00
##  Median : 47.99   h:25   Median :0.354   Median : 2.00
##  Mean   : 54.72   s:25   Mean   :0.461   Mean   : 2.71
##  3rd Qu.: 69.98          3rd Qu.:0.703   3rd Qu.: 3.00
##  Max.   :174.79          Max.   :2.306   Max.   :16.00
``````
``````hist(dat2\$x1)
`````` ``````m1.2 <- glm(y ~ x1 + x2 + x3, data = dat2, family = "poisson")

nd3 <- newdata(m1.2, predVals = list(x1 = "quantile", x2 = "l"), n = 5)
nd3\$fit <- predict(m1.2, nd3, type = "response")

par(mfcol = c(2,1))

plot(y ~ x1, data = dat2, col = "gray70")
lines(fit ~ x1, data = nd3, type = "b")

nd3\$q1 <- c("min", "25%", "median", "75%", "max")
nd3\$q2 <- seq(min(nd3\$x1), max(nd3\$x1), length.out = 5)

plot(fit ~ q2, data = nd3, type = "b", xlab = "x1 quantiles", axes = FALSE, col = "gray70")
axis(2)
axis(1, at = nd3\$q2, labels = nd3\$q1)
box(bty = "L")
lines(fit ~ q2, data = nd3, type = "b")
`````` ``````par(par.orig)

## What are the dangers here?

## Why might one want to do this?

## Detour 2. What if the data has only 2 or 3 unique values of a predictor?

set.seed(234234)
dat <- data.frame(x1 = rnorm(100, m = 50, sd = 40),
x2 = gl(4, 25, labels = c("l", "m", "h", "s")),
x3 = rgamma(100, 1, 2))
dat\$x1f <- cut(dat\$x1, c(-100, 50, 80, 200), labels = c("low", "med", "high"))
dat\$x1n <- as.numeric(dat\$x1f)
dat\$y <- rpois(100, lambda = exp(0.005 + 0.014 * dat\$x1 - 0.1 * dat\$x3))

m1 <- glm(y ~ x1n + x2 + x3, data = dat, family = "poisson")
summary(m1)
``````
``````##
## Call:
## glm(formula = y ~ x1n + x2 + x3, family = "poisson", data = dat)
##
## Deviance Residuals:
##     Min       1Q   Median       3Q      Max
## -2.7697  -0.9699  -0.0865   0.5910   3.0815
##
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)
## (Intercept)  -0.2765     0.2315   -1.19    0.232
## x1n           0.5952     0.0888    6.70    2e-11 ***
## x2m          -0.1629     0.1841   -0.88    0.376
## x2h          -0.4082     0.2019   -2.02    0.043 *
## x2s          -0.1977     0.1985   -1.00    0.319
## x3            0.2404     0.1273    1.89    0.059 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
##     Null deviance: 167.37  on 99  degrees of freedom
## Residual deviance: 114.33  on 94  degrees of freedom
## AIC: 340.7
##
## Number of Fisher Scoring iterations: 5
``````
``````nd4 <- newdata(m1, predVals = list(x1n = "quantile", x2 = "l"), n = 5)
nd4\$fit <- predict(m1, newdata = nd4, type = "response")
nd4
``````
``````##   x1n x2     x3   fit
## 1   1  l 0.4613 1.537
## 2   2  l 0.4613 2.787
## 3   3  l 0.4613 5.053
``````
``````plot(fit ~ x1n, data = nd4, type = "b")
`````` ``````## Question: Should newdata refuse to give back more possible values of
## x1n than there actually are?

## I'll fake it by manufacturing data points at 1.5 and 2.5
x1quant <- quantile(dat\$x1n)
nd5 <- expand.grid(x1n = seq(1, 3, length.out = 5), x2 = "l", x3 = mean(dat\$x3))
nd5\$fit <- predict(m1, newdata = nd5, type = "response")
nd5
``````
``````##   x1n x2     x3   fit
## 1 1.0  l 0.4613 1.537
## 2 1.5  l 0.4613 2.069
## 3 2.0  l 0.4613 2.787
## 4 2.5  l 0.4613 3.752
## 5 3.0  l 0.4613 5.053
``````
``````plot(fit ~ x1n, data = nd5, type = "b", xlab = "phony (?) x1 quantiles", axes = FALSE)
axis(2)
axis(1, at = nd5\$x1n, labels = nd5\$q1)
box(bty = "L")
`````` ``````## ## 7. More advanced: what's the problem with getting confidence
## ## intervals on that predicted value line?

## ## Can you explain why I feel like a cheater when I wrote this function?

## plotSlopes(m1, plotx = x1, interval = "conf")

## ## If you did not ever try to work with predicted value plots, I
## ## eagerly suggest you try it.

## Those calculations rely on a very simplistic "Wald approximation"
## to calculate confidence intervals on predictions. These are also
## supported in SAS, but the literature is generally calling for more
## elaborate methods none of which are easily available.

## Withers, Christopher S. and Saralees Nadarah. 2012. Improved
## confidence regions based on Edgeworth expansions. Computational
## Statistics and Data Analysis 56: 4366-4380.

## "The number of papers proposing confidence regions is too many to
## cite. Confidence regions can be constructed using many different
## methods. Some of these involve: the exact distribution theory,
## asymptotic chi-square distributions, asymptotic normal distributions,
## record values, bootstrapping, empirical likelihood method, likelihood
## ratio statistics, Bayesian test statistics, sign based test
## statistics, rank based test statistics, differential geometric
## frameworks, jackknife, least square statistics, sequential methods,
## Monte Carlo likelihood methods, adaptive estimation methods, the M
## estimation theory, empirical Bayes, calibration methods, score
## statistics, quadratic approximations, subsampling methods, Gibbs
## sampling methods, smoothing methods (for example, polynomial splines,
## kernel type smoothing, local adaptive smoothing and Granulometric
## smoothing), Kolmogorov–Smirnov statistics, highest posterior density
## methods, R estimation methods, stochastic inequalities (for example,
## Hajek–Rényi type inequalities), multiple comparison methods, maximum
## modulus methods, minimum area confidence methods, U statistics, Rényi
## statistics and statistics based on other entropy measures, generalized
## p values, James–Stein type statistics, power divergence statistics,
## importance sampling methods, Uusipaikka’s method, genetic algorithms
## and uniform minimum variance unbiased statistics."

## King, Gary, Michael Tomz, and Jason Wittenberg. 2000. Making the Most
## of Statistical Analyses: Improving Interpretation and
## Presentation. American Journal of Political Science 44: 341–355

## Sun, J., Loader, C., & William P., M. (2000). Confidence bands in
## generalized linear models. The Annals of Statistics, 28(2),
## 429–460. doi:10.1214/aos/1016218225

## Xu, Jun and J. Scott Long (2005) Confidence intervals for predicted outcomes in
## regression models for categorical outcomes. Stata Journal 5(4): 537-559.

## 6. Can you make a pleasant looking regression table for m1 and m2?
## If not, try R packages that can like memisc, texreg, stargazer,
## apsrtable, rockchalk, or ...

varLabs <- c(x1 = "Flowers", x2l = "x2: low", x2m = "x2: medium",
x2h = "x2: high", x2s = "x2: super", x3 = "Nearly Forgot x3's super long label")
outreg(list("Model 1" = m1, "Model 2" = m2),  varLabels = varLabs, type = "html")
``````
``````##
##  We are launching a browser to view that html file, which we have temporarily
##  saved in  /tmp/RtmpFUk9T8/file10705954650b.html
##  You can copy that temp file, or create
##  one of your own with R's cat function.  Like this:
## myreg <-
## outreg(modelList = list(`Model 1` = m1, `Model 2` = m2), type = "html",
##     varLabels = varLabs)
##
##  cat(myreg, file = "reg.html")
##  Then open reg.html in a word processor or web browser.
``````