postscript("l13.ps");options(width=68)
#*************************************************************/
# Mark A: Exact logistic regression.                         */
#*************************************************************/
#*************************************************************/
# Data obtained from http://www-fars.nhtsa.dot.gov by re-    */
# questing violations, accident severity, body type, speed   */
# limit, and restraint use for NJ, 1999.  Variables are id,  */
# id=line in the data set                                    */
# state=34 for New Jersey                                    */
# case=a numberof the accident                               */
# vh=number of the vehicle within accident                   */
# pn=number of passenger within vehicle                      */
# sl=speed limit where the accident occurred                 */
# bt=body type of vehicle                                    */
# v1-v3=codes for violations cited.                          */
# rs=whether proper restraints (ex. seat belts) were used    */
# is=injury severity                                         */
# Each accident potentially involves multiple vehicles, and  */
# each vehicle potentially involves multiple passengers.     */
# Select one vehicle (vh=1) from each crash.  Select one     */
# passenger (pn=1) from  this vehicle.                       */
#*************************************************************/
pcrash<-read.table("crash.dat",col.names=c("id","state","case",
  "vh","pn","sl","bt","v1","v2","v3","rs","is"),as.is=TRUE)
pcrash$rs<-as.numeric(pcrash$rs)
pcrash$fatal<-c("Not Fatal","Fatal")[(pcrash$is==4)+1]
pcrash$vt<-c("Car","Light truck","Van","Truck")[
  1+(pcrash$bt>10)+(pcrash$bt>20)+(pcrash$bt>30)]
pcrash$beltuse<-paste(c("Proper","Improper or no"), 
  "restraint use")[1+((pcrash$rs>12)|(pcrash$rs==0))]
pcrash<-pcrash[((pcrash$vh==1)|(pcrash$vh==2))&(pcrash$pn==1),
  c("case","vh","vt","fatal","beltuse")]
crash<-pcrash[pcrash$vh==1,c("vt","fatal","beltuse")]
#*************************************************************/
# Examine the data.  Note single-vehicle starata.            */
#*************************************************************/
#*************************************************************/
# A stratified approach is automatically conditional, but not*/
# exact.  This technique cannot give an estimate for Speed   */
# limit, since speeed limit is constant over each stratum.   */
#*************************************************************/
cat('\n Stratified results with 1 entry strata left in  \n')
library(survival)
pcrash$fatall<-pcrash$fatal!="Fatal"
clogit(fatall~beltuse+strata(case),data=pcrash)
#*************************************************************/
# Remove single-vehicle strata, and note identical results.  */
#*************************************************************/
cat('\n Stratified results with 1 entry strata removed   \n')
ok<-apply(table(pcrash$fatal,pcrash$case),2,prod)==1
discordantcases<-as.numeric(names(ok)[ok])
clogit(fatall~beltuse+strata(case),data=pcrash,
  subset=(pcrash$case%in%discordantcases))
#*************************************************************/
# The above analysis gave normal theory results.  Force exact*/
# results.                                                   */
#*************************************************************/
cat('\n Stratified exact results  \n')
library(elrm)
pcrash$n<-1
#Quite slow in R.  I didn't wait for it to finish.
#elrm(fatall/n~beltuse+strata(case),~beltuse,dataset=pcrash)
#*************************************************************/
# Stratification as above indicates to SAS that it needs to  */
# remove the intercept via conditioning.  SAS can use the in-*/
# formation that data are in strata to more efficiently do   */
# the exact algorithm.  The exact algorithm could also be ap-*/
# lied directly via a class variable, but this is slow.  Om- */
# mission of strata forces construction of agreement table   */
# including pairs, which is not what we need.  Estimates are */
# identical to asymptotic stratified approach, and exact     */
# stratified approach.  Tests are the same as the exact      */
# stratified approach.                                       */
#*************************************************************/
cat('\n General regression exact results  \n')
#Quite slow in R.  I didn't wait for it to finish.
#elrm(fatall/n~beltuse+strata(case),~beltuse,dataset=pcrash)