Day 8: Panel Models for Limited Outcomes

panel data
code
analysis
Author

Robert W. Walker

Published

August 17, 2022

Slides

Two core issues remain: the spatial and the temporal. There are fewer well worked diagnostics for this than for the pure linear case owing to partial observability and/or coarse measurement in some fashion or other.

Panel GLM Computing

R

All of the computing will take place in the confines of pglm, a package in the spirit of plm for linear models. The help file provides an example of the familiar models. They are specified to match up with the syntax of R’s glm functions with family and link arguments.

# install.packages("pglm")
library(pglm)
## an ordered probit example
data('Fairness', package = 'pglm')
Parking <- subset(Fairness, good == 'parking')
op <- pglm(as.numeric(answer) ~ education + rule,
Parking[1:105, ],
family = ordinal('probit'), R = 5, print.level = 3,
method = 'bfgs', index = 'id', model = "random")
## a binomial (probit) example
data('UnionWage', package = 'pglm')
anb <- pglm(union ~ wage + exper + rural, UnionWage, family = binomial('probit'),
model = "pooling", method = "bfgs", print.level = 3, R = 5)
## a gaussian example on unbalanced panel data
data(Hedonic, package = "plm")
ra <- pglm(mv ~ crim + zn + indus + nox + age + rm, Hedonic, family = gaussian,
model = "random", print.level = 3, method = "nr", index = "townid")
## some count data models
data("PatentsRDUS", package="pglm")
la <- pglm(patents ~ lag(log(rd), 0:5) + scisect + log(capital72) + factor(year), PatentsRDUS,
family = negbin, model = "within", print.level = 3, method = "nr",
index = c('cusip', 'year'))
la <- pglm(patents ~ lag(log(rd), 0:5) + scisect + log(capital72) + factor(year), PatentsRDUS,
family = poisson, model = "pooling", index = c("cusip", "year"),
print.level = 0, method="nr")
## a tobit example
data("HealthIns", package="pglm")
HealthIns$med2 <- HealthIns$med / 1000
HealthIns2 <- HealthIns[-2209, ]
set.seed(2)
subs <- sample(1:20186, 200, replace = FALSE)
HealthIns2 <- HealthIns2[subs, ]
la <- pglm(med ~ mdu + disease + age, HealthIns2,
model = 'random', family = 'tobit', print.level = 0,
method = 'nr', R = 5)

Stata

Stata separates out the commands by family/link combination. So there are xtprobit and xtlogit for binary. xtoprobit for ordered probit models and xtologit [with random effects] for ordered outcomes.

Events

A Stata do file on the box is replicated here. First, a use of regular old glm.

library(haven)
Patents <- read_dta(url("https://github.com/robertwwalker/Essex-Data/raw/main/Patents.dta"))
Mod.1 <- glm(PAT~LOGR+LOGR1+LOGR2+LOGR3+LOGR4+LOGR5+dyear2+dyear3+dyear4+dyear5, data=Patents, family="poisson")
summary(Mod.1)

Call:
glm(formula = PAT ~ LOGR + LOGR1 + LOGR2 + LOGR3 + LOGR4 + LOGR5 + 
    dyear2 + dyear3 + dyear4 + dyear5, family = "poisson", data = Patents)

Deviance Residuals: 
     Min        1Q    Median        3Q       Max  
-24.1602   -3.0131   -1.6689    0.2972   23.3629  

Coefficients:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept)  1.82831    0.01253 145.861  < 2e-16 ***
LOGR         0.15251    0.02939   5.189 2.11e-07 ***
LOGR1        0.02198    0.04083   0.538 0.590297    
LOGR2        0.04369    0.03802   1.149 0.250557    
LOGR3        0.08272    0.03568   2.318 0.020427 *  
LOGR4        0.10397    0.03176   3.274 0.001061 ** 
LOGR5        0.30109    0.02063  14.594  < 2e-16 ***
dyear2      -0.04399    0.01320  -3.332 0.000862 ***
dyear3      -0.06043    0.01343  -4.499 6.81e-06 ***
dyear4      -0.18925    0.01368 -13.830  < 2e-16 ***
dyear5      -0.22979    0.01363 -16.854  < 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: 144500  on 1729  degrees of freedom
Residual deviance:  33745  on 1719  degrees of freedom
AIC: 39807

Number of Fisher Scoring iterations: 5
library(MASS)
Mod.2 <- glm.nb(PAT~LOGR+LOGR1+LOGR2+LOGR3+LOGR4+LOGR5+dyear2+dyear3+dyear4+dyear5, data=Patents)
summary(Mod.2)

Call:
glm.nb(formula = PAT ~ LOGR + LOGR1 + LOGR2 + LOGR3 + LOGR4 + 
    LOGR5 + dyear2 + dyear3 + dyear4 + dyear5, data = Patents, 
    init.theta = 1.276776486, link = log)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-2.6541  -1.0412  -0.3894   0.2827   4.3518  

Coefficients:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept)  1.19343    0.05851  20.396  < 2e-16 ***
LOGR         0.40148    0.11591   3.464 0.000533 ***
LOGR1       -0.11198    0.15672  -0.714 0.474920    
LOGR2        0.11039    0.14742   0.749 0.453957    
LOGR3        0.08770    0.13696   0.640 0.521938    
LOGR4        0.24389    0.12259   1.989 0.046657 *  
LOGR5        0.17263    0.08478   2.036 0.041729 *  
dyear2      -0.05346    0.07726  -0.692 0.488943    
dyear3      -0.05482    0.07790  -0.704 0.481611    
dyear4      -0.11764    0.07768  -1.514 0.129941    
dyear5      -0.22006    0.07764  -2.835 0.004589 ** 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

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

    Null deviance: 7202.0  on 1729  degrees of freedom
Residual deviance: 1967.8  on 1719  degrees of freedom
AIC: 11594

Number of Fisher Scoring iterations: 1

              Theta:  1.2768 
          Std. Err.:  0.0558 

 2 x log-likelihood:  -11570.1870 
library(pglm)
Loading required package: maxLik
Loading required package: miscTools

Please cite the 'maxLik' package as:
Henningsen, Arne and Toomet, Ott (2011). maxLik: A package for maximum likelihood estimation in R. Computational Statistics 26(3), 443-458. DOI 10.1007/s00180-010-0217-1.

If you have questions, suggestions, or comments regarding the 'maxLik' package, please use a forum or 'tracker' at maxLik's R-Forge site:
https://r-forge.r-project.org/projects/maxlik/
Loading required package: plm
Patents.pdata <- pdata.frame(Patents, index=c("id","YEAR"))
Mod.3A <- pglm(PAT~LOGR+LOGR1+LOGR2+LOGR3+LOGR4+LOGR5+dyear2+dyear3+dyear4+dyear5+LOGK+SCISECT, Patents.pdata, effect="individual", model="random", family="poisson")
summary(Mod.3A)
--------------------------------------------
Maximum Likelihood estimation
Newton-Raphson maximisation, 4 iterations
Return code 8: successive function values within relative tolerance limit (reltol)
Log-Likelihood: -5234.927 
14  free parameters
Estimates:
            Estimate Std. error t value  Pr(> t)    
(Intercept)  0.41079    0.14674   2.799 0.005121 ** 
LOGR         0.40345    0.04350   9.274  < 2e-16 ***
LOGR1       -0.04618    0.04822  -0.958 0.338277    
LOGR2        0.10792    0.04471   2.414 0.015788 *  
LOGR3        0.02977    0.04132   0.720 0.471221    
LOGR4        0.01070    0.03771   0.284 0.776679    
LOGR5        0.04061    0.03157   1.286 0.198364    
dyear2      -0.04496    0.01313  -3.425 0.000616 ***
dyear3      -0.04839    0.01340  -3.610 0.000306 ***
dyear4      -0.17416    0.01397 -12.467  < 2e-16 ***
dyear5      -0.22590    0.01466 -15.404  < 2e-16 ***
LOGK         0.29169    0.03934   7.415 1.21e-13 ***
SCISECT      0.25700    0.11227   2.289 0.022074 *  
sigma        1.16969    0.09471  12.350  < 2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
--------------------------------------------
Mod.3B <- pglm(PAT~LOGR+LOGR1+LOGR2+LOGR3+LOGR4+LOGR5+dyear2+dyear3+dyear4+dyear5, Patents.pdata, effect="individual", model="random", family="poisson")
summary(Mod.3B)
--------------------------------------------
Maximum Likelihood estimation
Newton-Raphson maximisation, 5 iterations
Return code 1: gradient close to zero (gradtol)
Log-Likelihood: -5263.611 
12  free parameters
Estimates:
             Estimate Std. error t value  Pr(> t)    
(Intercept)  1.402861   0.067052  20.922  < 2e-16 ***
LOGR         0.476575   0.042263  11.276  < 2e-16 ***
LOGR1       -0.007708   0.047934  -0.161 0.872242    
LOGR2        0.136413   0.044734   3.049 0.002293 ** 
LOGR3        0.059191   0.041279   1.434 0.151598    
LOGR4        0.027518   0.037601   0.732 0.464267    
LOGR5        0.082547   0.030977   2.665 0.007704 ** 
dyear2      -0.046879   0.013131  -3.570 0.000357 ***
dyear3      -0.056088   0.013362  -4.198  2.7e-05 ***
dyear4      -0.190312   0.013789 -13.802  < 2e-16 ***
dyear5      -0.252677   0.014197 -17.798  < 2e-16 ***
sigma        1.154854   0.094192  12.261  < 2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
--------------------------------------------
Mod.3C <- pglm(PAT~LOGR+LOGR1+LOGR2+LOGR3+LOGR4+LOGR5+dyear2+dyear3+dyear4+dyear5+LOGK+SCISECT-1, Patents.pdata, effect="individual", model="random", family="poisson")
summary(Mod.3C)
--------------------------------------------
Maximum Likelihood estimation
Newton-Raphson maximisation, 4 iterations
Return code 1: gradient close to zero (gradtol)
Log-Likelihood: -5238.845 
13  free parameters
Estimates:
         Estimate Std. error t value  Pr(> t)    
LOGR     0.382042   0.042877   8.910  < 2e-16 ***
LOGR1   -0.055620   0.048157  -1.155 0.248098    
LOGR2    0.100308   0.044621   2.248 0.024576 *  
LOGR3    0.020949   0.041203   0.508 0.611146    
LOGR4    0.005334   0.037709   0.141 0.887523    
LOGR5    0.028364   0.031336   0.905 0.365383    
dyear2  -0.042732   0.013109  -3.260 0.001115 ** 
dyear3  -0.044442   0.013332  -3.333 0.000858 ***
dyear4  -0.167871   0.013796 -12.168  < 2e-16 ***
dyear5  -0.216749   0.014303 -15.154  < 2e-16 ***
LOGK     0.388567   0.019946  19.481  < 2e-16 ***
SCISECT  0.419111   0.097521   4.298 1.73e-05 ***
sigma    1.117773   0.088520  12.627  < 2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
--------------------------------------------
Mod.4A <- pglm(PAT~LOGR+LOGR1+LOGR2+LOGR3+LOGR4+LOGR5+dyear2+dyear3+dyear4+dyear5+LOGK+SCISECT, Patents.pdata, effect="individual", model="random", family="negbin")
summary(Mod.4A)
--------------------------------------------
Maximum Likelihood estimation
Newton-Raphson maximisation, 6 iterations
Return code 1: gradient close to zero (gradtol)
Log-Likelihood: -4948.494 
15  free parameters
Estimates:
             Estimate Std. error t value  Pr(> t)    
(Intercept)  0.899562   0.168111   5.351 8.75e-08 ***
LOGR         0.350312   0.065282   5.366 8.04e-08 ***
LOGR1       -0.003032   0.075092  -0.040 0.967796    
LOGR2        0.104988   0.068849   1.525 0.127284    
LOGR3        0.016352   0.063638   0.257 0.797210    
LOGR4        0.035942   0.058716   0.612 0.540445    
LOGR5        0.071832   0.048289   1.488 0.136867    
dyear2      -0.043674   0.021343  -2.046 0.040734 *  
dyear3      -0.055660   0.021857  -2.547 0.010881 *  
dyear4      -0.183105   0.022718  -8.060 7.64e-16 ***
dyear5      -0.230044   0.023152  -9.936  < 2e-16 ***
LOGK         0.161937   0.041787   3.875 0.000107 ***
SCISECT      0.117642   0.106616   1.103 0.269848    
a            2.685210   0.258163  10.401  < 2e-16 ***
b            2.015688   0.217631   9.262  < 2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
--------------------------------------------
Mod.4B <- pglm(PAT~LOGR+LOGR1+LOGR2+LOGR3+LOGR4+LOGR5+dyear2+dyear3+dyear4+dyear5, Patents.pdata, effect="individual", model="random", family="negbin")
summary(Mod.4B)
--------------------------------------------
Maximum Likelihood estimation
Newton-Raphson maximisation, 6 iterations
Return code 1: gradient close to zero (gradtol)
Log-Likelihood: -4956.259 
13  free parameters
Estimates:
            Estimate Std. error t value  Pr(> t)    
(Intercept)  1.39944    0.10060  13.910  < 2e-16 ***
LOGR         0.38762    0.06492   5.970 2.37e-09 ***
LOGR1        0.01833    0.07597   0.241  0.80935    
LOGR2        0.12909    0.06972   1.852  0.06408 .  
LOGR3        0.02469    0.06446   0.383  0.70172    
LOGR4        0.05775    0.05873   0.983  0.32544    
LOGR5        0.10003    0.04750   2.106  0.03521 *  
dyear2      -0.04671    0.02180  -2.142  0.03219 *  
dyear3      -0.06512    0.02223  -2.930  0.00339 ** 
dyear4      -0.20006    0.02282  -8.765  < 2e-16 ***
dyear5      -0.25096    0.02298 -10.921  < 2e-16 ***
a            2.64396    0.25190  10.496  < 2e-16 ***
b            2.02183    0.21889   9.237  < 2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
--------------------------------------------
Mod.4C <- pglm(PAT~LOGR+LOGR1+LOGR2+LOGR3+LOGR4+LOGR5+dyear2+dyear3+dyear4+dyear5+LOGK+SCISECT-1, Patents.pdata, effect="individual", model="random", family="negbin")
summary(Mod.4C)
--------------------------------------------
Maximum Likelihood estimation
Newton-Raphson maximisation, 6 iterations
Return code 1: gradient close to zero (gradtol)
Log-Likelihood: -4963.305 
14  free parameters
Estimates:
         Estimate Std. error t value  Pr(> t)    
LOGR     0.296803   0.066943   4.434 9.26e-06 ***
LOGR1   -0.005170   0.078300  -0.066   0.9474    
LOGR2    0.085684   0.071482   1.199   0.2307    
LOGR3    0.004177   0.065959   0.063   0.9495    
LOGR4    0.010518   0.060898   0.173   0.8629    
LOGR5    0.048459   0.050342   0.963   0.3357    
dyear2  -0.032498   0.021569  -1.507   0.1319    
dyear3  -0.037854   0.021894  -1.729   0.0838 .  
dyear4  -0.159944   0.022586  -7.082 1.42e-12 ***
dyear5  -0.201171   0.022791  -8.827  < 2e-16 ***
LOGK     0.343399   0.027143  12.652  < 2e-16 ***
SCISECT  0.427951   0.089606   4.776 1.79e-06 ***
a        2.435387   0.225784  10.786  < 2e-16 ***
b        2.154018   0.240574   8.954  < 2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
--------------------------------------------

The GEE

library(geepack)
gee.Mod <- geeglm(PAT~LOGR+LOGR1+LOGR2+LOGR3+LOGR4+LOGR5+dyear2+dyear3+dyear4+dyear5, data=Patents, id=id, corstr="exchangeable", family=poisson)
summary(gee.Mod)

Call:
geeglm(formula = PAT ~ LOGR + LOGR1 + LOGR2 + LOGR3 + LOGR4 + 
    LOGR5 + dyear2 + dyear3 + dyear4 + dyear5, family = poisson, 
    data = Patents, id = id, corstr = "exchangeable")

 Coefficients:
             Estimate   Std.err    Wald Pr(>|W|)    
(Intercept)  1.851156  0.122963 226.638  < 2e-16 ***
LOGR         0.392713  0.064716  36.824 1.29e-09 ***
LOGR1        0.001683  0.057782   0.001  0.97677    
LOGR2        0.132732  0.048185   7.588  0.00588 ** 
LOGR3        0.039535  0.064560   0.375  0.54029    
LOGR4        0.036056  0.049440   0.532  0.46582    
LOGR5        0.081190  0.047667   2.901  0.08851 .  
dyear2      -0.048400  0.017823   7.375  0.00661 ** 
dyear3      -0.056190  0.027314   4.232  0.03967 *  
dyear4      -0.191152  0.042825  19.923 8.06e-06 ***
dyear5      -0.248547  0.041809  35.341 2.77e-09 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Correlation structure = exchangeable 
Estimated Scale Parameters:

            Estimate Std.err
(Intercept)    22.57   4.102
  Link = identity 

Estimated Correlation Parameters:
      Estimate Std.err
alpha   0.8983 0.06245
Number of clusters:   346  Maximum cluster size: 5 

gee.Mod <- geeglm(PAT~LOGR+LOGR1+LOGR2+LOGR3+LOGR4+LOGR5+dyear2+dyear3+dyear4+dyear5, data=Patents, id=id, corstr="ar1", family=poisson)
summary(gee.Mod)

Call:
geeglm(formula = PAT ~ LOGR + LOGR1 + LOGR2 + LOGR3 + LOGR4 + 
    LOGR5 + dyear2 + dyear3 + dyear4 + dyear5, family = poisson, 
    data = Patents, id = id, corstr = "ar1")

 Coefficients:
            Estimate Std.err   Wald Pr(>|W|)    
(Intercept)   1.8121  0.1270 203.72  < 2e-16 ***
LOGR          0.3096  0.0554  31.17  2.4e-08 ***
LOGR1         0.0254  0.0549   0.21   0.6435    
LOGR2         0.1624  0.0518   9.82   0.0017 ** 
LOGR3         0.0801  0.0564   2.02   0.1555    
LOGR4         0.0793  0.0446   3.16   0.0756 .  
LOGR5         0.0408  0.0411   0.99   0.3200    
dyear2       -0.0494  0.0175   8.01   0.0047 ** 
dyear3       -0.0542  0.0266   4.14   0.0419 *  
dyear4       -0.1817  0.0427  18.08  2.1e-05 ***
dyear5       -0.2333  0.0409  32.48  1.2e-08 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Correlation structure = ar1 
Estimated Scale Parameters:

            Estimate Std.err
(Intercept)     22.3    4.02
  Link = identity 

Estimated Correlation Parameters:
      Estimate Std.err
alpha    0.948  0.0265
Number of clusters:   346  Maximum cluster size: 5 

And an anova…

anova(gee.Mod)
Analysis of 'Wald statistic' Table
Model: poisson, link: log
Response: PAT
Terms added sequentially (first to last)

       Df  X2 P(>|Chi|)    
LOGR    1 368   < 2e-16 ***
LOGR1   1  23   1.9e-06 ***
LOGR2   1  19   1.4e-05 ***
LOGR3   1  16   5.8e-05 ***
LOGR4   1   7   0.00708 ** 
LOGR5   1   1   0.28590    
dyear2  1   4   0.05740 .  
dyear3  1  11   0.00099 ***
dyear4  1   0   0.67359    
dyear5  1  32   1.2e-08 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

BKT and Carter-Signorino

Load the data

BKT.data <- read_dta(url("https://github.com/robertwwalker/Essex-Data/raw/main/bkt98ajps.dta"))
Mod.1 <- glm(dispute~ dem+growth+allies+contig+capratio+trade, family=binomial("logit"), data=BKT.data)
summary(Mod.1)

Call:
glm(formula = dispute ~ dem + growth + allies + contig + capratio + 
    trade, family = binomial("logit"), data = BKT.data)

Deviance Residuals: 
   Min      1Q  Median      3Q     Max  
-0.748  -0.319  -0.261  -0.152   4.096  

Coefficients:
             Estimate Std. Error z value Pr(>|z|)    
(Intercept) -3.29e+00   7.92e-02  -41.54  < 2e-16 ***
dem         -4.97e-02   7.38e-03   -6.73  1.7e-11 ***
growth      -2.23e-02   8.52e-03   -2.62   0.0088 ** 
allies      -8.21e-01   8.00e-02  -10.26  < 2e-16 ***
contig       1.31e+00   7.96e-02   16.49  < 2e-16 ***
capratio    -3.07e-03   4.17e-04   -7.37  1.8e-13 ***
trade       -6.61e+01   1.34e+01   -4.92  8.6e-07 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 7719.2  on 20989  degrees of freedom
Residual deviance: 6955.1  on 20983  degrees of freedom
AIC: 6969

Number of Fisher Scoring iterations: 8
Mod.1B <- glm(dispute~ dem+growth+allies+contig+capratio+trade+as.factor(py), family=binomial("logit"), data=BKT.data)
Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
summary(Mod.1B)

Call:
glm(formula = dispute ~ dem + growth + allies + contig + capratio + 
    trade + as.factor(py), family = binomial("logit"), data = BKT.data)

Deviance Residuals: 
   Min      1Q  Median      3Q     Max  
-1.362  -0.225  -0.140  -0.080   3.988  

Coefficients:
                 Estimate Std. Error z value Pr(>|z|)    
(Intercept)     -9.43e-01   9.28e-02  -10.16  < 2e-16 ***
dem             -5.47e-02   8.00e-03   -6.84  8.0e-12 ***
growth          -1.15e-02   9.24e-03   -1.24  0.21492    
allies          -4.71e-01   9.00e-02   -5.23  1.7e-07 ***
contig           6.99e-01   8.93e-02    7.82  5.3e-15 ***
capratio        -3.03e-03   4.17e-04   -7.28  3.3e-13 ***
trade           -1.27e+01   1.05e+01   -1.21  0.22735    
as.factor(py)2  -2.04e+00   1.41e-01  -14.48  < 2e-16 ***
as.factor(py)3  -2.15e+00   1.54e-01  -13.92  < 2e-16 ***
as.factor(py)4  -2.37e+00   1.74e-01  -13.64  < 2e-16 ***
as.factor(py)5  -2.62e+00   2.01e-01  -13.04  < 2e-16 ***
as.factor(py)6  -3.09e+00   2.52e-01  -12.25  < 2e-16 ***
as.factor(py)7  -3.08e+00   2.52e-01  -12.21  < 2e-16 ***
as.factor(py)8  -3.00e+00   2.46e-01  -12.23  < 2e-16 ***
as.factor(py)9  -2.95e+00   2.46e-01  -11.98  < 2e-16 ***
as.factor(py)10 -3.24e+00   2.86e-01  -11.32  < 2e-16 ***
as.factor(py)11 -3.45e+00   3.24e-01  -10.66  < 2e-16 ***
as.factor(py)12 -3.32e+00   3.10e-01  -10.72  < 2e-16 ***
as.factor(py)13 -3.92e+00   4.14e-01   -9.46  < 2e-16 ***
as.factor(py)14 -4.56e+00   5.82e-01   -7.84  4.6e-15 ***
as.factor(py)15 -4.58e+00   5.82e-01   -7.87  3.5e-15 ***
as.factor(py)16 -3.45e+00   3.41e-01  -10.12  < 2e-16 ***
as.factor(py)17 -4.22e+00   5.05e-01   -8.35  < 2e-16 ***
as.factor(py)18 -5.59e+00   1.00e+00   -5.57  2.5e-08 ***
as.factor(py)19 -3.61e+00   3.85e-01   -9.36  < 2e-16 ***
as.factor(py)20 -3.95e+00   4.53e-01   -8.71  < 2e-16 ***
as.factor(py)21 -4.80e+00   7.11e-01   -6.74  1.5e-11 ***
as.factor(py)22 -4.31e+00   5.82e-01   -7.39  1.4e-13 ***
as.factor(py)23 -4.13e+00   5.83e-01   -7.08  1.4e-12 ***
as.factor(py)24 -3.99e+00   5.83e-01   -6.84  7.7e-12 ***
as.factor(py)25 -1.76e+01   3.31e+02   -0.05  0.95752    
as.factor(py)26 -2.97e+00   4.19e-01   -7.10  1.2e-12 ***
as.factor(py)27 -1.76e+01   3.85e+02   -0.05  0.96353    
as.factor(py)28 -1.76e+01   3.83e+02   -0.05  0.96347    
as.factor(py)29 -3.92e+00   7.14e-01   -5.49  3.9e-08 ***
as.factor(py)30 -3.42e+00   5.86e-01   -5.84  5.2e-09 ***
as.factor(py)31 -4.41e+00   1.01e+00   -4.39  1.1e-05 ***
as.factor(py)32 -3.57e+00   7.15e-01   -4.99  6.0e-07 ***
as.factor(py)33 -4.07e+00   1.01e+00   -4.05  5.2e-05 ***
as.factor(py)34 -3.31e+00   7.16e-01   -4.61  3.9e-06 ***
as.factor(py)35 -3.83e+00   1.01e+00   -3.81  0.00014 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 7719.2  on 20989  degrees of freedom
Residual deviance: 5109.4  on 20949  degrees of freedom
AIC: 5191

Number of Fisher Scoring iterations: 17
Mod.1C <- glm(dispute~ dem+growth+allies+contig+capratio+trade+py+pys1+pys2+pys3, family=binomial("logit"), data=BKT.data)
summary(Mod.1C)

Call:
glm(formula = dispute ~ dem + growth + allies + contig + capratio + 
    trade + py + pys1 + pys2 + pys3, family = binomial("logit"), 
    data = BKT.data)

Deviance Residuals: 
   Min      1Q  Median      3Q     Max  
-1.349  -0.221  -0.138  -0.092   4.063  

Coefficients:
             Estimate Std. Error z value Pr(>|z|)    
(Intercept)  8.50e-01   1.60e-01    5.31  1.1e-07 ***
dem         -5.46e-02   7.97e-03   -6.84  7.7e-12 ***
growth      -1.15e-02   9.20e-03   -1.25     0.21    
allies      -4.70e-01   8.96e-02   -5.25  1.5e-07 ***
contig       6.94e-01   8.90e-02    7.79  6.6e-15 ***
capratio    -3.04e-03   4.17e-04   -7.29  3.0e-13 ***
trade       -1.29e+01   1.05e+01   -1.23     0.22    
py          -1.82e+00   1.11e-01  -16.30  < 2e-16 ***
pys1        -2.45e-01   2.61e-02   -9.37  < 2e-16 ***
pys2         7.92e-02   1.09e-02    7.24  4.5e-13 ***
pys3        -1.09e-02   2.76e-03   -3.96  7.5e-05 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 7719.2  on 20989  degrees of freedom
Residual deviance: 5165.8  on 20979  degrees of freedom
AIC: 5188

Number of Fisher Scoring iterations: 8
Mod.2A <- glm(dispute~ dem+growth+allies+contig+capratio+trade, family=binomial("logit"), subset = contdisp!=1, data=BKT.data)
summary(Mod.2A)

Call:
glm(formula = dispute ~ dem + growth + allies + contig + capratio + 
    trade, family = binomial("logit"), data = BKT.data, subset = contdisp != 
    1)

Deviance Residuals: 
   Min      1Q  Median      3Q     Max  
-0.520  -0.218  -0.169  -0.120   3.868  

Coefficients:
             Estimate Std. Error z value Pr(>|z|)    
(Intercept) -4.33e+00   1.15e-01  -37.78  < 2e-16 ***
dem         -4.01e-02   1.01e-02   -3.99  6.7e-05 ***
growth      -3.43e-02   1.25e-02   -2.74   0.0062 ** 
allies      -4.80e-01   1.13e-01   -4.25  2.1e-05 ***
contig       1.35e+00   1.21e-01   11.20  < 2e-16 ***
capratio    -1.96e-03   5.01e-04   -3.92  9.0e-05 ***
trade       -2.11e+01   1.13e+01   -1.86   0.0623 .  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 3978.5  on 20447  degrees of freedom
Residual deviance: 3693.8  on 20441  degrees of freedom
AIC: 3708

Number of Fisher Scoring iterations: 9
Mod.2B <- glm(dispute~ dem+growth+allies+contig+capratio+trade+py+pys1+pys2+pys3, family=binomial("logit"), subset = contdisp!=1, data=BKT.data)
summary(Mod.2B)

Call:
glm(formula = dispute ~ dem + growth + allies + contig + capratio + 
    trade + py + pys1 + pys2 + pys3, family = binomial("logit"), 
    data = BKT.data, subset = contdisp != 1)

Deviance Residuals: 
   Min      1Q  Median      3Q     Max  
-0.734  -0.210  -0.135  -0.098   3.732  

Coefficients:
             Estimate Std. Error z value Pr(>|z|)    
(Intercept) -3.956674   0.297828  -13.29  < 2e-16 ***
dem         -0.038809   0.010046   -3.86  0.00011 ***
growth      -0.040144   0.012502   -3.21  0.00132 ** 
allies      -0.365957   0.114272   -3.20  0.00136 ** 
contig       0.993251   0.124133    8.00  1.2e-15 ***
capratio    -0.002289   0.000523   -4.38  1.2e-05 ***
trade       -3.809220   9.679795   -0.39  0.69393    
py           0.390965   0.157920    2.48  0.01330 *  
pys1         0.089719   0.031787    2.82  0.00477 ** 
pys2        -0.028888   0.012634   -2.29  0.02222 *  
pys3         0.003314   0.002974    1.11  0.26513    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 3978.5  on 20447  degrees of freedom
Residual deviance: 3502.8  on 20437  degrees of freedom
AIC: 3525

Number of Fisher Scoring iterations: 9
Mod.3A <- glm(dispute~ dem+growth+allies+contig+capratio+trade+py+pys1+pys2+pys3, family=binomial("logit"), subset = prefail<1, data=BKT.data)
summary(Mod.3A)

Call:
glm(formula = dispute ~ dem + growth + allies + contig + capratio + 
    trade + py + pys1 + pys2 + pys3, family = binomial("logit"), 
    data = BKT.data, subset = prefail < 1)

Deviance Residuals: 
   Min      1Q  Median      3Q     Max  
-0.682  -0.159  -0.109  -0.081   3.856  

Coefficients:
             Estimate Std. Error z value Pr(>|z|)    
(Intercept) -2.132791   0.369090   -5.78  7.5e-09 ***
dem         -0.046143   0.013330   -3.46  0.00054 ***
growth      -0.022954   0.017754   -1.29  0.19604    
allies      -0.423113   0.163691   -2.58  0.00974 ** 
contig       1.105361   0.172537    6.41  1.5e-10 ***
capratio    -0.001916   0.000601   -3.19  0.00142 ** 
trade       -3.553603  11.726412   -0.30  0.76186    
py          -1.081689   0.235688   -4.59  4.4e-06 ***
pys1        -0.183730   0.050096   -3.67  0.00024 ***
pys2         0.069295   0.019788    3.50  0.00046 ***
pys3        -0.013987   0.004484   -3.12  0.00181 ** 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 2183.3  on 16990  degrees of freedom
Residual deviance: 1928.1  on 16980  degrees of freedom
AIC: 1950

Number of Fisher Scoring iterations: 9

Carter and Signorino

Use a Taylor series to clean this up with a suggestion that no baseline hazard is more than cubic [two inflection points].

Mod.CS <- glm(dispute~ dem+growth+allies+contig+capratio+trade+prefail+py+I(py^2)+I(py^3), family=binomial("logit"), data=BKT.data)
summary(Mod.CS)

Call:
glm(formula = dispute ~ dem + growth + allies + contig + capratio + 
    trade + prefail + py + I(py^2) + I(py^3), family = binomial("logit"), 
    data = BKT.data)

Deviance Residuals: 
   Min      1Q  Median      3Q     Max  
-2.576  -0.192  -0.128  -0.092   3.754  

Coefficients:
             Estimate Std. Error z value Pr(>|z|)    
(Intercept) -1.14e+00   1.18e-01   -9.64  < 2e-16 ***
dem         -4.09e-02   8.23e-03   -4.97  6.7e-07 ***
growth      -2.44e-02   9.55e-03   -2.56   0.0105 *  
allies      -2.67e-01   9.20e-02   -2.90   0.0038 ** 
contig       6.62e-01   9.20e-02    7.19  6.4e-13 ***
capratio    -2.00e-03   3.71e-04   -5.40  6.6e-08 ***
trade       -1.06e+01   1.04e+01   -1.03   0.3049    
prefail      1.75e-01   8.96e-03   19.53  < 2e-16 ***
py          -7.48e-01   4.14e-02  -18.08  < 2e-16 ***
I(py^2)      4.22e-02   3.72e-03   11.34  < 2e-16 ***
I(py^3)     -7.18e-04   8.76e-05   -8.19  2.5e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 7719.2  on 20989  degrees of freedom
Residual deviance: 4909.5  on 20979  degrees of freedom
AIC: 4932

Number of Fisher Scoring iterations: 8

As a technical matter, none of this is exactly equivalent to a Cox model. That requires the cloglog.

Mod.CS <- glm(dispute~ dem+growth+allies+contig+capratio+trade+prefail+py+I(py^2)+I(py^3), family=binomial("cloglog"), data=BKT.data)
summary(Mod.CS)

Call:
glm(formula = dispute ~ dem + growth + allies + contig + capratio + 
    trade + prefail + py + I(py^2) + I(py^3), family = binomial("cloglog"), 
    data = BKT.data)

Deviance Residuals: 
   Min      1Q  Median      3Q     Max  
-3.396  -0.194  -0.134  -0.097   3.700  

Coefficients:
             Estimate Std. Error z value Pr(>|z|)    
(Intercept) -1.13e+00   1.07e-01  -10.62  < 2e-16 ***
dem         -4.16e-02   7.50e-03   -5.54  2.9e-08 ***
growth      -1.97e-02   8.00e-03   -2.47   0.0137 *  
allies      -2.33e-01   7.99e-02   -2.92   0.0035 ** 
contig       5.16e-01   7.86e-02    6.56  5.4e-11 ***
capratio    -2.13e-03   3.60e-04   -5.93  3.0e-09 ***
trade       -1.15e+01   1.01e+01   -1.13   0.2567    
prefail      1.15e-01   6.46e-03   17.80  < 2e-16 ***
py          -7.03e-01   3.84e-02  -18.32  < 2e-16 ***
I(py^2)      3.89e-02   3.51e-03   11.07  < 2e-16 ***
I(py^3)     -6.55e-04   8.33e-05   -7.85  4.1e-15 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 7719.2  on 20989  degrees of freedom
Residual deviance: 4944.0  on 20979  degrees of freedom
AIC: 4966

Number of Fisher Scoring iterations: 9