John Snow Project - DiD Regressions

Difference in Differences Regressions - Data from Snow 1855

See “Causality in the Time of Cholera” working paper at https://papers.ssrn.com/abstract=3262234 and my John Snow project website

This notebook is licensed under the BSD 2-Clause License

Introduction

This notebook performs linear and count regressions for 1849 versus 1854 data and builds on the discussion and data in the notebook “Snow1855_SimpleDiD_QRCT” - read that first.

For a brief introduction to Snow’s work, see:

• Snow’s original 1855 monograph (it is masterful): Snow, John. 1855. On the Mode of Communication of Cholera. 2nd ed. London: John Churchill. http://archive.org/details/b28985266.
• The best popular exposition I have found: Johnson, Steven. 2007. The Ghost Map: The Story of London’s Most Terrifying Epidemic–and How It Changed Science, Cities, and the Modern World. Reprint edition. New York: Riverhead Books.
• Another good popular version: Hempel, Sandra. 2007. The Strange Case of the Broad Street Pump: John Snow and the Mystery of Cholera. First edition. Berkeley: University of California Press.
• Tufte’s classic discussion of Snow’s mapping (a topic I don’t cover here): Tufte, Edward R. 1997. Visual Explanations: Images and Quantities, Evidence and Narrative. 1st edition. Graphics Press.
• Biography: Vinten-Johansen, Peter, Howard Brody, Nigel Paneth, Stephen Rachman, and Michael Russell Rip. 2003. Cholera, Chloroform and the Science of Medicine: A Life of John Snow. Oxford; New York: Oxford University Press. Linked on-line resources https://johnsnow.matrix.msu.edu/snowworks.php

This is an R Markdown Notebook. When you execute code within the notebook, the results appear beneath the code. The results are also saved in a self-contained html document with the suffix .nb.html. If you want pure r code (for example to run outside RStudio) you can easily extract code with the command knit(‘notebook.Rmd’,tangle=TRUE) which will save a file ‘notebook.R’ under your working directory.

Try executing the chunk below by clicking the Run button within the chunk or by placing your cursor inside it and pressing Cmd+Shift+Enter.

NB: We need the data and some of the calculations from the notebook ‘Snow1855_SimpleDiD_QRCT.Rmd’ so we execute this little hack to extract the R code from that sheet, write a new .R file, and then source that here

# Copyright (c) 2019, Thomas Coleman
#
#  -------  Licensed under BSD 2-Clause "Simplified" License  -------
#
# Results and discussion in "Causality in the Time of Cholera: John Snow as a Prototype 
# for Causal Inference (Working Paper)" available at SSRN: https://papers.ssrn.com/abstract=3262234
library(knitr)
options(scipen=5)
knit('Snow1855_SimpleDiD_QRCT.Rmd', tangle=TRUE)
source('Snow1855_SimpleDiD_QRCT.R') 
# The following libraries are used for the Negative Binomial regression and the robust standard error analysis
#install.packages("sandwich")
#install.packages("lmtest")
library("MASS")
library("sandwich") 
library("lmtest") 

Difference-in-Differences Framework

The basic difference-in-differences structure can be displayed in a table like:

Sub-districts	1849 Mortality (rate)	1854 Mortality (rate)	Diff 1854 less 1849
First 12 Southwark-supplied	\(\mu\)	\(\mu + \delta 54\)	\(\delta 54\)
Next 16 jointly-supplied	\(\mu + \gamma\)	\(\mu + \delta 54 + \gamma + \beta\)	\(\delta 54 + \beta\)
Diff next-16 less first-12	\(\gamma\)	\(\gamma + \beta\)	\(\beta\)

# The diff-in-diffs table (levels) from the notebook 'Snow1855_SimpleDiD_QRCT.Rmd'
kable(diff_in_diffs_Level,digits=2,caption="Diff-in-Diffs, Levels (rate per 10,000)",format='pandoc')

Diff-in-Diffs, Levels (rate per 10,000)

	level 1849	level 1854	col diffs
first 12 Southwark-only	134.86	146.61	11.75
next 16 Jointly-supplied	130.10	84.86	-45.24
row diffs	-4.76	-61.75	-56.99

kable(diff_in_diffs_Log,digits=3,caption="Diff-in-Diffs, Levels (rate per 10,000)",format='pandoc')

Diff-in-Diffs, Levels (rate per 10,000)

	log 1849	log 1854	col diffs
first 12 Southwark-only	-4.306	-4.223	0.084
next 16 Jointly-supplied	-4.342	-4.769	-0.427
row diffs	-0.036	-0.547	-0.511

This table looks very simple, and with a population of over 450,000 (167654 in the “first-12” and 300149 in the “next-16”) we would think the estimate of -0.511 would be precisely estimated. In fact, we can run a Poisson regrssion, a count regression which seems the right way to fit this table and test for the statistical significance of the diff-in-diffs coefficient.

regdatasimple <- data.frame("deaths" = c(comparecounts1849_1854[1:2,]), "pop1851" = c(comparepop1849_1854[1:2,]))
regdatasimple$supplier <- c("SouthwarkVauxhall","SouthwarkVauxhall_Lambeth","SouthwarkVauxhall","SouthwarkVauxhall_Lambeth")
regdatasimple$year <- c("1849","1849","1854","1854")
pois1simple <- glm(deaths ~ supplier * year + offset(log(pop1851)), family=poisson, data=regdatasimple) 
summary(pois1simple)

Call:
glm(formula = deaths ~ supplier * year + offset(log(pop1851)), 
    family = poisson, data = regdatasimple)

Deviance Residuals: 
[1]  0  0  0  0

Coefficients:
                                           Estimate Std. Error  z value Pr(>|z|)    
(Intercept)                                -4.30610    0.02103 -204.755  < 2e-16 ***
supplierSouthwarkVauxhall_Lambeth          -0.03593    0.02643   -1.359  0.17400    
year1854                                    0.08354    0.02914    2.867  0.00414 ** 
supplierSouthwarkVauxhall_Lambeth:year1854 -0.51088    0.03870  -13.201  < 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: 4.9099e+02  on 3  degrees of freedom
Residual deviance: 8.4466e-13  on 0  degrees of freedom
AIC: 46.995

Number of Fisher Scoring iterations: 2

The diff-in-diffs coffficient is supplierSouthwarkVauxhall_Lambeth:year1854 with the same value (in log terms) as the table above: -0.511 (in ratio terms - the increase in mortality for Southwark & Vauxhall customers versus Lambeth customers - this is 1.667). This appears to be very significant, with z value -13.2. But we would be wrong to rely on this regression.

We can see the reason when we examine the original data (Snow’s Table XII): there is considerable variation across sub-districts. Looking only at the first two sub-districts in 1849, for example, St. Saviour, Southwark has mortality of 143.6 while St. Olave, Southwark has mortality of 195.9. We cannot simply highlight and exploit the difference between the aggregate “first-12” and “next-16” sub-districts (the diff-in-diffs coefficient -0.5109) while ignoring the variability across and within sub-districts. We need some statistical framework to incorporate and analyze the variability both within and across sub-districts.

tablexii[,c(1,2,3,4,6,7)]

For a reasonable statistical framework we can consider the following equation, which puts the difference-in-differences table above into equation form but also allows multiple sub-districts - each sub-district having its own rate, etc.:

\(ln(Rate) = \mu + \delta 54*I(54) + \gamma*I(joint) + \beta*I(54)*I(joint) + \epsilon\)

• an overall constant (\(\mu\))
• a difference for 1854 (\(\delta54\))
• a difference for joint “next 16” region (\(\gamma\))
• an interaction for 1854 and joint (\(\beta\))

For various reasons we want to run this as linear-in-logs rather than linear-in-rates.

Formatting the regression data

To run this as a regression we need to take the death counts from Table XII (deaths by sub-district), combine with population from Table VIII for calculating rates, incorporate proper indicator variables, and then stack the data. Our original data (the .csv files read in above) have a variable “supplier” that acts as indicator for “first-12” (Southwark-only supply) versus “next-16” (jointly-supplied) sub-districts and another (“lambethdegree”) that serves as an indicator for “Southwark-only”, “less-Lambeth” and “more-Lambeth” supply. We need to construct a year factor (1849 versus 1854). The stacking an indicator construction is executed with the following code chunk.

x1 <- subset(tableviii,supplier == "SouthwarkVauxhall" | supplier == "SouthwarkVauxhall_Lambeth")
x1849 <- x1[c("subDistrict","pop1851","supplier","lambethdegree")]
x1 <- subset(tablexii,supplier == "SouthwarkVauxhall" | supplier == "SouthwarkVauxhall_Lambeth")
x1849$deaths <- x1$deaths1849
x1849$rate <- 10000 * x1$deaths1849 / x1849$pop1851
x1849$seq <- c(seq(1,length(x1849$deaths)))
#x1849$dum1854 <- 0
xyear <- factor(c(rep(1849,28),rep(1854,28)))
x1849$year <- xyear[1:28]
x1854 <- x1849
#x1849$lambethdegree <- "dirty"
x1854$deaths <- x1$deaths1854
x1854$rate <- 10000 * x1$deaths1854 / x1849$pop1851
x1854$seq <- c(seq(1,length(x1849$deaths)))
#x1854$dum1854 <- 1
x1854$year <- xyear[29:56]
regdata <- rbind(x1849,x1854)
regdata

Linear Regressions

We can start out by running a linear regression (linear in logs). The following code chunk executes three linear regressions but only displays the first: with Lambeth-supplier (“supplier” or “first-12” vs “next-16”) and year indicators.

# Linear-in-logs with Lambeth ("supplier" or "first-12" vs "next-16") and year indicators
lin1single <- lm(log(rate/10000) ~ supplier * year, data=regdata) 
lin1singlerobustse <- coeftest(lin1single, vcov = vcovHC(lin1single))
summary(lin1single)

Call:
lm(formula = log(rate/10000) ~ supplier * year, data = regdata)

Residuals:
    Min      1Q  Median      3Q     Max 
-1.9897 -0.1435  0.1538  0.3184  0.8443 

Coefficients:
                                           Estimate Std. Error t value Pr(>|t|)    
(Intercept)                                -4.50256    0.16643 -27.053   <2e-16 ***
supplierSouthwarkVauxhall_Lambeth           0.07770    0.22017   0.353   0.7256    
year1854                                    0.07392    0.23537   0.314   0.7547    
supplierSouthwarkVauxhall_Lambeth:year1854 -0.56665    0.31137  -1.820   0.0745 .  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.5765 on 52 degrees of freedom
Multiple R-squared:  0.1288,    Adjusted R-squared:  0.07852 
F-statistic: 2.562 on 3 and 52 DF,  p-value: 0.06473

# Linear-in-logs with two Lambeth indicators ("first-12" vs "next-16 split into less_lambeth" and "more_lambeth") and year indicators
lin1both <- lm(log(rate/10000) ~ lambethdegree * year, data=regdata) 
lin1bothrobustse <- coeftest(lin1both, vcov = vcovHC(lin1both))
# Linear-in-logs with sub-district fixed effects and Lambeth ("supplier" or "first-12" vs "next-16") and year indicators
lin2both <- lm(log(rate/10000) ~ subDistrict + lambethdegree * year, data=regdata) 
lin2bothrobustse <- coeftest(lin2both, vcov = vcovHC(lin2both))

The coefficient we are interested in is the last: supplierSouthwarkVauxhall_Lambeth:year1854: the coefficient \(\beta\) in the equation above. This calculates the average Lambeth effect, averaged across the sub-districts. The standard errors measure the precision of our estimated average based on the variation across sub-districts.

But here is an important puzzle: why is the coefficient -0.567 not the same as the calculated Lambeth effect in the table above (-0.511)? The reason is that the linear regressios are not correct, and this leads to our next topic.

Poisson Count Regressions

The linear-in-logs regression above looks correct but it is not: we have implicitely assumed that the error term \(\epsilon\) is normal and this cannot be correct. To see why, re-write the regression equation (using the fact that Rate = Count/Population so that ln(Rate) = ln(Count)-ln(Population)) as

\(ln(Count) = \mu + \delta 54*I(54) + \gamma*I(joint) + \beta*I(54)*I(joint) + \epsilon + ln(Population)\)

Count is the number of deaths and can only take non-negative integer values, so the error \(\epsilon\) cannot be Normal. It can, however, be Poisson. (For a discussion of the Poisson distribution and Poisson and Negative Binomial regression see my working paper at https://papers.ssrn.com/abstract=3262234 and the references therein.) There are various routines in R and other statstical packages for running regression assuming that the error \(\epsilon\) is Poisson (or other count distributions such as Negative Binomial).

The following code chunk runs Poisson regressions. The various combinations of single versus two Lambeth (treatment) effects and fixed effects are discussed in my working paper at https://papers.ssrn.com/abstract=3262234

# Poisson with single "Lambeth effect" and same rate for all sub-districts (no sub-district fixed effects)
pois1single <- glm(deaths ~ supplier * year 
    + offset(log(pop1851)), family=poisson, data=regdata) 
pois1singlerobustse <- coeftest(pois1single, vcov = vcovHC(pois1single))
#summary(pois1single)
#logLik(pois1single)
# Poisson with two "Lambeth effects" and same rate for all sub-districts (no sub-district fixed effects)
pois1both <- glm(deaths ~ lambethdegree * year 
    + offset(log(pop1851)), family=poisson, data=regdata) 
pois1bothrobustse <- coeftest(pois1both, vcov = vcovHC(pois1both))
#summary(pois1both)
#logLik(pois1both)
# Poisson with single "Lambeth effect" and different rates by sub-district (fixed effects)
pois2single <- glm(deaths ~ subDistrict + supplier * year 
    + offset(log(pop1851)), family=poisson, data=regdata) 
pois2singlerobustse <- coeftest(pois2single, vcov = vcovHC(pois2single))
#summary(pois2single)
#logLik(pois2single)
# Poisson with two "Lambeth effects" and different rates by sub-district (fixed effects)
pois2both <- glm(deaths ~ subDistrict + supplier * year 
    + offset(log(pop1851)), family=poisson, data=regdata) 
pois2bothrobustse <- coeftest(pois2both, vcov = vcovHC(pois2both))
#summary(pois2both)
#logLik(pois2both)
# NB: The robust standard errors approach I found at:
#  https://stat.ethz.ch/pipermail/r-help/2008-May/161591.html
#Then click on "Next message:" for the second. Unfortunately, 
#at that point the thread broke 
# So to continue, next take: 
#  https://stat.ethz.ch/pipermail/r-help/2008-May/161640.html
#and thereafter continue to click on "Next message:" until the 
#thread runs out. 
# Also cf https://stats.stackexchange.com/questions/117052/replicating-statas-robust-option-in-r
# NB: STATA seems to use HC1 for robust BUT with n-1 instead of n-k (see "sandwich.pdf" and 
#   https://stats.stackexchange.com/questions/89999/how-to-replicate-statas-robust-binomial-glm-for-proportion-data-in-r)
# Display the regression for the single Lambeth effect here:
summary(pois1single)

Call:
glm(formula = deaths ~ supplier * year + offset(log(pop1851)), 
    family = poisson, data = regdata)

Deviance Residuals: 
     Min        1Q    Median        3Q       Max  
-12.1542   -3.0953    0.0306    2.7398   10.2518  

Coefficients:
                                           Estimate Std. Error  z value Pr(>|z|)    
(Intercept)                                -4.30610    0.02103 -204.755  < 2e-16 ***
supplierSouthwarkVauxhall_Lambeth          -0.03593    0.02643   -1.359  0.17400    
year1854                                    0.08354    0.02914    2.867  0.00414 ** 
supplierSouthwarkVauxhall_Lambeth:year1854 -0.51088    0.03870  -13.201  < 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: 2032.6  on 55  degrees of freedom
Residual deviance: 1541.6  on 52  degrees of freedom
AIC: 1937.4

Number of Fisher Scoring iterations: 4

This regession is for a single Lambeth effect. Now the coefficient for the Lambeth effect (-0.511) matches exactly the value for \(\beta\) calculated in the table above.

It appears that this estimated Lambeth or treatment effect is very highly significant: z value = -13.2. But this is wrong. This is discussed in more detail in another notebook (and my working paper), but suffice to say here that the Poisson regression does not fit the data very well. In the statistics literature this is termed “overdispersion”. The diagnostic tool we can use is the “Residual Deviance”. It is approximately chi-squared distributed and the probability of observing a value as large as we see here (1542) is tiny.

The result of the overdispersion for the Poisson regression is that the true standard errors are much larger than reported in the regression. One way to handle this is to use robust standard errors. The code chunk above calculates robust standard errors and the comments in the code chunk have links to more discussion.

Negative Binomial Count Regressions

Another way to handle the overdispersion (again, discussed in more detail in another notebook and my working paper) is to allow the counts (specifically the error term \(\epsilon\)) to be Negative Binomial distributed. Once again, R has routines for this, in the MASS package.

# Negative Binomial with single "Lambeth effect" 
nb1single <- glm.nb(deaths ~ supplier * year 
    + offset(log(pop1851)), data=regdata) 
nb1singlerobustse <- coeftest(nb1single, vcov = vcovHC(nb1single))
#summary(nb1single)
#logLik(nb1single)
#print(coeftest(nb1single, vcov = vcovHC(nb1single)))
# Negative Binomial with two "Lambeth effects" 
nb1both <- glm.nb(deaths ~ lambethdegree * year 
    + offset(log(pop1851)), data=regdata) 
nb1bothrobustse <- coeftest(nb1both, vcov = vcovHC(nb1both))
#summary(nb1both)
#logLik(nb1both)
#print(coeftest(nb1both, vcov = vcovHC(nb1both)))
# Negative Binomial with two "Lambeth effects" and sub-district fixed effects
nb2both <- glm.nb(deaths ~ subDistrict + lambethdegree * year 
    + offset(log(pop1851)), data=regdata) 
nb2bothrobustse <- coeftest(nb2both, vcov = vcovHC(nb2both))
#summary(nb2both)
#logLik(nb2both)
#print(coeftest(nb2both, vcov = vcovHC(nb2both)))
# Show results for the Negative Binomial with one Lambeth effect
summary(nb1single)

Call:
glm.nb(formula = deaths ~ supplier * year + offset(log(pop1851)), 
    data = regdata, init.theta = 4.956204025, link = log)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-3.2718  -0.5282   0.0652   0.4958   1.8326  

Coefficients:
                                           Estimate Std. Error z value Pr(>|z|)    
(Intercept)                                -4.33237    0.13169 -32.898   <2e-16 ***
supplierSouthwarkVauxhall_Lambeth          -0.03182    0.17386  -0.183   0.8548    
year1854                                    0.05734    0.18616   0.308   0.7581    
supplierSouthwarkVauxhall_Lambeth:year1854 -0.50027    0.24612  -2.033   0.0421 *  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for Negative Binomial(4.9562) family taken to be 1)

    Null deviance: 71.825  on 55  degrees of freedom
Residual deviance: 59.764  on 52  degrees of freedom
AIC: 657.19

Number of Fisher Scoring iterations: 1

              Theta:  4.956 
          Std. Err.:  0.973 

 2 x log-likelihood:  -647.185 

The negative binomial is reasonably good at fitting the observed sub-district variability: the Residual Deviance is now 59.8 which is small enough to not reject this regression. As a result we can have some confidence that the standard errors reported are realistic. For this regression the single Lambeth effect is modestly significant, with a z value of -2.03.

Two “Lambeth effects”

Snow noted that four of the “next-16” jointly-supplied sub-districts had a higher proportion of Lambeth customers than the others: “In certain sub-districts, where I know that the supply of the Lambeth Water Company is more general than elsewhere, as Christchurch, London Road, Waterloo Road 1st, and Lambeth Church 1st, the decrease of mortality in 1854 as compared with 1849 is greatest, as might be expected.” (Snow 1855 p. 89) We can incorporate this by having two indicator variables, the variable “lambethdegree” with “less Lambeth” and “more Lambeth”.

# Show results for the Negative Binomial with two Lambeth effects
summary(nb1both)

Call:
glm.nb(formula = deaths ~ lambethdegree * year + offset(log(pop1851)), 
    data = regdata, init.theta = 5.57218113, link = log)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-3.4384  -0.5744   0.0754   0.5203   1.5608  

Coefficients:
                                   Estimate Std. Error z value Pr(>|z|)    
(Intercept)                        -4.33210    0.12443 -34.815  < 2e-16 ***
lambethdegreeless_Lambeth          -0.06426    0.17551  -0.366  0.71427    
lambethdegreemore_Lambeth           0.05911    0.24787   0.238  0.81151    
year1854                            0.05740    0.17589   0.326  0.74415    
lambethdegreeless_Lambeth:year1854 -0.33838    0.24839  -1.362  0.17310    
lambethdegreemore_Lambeth:year1854 -1.13213    0.35350  -3.203  0.00136 ** 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for Negative Binomial(5.5722) family taken to be 1)

    Null deviance: 80.122  on 55  degrees of freedom
Residual deviance: 59.991  on 50  degrees of freedom
AIC: 654.94

Number of Fisher Scoring iterations: 1

              Theta:  5.57 
          Std. Err.:  1.11 

 2 x log-likelihood:  -640.943 

The effect for “more Lambeth” is large (-1.132, ratio effect 3.102) and highly significant (z value -3.203). When we introduce the population data published in 1856 (another notebook and my working paper) we find an even stronger Lambeth effect.

Creating Tables for Results

To summarize the results I create tables in the following code chunk to summarize the regressions. I don’t display these tables but you can if you wish.

# Regression table for linear, with robust SEs for poisson
regtablelin <- matrix(0,nrow=14,ncol=5)
colnames(regtablelin) <- c("Linear Single","Poiss1 Single robust","NB1 Single","Linear Both","NB1 Both")
rownames(regtablelin) <- c("Treat less Lamb","SE", "z value","p-value","treat (ratio)",
    "Treat more Lamb","SE","z value", "p-value","treat (ratio)","theta","region1","region2","time")
# Populate the table from the regressions
regtablelin[c(1,2,3,4),1] <- summary(lin1single)$coefficients[4,c(1,2,3,4)]
regtablelin[c(1,2,3,4),2] <- summary(pois1single)$coefficients[4,c(1,2,3,4)]
regtablelin[c(1,2,3,4),3] <- summary(nb1single)$coefficients[4,c(1,2,3,4)]
regtablelin[c(1,2,3,4),4] <- summary(lin1both)$coefficients[5,c(1,2,3,4)]
regtablelin[c(6,7,8,9),4] <- summary(lin1both)$coefficients[6,c(1,2,3,4)]
regtablelin[c(1,2,3,4),5] <- summary(nb1both)$coefficients[5,c(1,2,3,4)]
regtablelin[c(6,7,8,9),5] <- summary(nb1both)$coefficients[6,c(1,2,3,4)]
regtablelin[11,3] <- nb1single$theta
regtablelin[11,5] <- nb1both$theta
# Calculate the treatment effect as a ratio (exponentiate)
regtablelin[5,] <- exp(-regtablelin[1,])
regtablelin[10,c(4,5)] <- exp(-regtablelin[6,c(4,5)])
# Now the control (region & time) effects
regtablelin[c(12,14),1] <- lin1single$coefficients[c(2,3)]
regtablelin[c(12,14),2] <- pois1singlerobustse[c(2,3),1]
regtablelin[c(12,14),3] <- nb1single$coefficients[c(2,3)]
regtablelin[c(12,13,14),4] <- lin1both$coefficients[c(2,3,4)]
regtablelin[c(12,13,14),5] <- nb1both$coefficients[c(2,3,4)]
# Create table with regression results (diff-in-diffs estimate) using robust SEs
# for all except the NB1 (Negative Binomial without sub-district fixed effects)
regtablepoiss <- matrix(0,nrow=19,ncol=5)
colnames(regtablepoiss) <- c("Poiss Single","Poiss Single FE","NB1 Single","NB1 Both","NB2 Both")
rownames(regtablepoiss) <- c("Treat less Lamb","SE","z value","p-value",
    "robust SE","z value","treat (ratio)",
    "Treat more Lamb","SE","z value","p-value",
    "robust SE","z value","treat (ratio)",
    "theta","resid dev","region1","region2","time")
# Populate the table from the regressions
regtablepoiss[c(1,2,3,4),1] <- summary(pois1single)$coefficients[4,c(1,2,3,4)] # coeff, SE, t-ratio, p-value
regtablepoiss[c(5,6),1] <- pois1singlerobustse[4,c(2,3)]             # robust SE, t-ratio
regtablepoiss[16,1] <- pois1single$deviance
# Now the control (region & time) effects
regtablepoiss[c(17,19),1] <- summary(pois1single)$coefficients[c(2,3),1] 
regtablepoiss[c(1,2,3,4),2] <- summary(pois2single)$coefficients[30,c(1,2,3,4)]  # Lambeth effect is at end
regtablepoiss[c(5,6),2] <- pois2singlerobustse[30,c(2,3)]     
regtablepoiss[16,2] <- pois2single$deviance
regtablepoiss[19,2] <- summary(pois2single)$coefficients[29,1]    # time effect
regtablepoiss[c(1,2,3,4),3] <- summary(nb1single)$coefficients[4,c(1,2,3,4)]
regtablepoiss[c(5,6),3] <- nb1singlerobustse[4,c(2,3)]
regtablepoiss[15,3] <- nb1single$theta
regtablepoiss[16,3] <- nb1single$deviance
regtablepoiss[c(17,19),3] <- nb1single$coefficients[c(2,3)]
regtablepoiss[c(1,2,3,4),4] <- summary(nb1both)$coefficients[5,c(1,2,3,4)] # coeff, SE, t-ratio, p-value for 1st Lambeth effect
regtablepoiss[c(5,6),4] <- nb1bothrobustse[5,c(2,3)]
regtablepoiss[c(8,9,10,11),4] <- summary(nb1both)$coefficients[6,c(1,2,3,4)] # coeff, SE, t-ratio, p-value for 2nd Lambeth effect
regtablepoiss[c(12,13),4] <- nb1bothrobustse[6,c(2,3)]
regtablepoiss[15,4] <- nb1both$theta
regtablepoiss[16,4] <- nb1both$deviance
regtablepoiss[c(17,18,19),4] <- nb1both$coefficients[c(2,3,4)]
regtablepoiss[c(1,2,3,4),5] <- summary(nb2both)$coefficients[30,c(1,2,3,4)]
regtablepoiss[c(5,6),5] <- nb2bothrobustse[30,c(2,3)]
regtablepoiss[c(8,9,10,11),5] <- summary(nb2both)$coefficients[31,c(1,2,3,4)]
regtablepoiss[c(12,13),5] <- nb2bothrobustse[31,c(2,3)]
regtablepoiss[15,5] <- nb2both$theta
regtablepoiss[16,5] <- nb2both$deviance
regtablepoiss[19,5] <- summary(nb2both)$coefficients[29]
# Calculate the treatment effect as a ratio (exponentiate)
regtablepoiss[7,] <- exp(-regtablepoiss[1,])
regtablepoiss[14,c(4,5)] <- exp(-regtablepoiss[8,c(4,5)])
regtablelin <- as.data.frame(regtablelin)
regtablepoiss <- as.data.frame(regtablepoiss)
#regtablelin
#regtablepoiss

---

---

Digression on Binomial, Poisson, and Negative Binomial

For more detail on count regressions see my working paper “Causality in the Time of Cholera” working paper at https://papers.ssrn.com/abstract=3262234 and references there. But here are some notes

For binomial:

• Rate = Count / n = Mean / n = p
• Mean = count = n*p
• Variance(count) = n*p*(1-p) => Std dev = sqrt(n*p*(1-p)) ~ sqrt(count) = sqrt(n*p) (since p is small for this application)
• Std dev (rate) ~ sqrt(count)/n = sqrt(p)*sqrt(n)/n = sqrt(p) / sqrt(n)

For Poisson (which is good approx to binomial for low intensities):

• Rate = count / n = rate
• Mean(count) = count = n*rate
• Variance(count) = n*rate => Std dev = sqrt(n*rate) = sqrt(count)
• Std dev (rate) = sqrt(count)/n = sqrt(rate)*sqrt(n)/n = sqrt(rate) / sqrt(n)

Mixing with Gamma(k,k):

• Gamma(a,b): Mean=a/b, Var=a/b^2
• Gamma(k,k): Mean=1, Var=1/k
• Poisson(mu) mixed with Gamma(k,k) -> Mean = mu & variance mu/k or mu + mu^2/k ?

It looks like for the R glm.nb the variance is mu + mu^2/theta ?? I think this is

• mean = rate * pop
• variance = mean + mean^2 / theta = rate*pop + (rate*pop)^2/theta
• => stdev of rate = sqrt(variance(count)) / pop = sqrt(rate/pop + rate^2/theta)

Example:

• Estimate neg binom const = -4.35 => rate=exp(-4.35)=.012907, theta = 5.51
• For pop = 20,000:
• Poisson stdev(rate) = sqrt(.0129/20000) = .0008 (8 per 10,000)
• Neg binom stdev(rate) = sqrt(.0129/20000 + .0129^2/5) = .005828 (58 per 10,000)