```{r setup, include=FALSE}
# set global chunk options
library(knitr)
opts_chunk$set(out.width='600px', dev.args=list(bg='white'), comment=NA, tidy=FALSE)
options(digits=3)
```

Female vs Male Hurricanes --- Commentary
====================================================================
John Maindonald
---------------

Dec 26 2014.

> The notes that follow comment on analyses that are reported in 
> Jung et al (2014), and in particular on the analysis that leads 
> to their Figure 1. Diagnostic plots reveal serious problems with 
> their analysis. Use of the Box-Cox methodology to identify a
> suitable power transformation leads to a model that fits the
> data well. It has `NDAM`$^{0.13}$ as the only explanatary variable,
> with addition of terms involving the pressure at landfall and
> the masculinity-feminity index giving no improvement to the fit.
> Diagnostics are very satisfactory, with no outliers or highly
> influential points identified. For the modeling reported here, 
> highly crucial roles attach both to preliminary graphical 
> exploration and to careful examination of model diagnostics.

> A further analysis works with the updated ICAT damage estimates
> that are in 2014 US dollars, and investigates inclusion of Audrey 
> and Katrina that were omitted from the Jung et al analysis.
> Also examined is the model fit to data from 1979 onwards, when 
> male and female names were alternated.

> A further comment is that the regression coefficients are not
> capable of sustaining the interpretation that Jung et al wish 
> to place on them.

> Aside from the main theme are comments on the variety of functions,
> in different R packages, that are available for the fitting of 
> regression models for data with a negative binomial error.

> The data for Audrey and Katrina is corrected relative to that used
> in the June 20 2014 version of this document.

### Preparatory Steps
#### Download and input data
```{r get-data}
fnam <- paste0("http://www.pnas.org/lookup/suppl/doi:10.1073/",
              "pnas.1402786111/-/DCSupplemental/pnas.1402786111.sd01.xlsx")
if(!file.exists('hurricane-names.xlsx')){
download.file(fnam, destfile='hurricane-names.xlsx', mode="wb")
}
library(xlsx)
hurricJ <- read.xlsx('~/r/_rn/Data/hurricane-names.xlsx', 
                         sheetIndex=1, rowIndex=1:93)
rnam <- with(hurricJ, paste0(as.character(Name), Year))
row.names(hurricJ) <- rnam
hurricJ[,"Gender_MF"] <- factor(hurricJ[,"Gender_MF"], labels=c("m","f"))
```

#### Check variables
keyVar <- hurricJ[, c("MasFem","alldeaths","MinPressure_before","Category","ZNDAM","Year")]
sapply(keyVar,mean)

#### Examine plots
Use a `log(alldeaths+0.5)` scale on the y-axis in order to separate points.
Smooth curves have been added, separately for hurricanes with male and
female names.  
```{r plots, fig.width=8, fig.height=4.25, out.width="600px", fig.cap="Figure 1: Plot of log(alldeaths+0.5) against log(NDAM)"}
library(lattice)
ypos <- seq(from=0, to=200, by=50)
sc <- list(x=list(relation="free"), y=list(at=sqrt(ypos), labels=paste(ypos)))
xyplot(log(alldeaths+0.5) ~ NDAM+log(NDAM), groups=Gender_MF, xlab="",
       par.settings=simpleTheme(pch=c(15,16),cex=1.1), data=hurricJ, 
       scale=sc, auto.key=list(columns=2), type=c('p',"smooth"))
```

Use of a logarithmic x-axis scale ameliorates severe leverage issues.
It also gives variables that are closer to linearly related.  Curvsture
in the smooths, at least for female names, suggests that a logarithmic 
transformation may be too extreme.  

Now check the range of values of `MinPressure_before`. 
```{r MinPressure_before-range}
with(hurricJ, range(MinPressure_before))
```
The ratio is 1.1 to 1, too small for a logarithmic or power transformation
to make much difference.

Now (Figure 2) plot against `MinPressure_before`
```{r MinPressure_before, fig.width=4.25, fig.height=4.25, out.width="350px", fig.cap="Figure 2: Plot of log(alldeaths+0.5) against MinPressure_before"}
xyplot(log(alldeaths+0.5) ~ MinPressure_before, 
       groups=Gender_MF, data=hurricJ, 
       par.settings=simpleTheme(pch=c(15,16)),
       scale=sc, auto.key=list(columns=2), type=c('p',"smooth"))
```

### Negative Binomial versus `lm` style Linear Models
Once the counts are large enough, there is little to distinguish a
model that uses `log(alldeaths+0.5)` as the outcpme variabke from a
Negative Binomial model that has `var[deaths]` proportional to 
`E[deaths]`$^2$.  For these data, with 41 of the counts 4 or less,
the approximation is rough.  It may nevertheless be useful for
exploratory purposes.

The following fits the preferred Jung at al model.
```{r as-JungEtAl}
library(MASS)
hurricJ0.nb <- glm.nb(alldeaths ~ 1, data=hurricJ)
hurricJ3.nb <- glm.nb(alldeaths ~ 
                    MasFem*(MinPressure_before+NDAM), 
                    data=hurricJ)
## MasFem:NDAM indicates that the product of MasFem and NDAM 
## is to be included as a variable in the model.
```

Now compare the two models.
```{r Jung-anova}
anova(hurricJ0.nb, hurricJ3.nb)
## Pearson chi-squared for Jung et al preferred model
sum(resid(hurricJ3.nb, type='pearson')^2)/hurricJ3.nb[["df.residual"]]
## 
## Check model coefficients
coef(summary(hurricJ3.nb))
```

Figure 3 shows the diagnostic plots.
```{r as-JungEtAl-diag, fig.width=7, fig.height=7, out.width="500px", fig.cap="Figure 3: Diagnostic plots for the preferred Jung et al model"}
opar <- par(mar=c(3.5,3.5,2.1,0.6), mgp=c(2.25,0.5,0), mfrow=c(2,2))
plot(hurricJ3.nb)
par(opar)
```

A small number of data points have very high leverage, relative 
to other points.  While Diane is the only hurricane that triggers
the standard Cook's distance criterion, there may well be
combined effects that are distorting the model.

Notice that the coefficient of `NDAM` is negative.  The 
interaction with `MasFem` is accounting for a major part of the
effect of `NDAM`.

Figures 1 and 3 are together good reason not to go with the Jung et al model.
Reasonable alternatives are then to work with `log(NDAM)`, or with a power
transformation of `NDAM`.  The right panel in Figure 1 suggests a
power transformation.

## Use of `NDAM`$^\lambda$ as explanatory variable

For checking for the most appropriate transformation, Cook (1998) argues
for examining the regressions of powers of the explanatory variable 
on the outcome variable.  This reverses the roles of NDAM` 
and `log(deaths+0.5)`, relative to their roles in the model fit. 
See further Weisberg (2014), Fox and Weisberg (2014), and
Maindonald and Braun (2012),

The following uses the function `powerTransform()` in the {\em car} package to
find a suitable transformation.  This assumes, of course, that `log(deaths+0.5)`
is a reasonable choice for dependent variable.  As already noted, this is roughly
equivalent to use of a log link function in a negative binomial model.
```{r powerTransform2}
library(car)
powerTransform(NDAM ~ log(alldeaths+0.5), data=hurricJ)
```
The function `powerTransform()` chooses a transformation so that
residuals are as close as possible to normal.
The same calculation with `Minpressure_before`
in place of `log(deaths+0.5)` suggests a power of 0.129, i.e., for
all practical purposes the same.

Working then with `NDAM`$^{0.13}$, we have:
```{r use-power0.13}
## Negative binomial
hurric.nb <- glm.nb(alldeaths ~ I(NDAM^0.13),
                     data=hurricJ)
coef(summary(hurric.nb))
## Calculate Pearson chi-squared
sum(resid(hurric.nb, type='pearson')^2)/hurric.nb[["df.residual"]]
## Roughly equivalent lm model, for comparison
hurric.lm <- lm(log(alldeaths+0.5) ~ I(NDAM^0.13),
                    data=hurricJ)
coef(summary(hurric.lm))
```
The estimates are about as similar as is reasonable to expect.
The diagnostic plots for the negative binomial model are:
```{r nb-diagnostics, fig.width=7, fig.height=7, out.width="500px", fig.cap="Figure 4: Diagnostic plots for the model that has NDAM$^{0.13}$ as the only explanatory variable"}
opar <- par(mar=c(3.5,3.5,2.1,0.6), mgp=c(2.25,0.5,0), mfrow=c(2,2))
plot(hurric.nb)
par(opar)
```

### Is there a case for adding further terms?

```{r try-further}
library(MASS)
hurric2.nb <- glm.nb(alldeaths ~ (I(NDAM^0.13)+MasFem)^2+MinPressure_before+
                 MasFem:MinPressure_before, maxit=100, 
                 data=hurricJ)
anova(hurric.nb, hurric2.nb)
```
Adding further terms to the model that has `NDAM`$^{0.13}$ as explanatory 
variable does not improve the explanatory power.

The plot method gives the following:
```{r plot-diags, fig.width=7, fig.height=7, out.width="500px"}
opar <- par(mar=c(3.5,3.5,2.1,0.6), mgp=c(2.25,0.5,0), mfrow=c(2,2))
plot(hurric.nb)
par(opar)
par(mfrow=c(1,1))
```

There are no outliers, and no highly influential points.

## Re-inclusion of Audrey (1957) and Katrina (2005)
Note that the following corrects mistakes in the figures for Audrey and
Katrina that were used in the June 20 2014 version of this document.

Wih the more acceptable model, can we now re-include the omitted hurricanes?
Unfortunately Jung et al left Audrey and Katrina out of their dataset. 
(It is good practice to report data for outliers.  They may turn out not to 
be outliers when a different model is used.  Or they may hold important 
information.)  

Data that included Audrey and Katrina, with damage estimates in 2014 dollars, 
were therefore obtained from http://www.icatdamageestimator.com/commonsearch .
NOAA Monthly Weather Reports (MWRs) supplied the numbers of deaths for most 
of the hurricanes in Jung et al's data.  Wikipedia wiki pages for individual
hurricanes provide what will often be better attested data, certainly for
more damaging hurricanes where time is required to pull together the evidence.
For present purposes, the number of deaths data from Jung et al's data set 
was supplemented with the numbers given on the wiki pages for 
Hurricane_Katrina and Hurricane_Audrey.

```{r bring-back}
if(!file.exists('hurric2014.csv')){
url <- "http://maths-people.anu.edu.au/~johnm/stats-issues/data/hurric2014.csv"
download.file(url, destfile="hurric2014.csv")
}
hurric2014 <- read.csv('hurric2014.csv', row.names=1)
hurric2014 <- within(hurric2014, Gender_MF <- factor(Gender_MF, levels=c("m","f")))
```

Figure 5 shows a scatterplot of `log(alldeaths+0.5)` against `log(NDAM2014)`
```{r plot-all, fig.width=7, fig.height=7, out.width="350px", fig.cap="Figure 5: Plot of log(alldeaths+0.5) against log(NDAM2014)"}
library(latticeExtra)
xpos <- c(1,200, 4*10^4, 8*10^6)
ypos <- c(0, 2, 6, 24, 100, 400, 1600)
sc <- list(y=list(at=log(ypos+0.5), labels=ypos),
           x=list(at=log(xpos), labels=xpos))
gph <- xyplot(log(alldeaths+0.5) ~ log(NDAM2014), groups=Gender_MF,
              xlab="Damage to property ($000,000)",
              ylab="Total deaths",
              par.settings=simpleTheme(pch=c(15,16),cex=1.35),
              data=hurric2014, scales=sc, auto.key=list(columns=2))
hnames <- c("Audrey1957", "Katrina2005")
nlabs <- match(hnames, row.names(hurric2014))
gph+layer(panel.text(x[nlabs], y[nlabs], labels=row.names(hurric2014)[nlabs], pos=2, ...))
```

### Check the appropriate power transformation
The following again uses the function `powerTransform()` in the {\em car} package
to check the transformation that might be used:
```{r powerTransform}
powerTransform(NDAM2014 ~ log(alldeaths+0.5), data=hurric2014)
```
As the change from 0.13 to 0.15 is of no consequence,
we will then work with `NDAM2014`$^{0.13}$, as before.

```{r use2-power0.13}
## Negative binomial
hurric2014.nb <- glm.nb(alldeaths ~ I(NDAM2014^0.13),
                     data=hurric2014)
coef(summary(hurric2014.nb))
## lm model as alternative
hurric2014.lm <- lm(log(alldeaths+0.5) ~ I(NDAM2014^0.13),
                    data=hurric2014)
coef(summary(hurric2014.lm))
```
The estimates are about as similar as is reasonable to expect.
The diagnostic plots for the negative binomial model are:
```{r nb-diagnostics2, fig.width=7, fig.height=7, out.width="500px", fig.cap="Figure 6: Diagnostic plots, model with Audrey and Katrina included."}
opar <- par(mar=c(3.5,3.5,2.1,0.6), mgp=c(2.25,0.5,0), mfrow=c(2,2))
plot(hurric2014.nb)
par(opar)
```
Audrey stands out very clearly as an outlier and as a high influence point.
Omitting Audrey gives:
```{r use-power0.13-omitAudrey, fig.width=7, fig.height=7, out.width="500px", fig.cap="Figure 7: Diagnostic plots, with Audrey excluded."}
## Negative binomial
hurric2014a.nb <- glm.nb(alldeaths ~ I(NDAM2014^0.13),
                     data=hurric2014[-13,])
coef(summary(hurric2014a.nb))
opar <- par(mar=c(3.5,3.5,2.1,0.6), mgp=c(2.25,0.5,0), mfrow=c(2,2))
plot(hurric2014a.nb)
par(opar)
```
Katrina now appears as a high influence point.

## Restriction to hurricanes from 1979 or later
Until 1979, all names used were female.  Since then, male and female 
names have alternated.  The following limits the analysis to data 
from 1979 onwards:

```{r fit-1979on}
hurric1979on.nb <- glm.nb(alldeaths ~ I(NDAM2014^0.13),
                     data=subset(hurric2014, Year>=1979))
coef(summary(hurric1979on.nb))
hurric1979onALL.nb <- glm.nb(alldeaths ~ I(NDAM2014^0.13)+
                      Minpressure_Updated.2014+
                    Gender_MF+ Gender_MF:log(NDAM2014)+ 
                      Gender_MF:Minpressure_Updated.2014,
                    data=subset(hurric2014, Year>=1979), maxit=100)
anova(hurric1979on.nb, hurric1979onALL.nb)
```

Now examine the diagnostics:
```{r lm-diag1979on, fig.width=7, fig.height=7, out.width="500px", fig.cap="Figure 8: Diagnostic plots, data from 1979 only"}
opar <- par(mar=c(3.5,3.5,2.1,0.6), mgp=c(2.25,0.5,0), mfrow=c(2,2))
plot(hurric1979on.nb)
par(opar)
```

Katrina has a large influence on the fit.  Hence, omit Katrina.
```{r fit-1979onOmitK, fig.width=7, fig.height=7, out.width="500px", fig.cap="Figure 9: Diagnostic plots from 1979 on, model that excludes Katrina"}
hurric1979onOmitK.nb <- glm.nb(alldeaths ~ I(NDAM2014^0.13),
                     data=subset(hurric2014[-84,], Year>=1979))
coef(summary(hurric1979onOmitK.nb))
## Now examine the diagnostics:
opar <- par(mar=c(3.5,3.5,2.1,0.6), mgp=c(2.25,0.5,0), mfrow=c(2,2))
plot(hurric1979onOmitK.nb)
par(opar)
```

The only point of note is a clear suggestion of increasing variance,
beyond what the model accommodates, with increasing fitted values.
A model that adequately accommodates data from 1979 on perhaps still
has to be found!  With a more appropriate model, perhaps Katrina would 
not be a point of high influence!

## Other Aspects of the Methodology
Jung et al's argument, as I understand it, is that the variables `NDAM`
(property damage in 2013 monetary values) and `MinPressure_before`
together account for the potential of hurricanes to cause human
deaths, and that `MasFem`, perhaps in interaction with `NDAM` and/or
`MinPressure_before`, then accounts for the difference between
potential deaths and observed deaths. The implied assumption that
`NDAM` and `MinPressure_before` are together an effective proxy for
the potential to cause human deaths may then be regarded with
scepticism.

Some damage to property, in some cases a large part. will have
occurred away from urban centres.  In effect, then, the analysis is
missing a variable that measures the extent of mismatch between the
combined effect of `NDAM` and `MinPressure_before` and the potential
for human damage.  Yule-Simpson type effects, involving continuous as
well as categorical variables, can in these circumstances bias or
induce effects in variables that are included in the model.  With
observational data, omission of a variable that has a substantial
effect will commonly induce such bias.  Look for example under the
index entries for *Lord's paradox* in Tu and Gilthorpe (2012).

An analysis of this data that passes standard checks has not
shown an effect from `MasFem` or from interactions with `MasFem`. The
issue of what interpretations can be placed on the regression coefficents
does not then arise, but remains a damaging criticism of the methodology.
See also the points raised in comments at the urls:

http://scatter.wordpress.com/2014/06/03/my-thoughts-on-that-hurricane-study/.

http://andrewgelman.com/2014/06/06/hurricanes-vs-himmicanes/

**Note:** This document was created from an R Markdown markup file.  The 
file is available as:
http://maths-people.anu.edu.au/~johnm/stats-issues/hurricnames.Rmd .
Under RStudio, click on "File", then on "Open File", and select the file.
Then click on "Knit HTML" to recreate this document.

## Additional Notes

### A dispersed Poisson model

This fits without problem, but fails the diagnostic tests!  This
model assumes a variance that is proportional to `E[deaths]`
rather than to `E[deaths]`$^2$.

```{r do-glm--disp, fig.width=7, fig.height=7, out.width="500px"}
library(MASS)
hurric.glm <- glm(alldeaths ~ log(NDAM)+MinPressure_before+MasFem+
                 I(MasFem*log(NDAM))+I(MasFem*MinPressure_before),
                 data=hurricJ, family=quasipoisson)
summary(hurric.glm)
opar <- par(mar=c(3.5,3.5,2.1,0.6), mgp=c(2.25,0.5,0), mfrow=c(2,2))
plot(hurric.glm)
par(opar)
```

The linear downward trend in the red line in the "Residuals vs Fitted"
plot presumably arises because it is a "smooth" that has been fitted
using a smoothing approach that is appropriate for a normal errors
model.  The scale-location plot is the big concern. The deviance
residuals adjust for a variance proportional to mean increase in
the variance.  Here, mean/sqrt(variance) increases very markedly with 
the mean.  The SEs from the negative binomial model are accordingly
suspect. 


### Functions in R Packages that fit negative binomial models

The above has benefitted from discussion on the ANZSTAT email list.
Negative binomial regression provides a somewhat extreme example of 
overlap between R packages, with notation that is not consistent 
between these different implementations.

In addition to `MASS::glm.nb()`, note `msme::nbinomial()`, 
`aod::negbin()` and `gamlss::gamlss()`. Mark Donaghue noted that 
`gamlss()` in the *gamlss* package fits two different types of 
NB model, either `family = NBI` (quadratic; var = mu(1 + sigma * mu)) 
as I think for all the functions above, or `family=NBII` 
(linear; var = mu(1 + sigma)).  

Note that the notation is not consistent across the different packages
and functions that allow a negative binomial error structure.

## References

> Cook, R (1998). Regression Graphics.  Wiley.

> Fox and Weisberg (2nd edn, 2014). An R Companion to Applied Regression. Sage.

> Jung, Shavitta, Viswanathana, and Hilbe (2013). 
> "Female hurricanes are deadlier than male hurricanes". 
> http://www.pnas.org/cgi/doi/10.1073/pnas.1402786111.

> Maindonald and Braun (3rd edn, 2010). Data Analysis and Graphics Using R. CUP.

> Tu and Gilthorpe (2012). "Statistical Thinking in Epidemiology".  CRC Press.

> Weisberg, S (4th edn, 2014).  "Applied Linear Regression". Wiley.



