R version 4.3.1 (2023-06-16) -- "Beagle Scouts"
Copyright (C) 2023 The R Foundation for Statistical Computing
Platform: x86_64-pc-linux-gnu (64-bit)

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.

  Natural language support but running in an English locale

R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.

Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.

> pdf("g06.pdf");options(width=64)
> #setwd("C:\\Users\\kolassa\\Class553")
> setwd("~/Taught1/960-553/Data")
> #**************************************************/
> #* Multi-Way Analyses                             */
> #**************************************************/
> #*************************************************************/
> #*  Prostate data                                            */
> #* From Andrews and Herzberg (1985), Data: A Collection of   */
> #* Problems from Many Fields for the Student and Research    */
> #* Worker, Table 46.  Observations represent subjects in a   */
> #* prostate cohort study, randomized to one of four dose     */
> #* levels of diethylstilbestrol.  Rx records dose in four    */
> #* ordered categories, with 1 being placebo.  Disease stage  */
> #* is 3 or 4.  monfol is months of followup.  surv is 0 if   */
> #* alive after 50 mo, and codes cause of death otherwise.    */
> #* http://lib.stat.cmu.edu/datasets/Andrews/T46.1            */
> #* The on-line version of the data set adds 3 fields before  */
> #* the first field in the book.  Variables of interest are   */
> #* stage, rx, monfol, and surv in fields 5, 6, 10, 11 of the */
> #* online version, resp.  Causes of death are given by var-  */
> #* ious positive integers in surv; I recode these to 1.  The */
> #* data file has more columns than we are interested in.  Put*/
> #* place-holding variables in for variables not of interest  */
> #* between variables of interest.  Data were previously pub- */
> #* lished by Byar and Corle (1977 Chronic Disease) and Byar  */
> #* and Green (1980 Bull. Cancer).  Lower value of dichoto-   */
> #* mized dose begins with blank to make it alphabetized      */
> #* before high.                                              */
> #*************************************************************/
> # If we omit the ones past where the data of interest stops
> # in the "what" list, and use flush=T, R will ignore them. 
> prostate<-read.table("T46.1")[,c(5,6,10,11)]
> dimnames(prostate)[[2]]<-c("stage","rx","monfol","surv")
> prostate$alive<-(prostate$surv==0)+0                      
> prostate$dose<-c(" low","high")[1+(prostate$rx>1)]         
> stage4<-prostate[prostate$stage==4,]
> #*******************************************************/
> #* Test whether survival is the same across treatments.*/
> #*******************************************************/
> # Run Poisson regression.  We need to change this into a data
> # frame with a variable representing counts.  
> # First fit additive model.
> newdata<-as.data.frame(xtabs(~rx+alive,data=stage4))
> summary(m1<-glm(Freq~rx+alive,family=poisson, data=newdata))

Call:
glm(formula = Freq ~ rx + alive, family = poisson, data = newdata)

Coefficients:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept)  3.68572    0.14285  25.801  < 2e-16 ***
rx2         -0.01905    0.19519  -0.098    0.922    
rx3          0.03704    0.19248   0.192    0.847    
rx4          0.01869    0.19336   0.097    0.923    
alive1      -1.11111    0.15836  -7.016 2.28e-12 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for poisson family taken to be 1)

    Null deviance: 63.094  on 7  degrees of freedom
Residual deviance:  5.908  on 3  degrees of freedom
AIC: 55.605

Number of Fisher Scoring iterations: 4

> # Now fit model with interactions.
> m2<-glm(Freq~rx*alive,family=poisson,data=newdata)
> # Do we need interactions?
> (anovaout<-anova(m1,m2))
Analysis of Deviance Table

Model 1: Freq ~ rx + alive
Model 2: Freq ~ rx * alive
  Resid. Df Resid. Dev Df Deviance
1         3      5.908            
2         0      0.000  3    5.908
> # Calculate p value
> 1-pchisq(anovaout$Deviance[2],anovaout$Df[2])
[1] 0.1161754
> # Compare to earlier Pearson test
> chisq.test(stage4$rx,stage4$alive)

	Pearson's Chi-squared test

data:  stage4$rx and stage4$alive
X-squared = 6.0722, df = 3, p-value = 0.1081

> #In the above call to glm, rx was treated as categorical, 
> #because, while it was originally numeric, it got changed
> #to a character string when it got used as the labels for
> #table rows, and this got carried over into the data frame
> #created from teh table.  If we want this treated as
> #continuous, convert back to numeric.
> summary(m3<-glm(Freq~as.numeric(rx)+alive,family=poisson, 
+    data=newdata))

Call:
glm(formula = Freq ~ as.numeric(rx) + alive, family = poisson, 
    data = newdata)

Coefficients:
               Estimate Std. Error z value Pr(>|z|)    
(Intercept)     3.66699    0.17275  21.228  < 2e-16 ***
as.numeric(rx)  0.01122    0.06115   0.183    0.854    
alive1         -1.11111    0.15836  -7.016 2.28e-12 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for poisson family taken to be 1)

    Null deviance: 63.0938  on 7  degrees of freedom
Residual deviance:  5.9678  on 5  degrees of freedom
AIC: 51.665

Number of Fisher Scoring iterations: 4

> #**************************************************************/
> #* Shipping Data: McCullagh and Nelder (1989) Generalized     */
> #* Linear Models provide data on losses by a shipping insurer.*/
> #* Data are grouped into putatively homogeneous classes, based*/
> #* on ship type, start of construction 5 year period, star of */
> #* observation 5 year period, ship months at risk, and count  */
> #* of losses.  Empty categories are removed.                  */
> #**************************************************************/
> ships<-read.table("ships.dat",
+    col.names=c("type","built", "period","smar","cases"))
> #You can save the result of a complicated model in R just as
> #you can save any other object.  Printing the result of a 
> #complicated model in R is generally not informative.
> glmo<-glm(cases~type+built,family=poisson,data=ships)
> cat("Result of printing glm output\n")
Result of printing glm output
> glmo

Call:  glm(formula = cases ~ type + built, family = poisson, data = ships)

Coefficients:
(Intercept)        typeB        typeC        typeD  
   1.456470     1.795720    -1.252763    -0.904456  
      typeE        built  
  -0.123139     0.005042  

Degrees of Freedom: 33 Total (i.e. Null);  28 Residual
Null Deviance:	    614.5 
Residual Deviance: 170.5 	AIC: 280.3
> #Instead, print the results of a summary.
> cat("Result of printing glm summary\n")
Result of printing glm summary
> summary(glmo)

Call:
glm(formula = cases ~ type + built, family = poisson, data = ships)

Coefficients:
             Estimate Std. Error z value Pr(>|z|)    
(Intercept)  1.456470   0.704790   2.067  0.03878 *  
typeB        1.795720   0.166620  10.777  < 2e-16 ***
typeC       -1.252763   0.327327  -3.827  0.00013 ***
typeD       -0.904456   0.287459  -3.146  0.00165 ** 
typeE       -0.123139   0.234900  -0.524  0.60013    
built        0.005042   0.010331   0.488  0.62551    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for poisson family taken to be 1)

    Null deviance: 614.54  on 33  degrees of freedom
Residual deviance: 170.47  on 28  degrees of freedom
AIC: 280.34

Number of Fisher Scoring iterations: 6

> summary(glmo1<-glm(cases~type+factor(built)+
+    offset(log(smar)),family=poisson,data=ships))

Call:
glm(formula = cases ~ type + factor(built) + offset(log(smar)), 
    family = poisson, data = ships)

Coefficients:
                Estimate Std. Error z value Pr(>|z|)    
(Intercept)     -6.26892    0.21241 -29.513  < 2e-16 ***
typeB           -0.55861    0.17767  -3.144  0.00167 ** 
typeC           -0.68388    0.32908  -2.078  0.03769 *  
typeD           -0.07032    0.29071  -0.242  0.80887    
typeE            0.30425    0.23581   1.290  0.19698    
factor(built)65  0.75347    0.14855   5.072 3.94e-07 ***
factor(built)70  0.96530    0.16369   5.897 3.70e-09 ***
factor(built)75  0.70828    0.22153   3.197  0.00139 ** 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for poisson family taken to be 1)

    Null deviance: 146.328  on 33  degrees of freedom
Residual deviance:  49.355  on 26  degrees of freedom
AIC: 163.22

Number of Fisher Scoring iterations: 5

> summary(glmo2<-glm(cases~type+factor(built)+factor(period)+
+    offset(log(smar)),family=poisson,data=ships))

Call:
glm(formula = cases ~ type + factor(built) + factor(period) + 
    offset(log(smar)), family = poisson, data = ships)

Coefficients:
                 Estimate Std. Error z value Pr(>|z|)    
(Intercept)      -6.40590    0.21744 -29.460  < 2e-16 ***
typeB            -0.54334    0.17759  -3.060  0.00222 ** 
typeC            -0.68740    0.32904  -2.089  0.03670 *  
typeD            -0.07596    0.29058  -0.261  0.79377    
typeE             0.32558    0.23588   1.380  0.16750    
factor(built)65   0.69714    0.14964   4.659 3.18e-06 ***
factor(built)70   0.81843    0.16977   4.821 1.43e-06 ***
factor(built)75   0.45343    0.23317   1.945  0.05182 .  
factor(period)75  0.38447    0.11827   3.251  0.00115 ** 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for poisson family taken to be 1)

    Null deviance: 146.328  on 33  degrees of freedom
Residual deviance:  38.695  on 25  degrees of freedom
AIC: 154.56

Number of Fisher Scoring iterations: 5

> #**************************************************/
> #* Models with Interaction                        */
> #**************************************************/
> summary(glmo3<-glm(cases~type*factor(built)+offset(log(smar)),
+    family=poisson,data=ships))

Call:
glm(formula = cases ~ type * factor(built) + offset(log(smar)), 
    family = poisson, data = ships)

Coefficients:
                        Estimate Std. Error z value Pr(>|z|)
(Intercept)             -23.7962  6467.9269  -0.004    0.997
typeB                    16.9799  6467.9269   0.003    0.998
typeC                    17.0329  6467.9269   0.003    0.998
typeD                    -0.5960  9075.0371   0.000    1.000
typeE                     0.6870 11432.1995   0.000    1.000
factor(built)65          18.0505  6467.9269   0.003    0.998
factor(built)70          18.4845  6467.9269   0.003    0.998
factor(built)75          18.4781  6467.9269   0.003    0.998
typeB:factor(built)65   -17.3238  6467.9269  -0.003    0.998
typeC:factor(built)65   -18.5713  6467.9270  -0.003    0.998
typeD:factor(built)65   -18.4210 11220.0344  -0.002    0.999
typeE:factor(built)65     0.5863 11432.1995   0.000    1.000
typeB:factor(built)70   -17.5544  6467.9269  -0.003    0.998
typeC:factor(built)70   -17.5542  6467.9270  -0.003    0.998
typeD:factor(built)70     1.1222  9075.0371   0.000    1.000
typeE:factor(built)70    -0.6491 11432.1995   0.000    1.000
typeB:factor(built)75   -17.6417  6467.9269  -0.003    0.998
typeC:factor(built)75   -17.3279  6467.9270  -0.003    0.998
typeD:factor(built)75    -0.3256  9075.0371   0.000    1.000
typeE:factor(built)75    -1.6641 11432.1995   0.000    1.000

(Dispersion parameter for poisson family taken to be 1)

    Null deviance: 146.328  on 33  degrees of freedom
Residual deviance:  25.208  on 14  degrees of freedom
AIC: 163.07

Number of Fisher Scoring iterations: 17

> #Do we need the interactions? 
> (anovaout<-anova(glmo2,glmo3))
Analysis of Deviance Table

Model 1: cases ~ type + factor(built) + factor(period) + offset(log(smar))
Model 2: cases ~ type * factor(built) + offset(log(smar))
  Resid. Df Resid. Dev Df Deviance
1        25     38.695            
2        14     25.208 11   13.487
> # Calculate p value
> 1-pchisq(anovaout$Deviance[2],anovaout$Df[2])
[1] 0.2627082
> #**************************************************************/
> #*          Death Penalty Example                             */
> #* Death penalty data from Radelet, M. (1981), Racial         */
> #* Characteristics and Imposition of the Death Penalty        */
> #* American Sociological Review, 46, 918-927. Data reflect    */
> #* outcomes of murder cases in which race of victim and       */
> #* principle suspect are known.  1st column is indicator of   */
> #* relationship between victim and suspect (P for close       */
> #* relationship, N for not), Race of victim and defendant     */
> #* (both either B or W; other races deleted), number of cases,*/
> #* number of indictments, and number of death sentences.      */
> #* Referenced by Moore (2000) Basic Practice of Statistics, p.*/
> #* Following Moore, we will only use the non-primary cases.   */
> #* http://stat.rutgers.edu/~kolassa/960-553/floridadeath.dat  */
> #**************************************************************/
> death1<-read.table("floridadeath.dat",
+    col.names=c("rela","rv","rd","nc","ni","nd"))
> death2<-death1; death1$we<-death1$nd;death1$lc<-"Death" 
> death2$we<-death1$nc-death1$nd;death2$lc<-"Live"
> death<-rbind(death1,death2)
> death$rdc<-c("White defend.","Black defend.")[(death$rd=="B")+1]
> death$rvc<-c("White victim","Black victim")[(death$rv=="B")+1]
> nonprim<-death[death$rela=="N",]                       
> # Fit model ignoring race of victim.  The * indicates 
> # interactions, which forces main effects.
> summary(d1<-glm(we~rdc*lc,data=nonprim,family=poisson))

Call:
glm(formula = we ~ rdc * lc, family = poisson, data = nonprim)

Coefficients:
                        Estimate Std. Error z value Pr(>|z|)
(Intercept)               2.1401     0.2425   8.824   <2e-16
rdcWhite defend.          0.1112     0.3338   0.333    0.739
lcLive                    2.1707     0.2560   8.479   <2e-16
rdcWhite defend.:lcLive  -0.1664     0.3539  -0.470    0.638
                           
(Intercept)             ***
rdcWhite defend.           
lcLive                  ***
rdcWhite defend.:lcLive    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for poisson family taken to be 1)

    Null deviance: 395.92  on 7  degrees of freedom
Residual deviance: 170.16  on 4  degrees of freedom
AIC: 213.85

Number of Fisher Scoring iterations: 5

> # Allow for different marginal probability for each race
> # of victim, but odds ratio must be the same.
> summary(d2<-glm(we~ rdc*lc+ rdc*rvc+ lc*rvc,data=nonprim,
+    family=poisson))

Call:
glm(formula = we ~ rdc * lc + rdc * rvc + lc * rvc, family = poisson, 
    data = nonprim)

Coefficients:
                                 Estimate Std. Error z value
(Intercept)                        1.7360     0.4092   4.242
rdcWhite defend.                  -2.8579     0.5195  -5.502
lcLive                             2.8421     0.4203   6.761
rvcWhite victim                    0.6911     0.4947   1.397
rdcWhite defend.:lcLive            0.4402     0.4009   1.098
rdcWhite defend.:rvcWhite victim   3.3580     0.3820   8.791
lcLive:rvcWhite victim            -1.3242     0.5193  -2.550
                                 Pr(>|z|)    
(Intercept)                      2.21e-05 ***
rdcWhite defend.                 3.76e-08 ***
lcLive                           1.37e-11 ***
rvcWhite victim                    0.1624    
rdcWhite defend.:lcLive            0.2722    
rdcWhite defend.:rvcWhite victim  < 2e-16 ***
lcLive:rvcWhite victim             0.0108 *  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for poisson family taken to be 1)

    Null deviance: 395.91531  on 7  degrees of freedom
Residual deviance:   0.70074  on 1  degrees of freedom
AIC: 50.382

Number of Fisher Scoring iterations: 4

> #In mantelhaen.test last entry in model gives strata, and so
> #estimated OR measures association betw rdc and lc.  Est-   
> #imator is MH estimator.                                    
> #***********************************************************/
> #*Compare OR estimate to exp of rdc*lc interaction in      */
> #*Poisson regression d2 with all 2-way interactions above. */
> #***********************************************************/
> mantelhaen.test(mm<-xtabs(we~rdc+lc+rvc,data=nonprim))

	Mantel-Haenszel chi-squared test with continuity
	correction

data:  mm <- xtabs(we ~ rdc + lc + rvc, data = nonprim)
Mantel-Haenszel X-squared = 0.79633, df = 1, p-value =
0.3722
alternative hypothesis: true common odds ratio is not equal to 1
95 percent confidence interval:
 0.7096435 3.4914549
sample estimates:
common odds ratio 
         1.574067 

> # Allow different odds ratio for relationship between race of
> # the defendant and death penalty for each race of victim.
> summary(d3<-glm(we~ rdc*lc*rvc,data=nonprim,
+    family=poisson))

Call:
glm(formula = we ~ rdc * lc * rvc, family = poisson, data = nonprim)

Coefficients:
                                          Estimate Std. Error
(Intercept)                                 1.7918     0.4082
rdcWhite defend.                          -24.0943 42247.1657
lcLive                                      2.7830     0.4207
rvcWhite victim                             0.6061     0.5075
rdcWhite defend.:lcLive                    21.7169 42247.1657
rdcWhite defend.:rvcWhite victim           24.6409 42247.1657
lcLive:rvcWhite victim                     -1.2296     0.5358
rdcWhite defend.:lcLive:rvcWhite victim   -21.3318 42247.1657
                                        z value Pr(>|z|)    
(Intercept)                               4.389 1.14e-05 ***
rdcWhite defend.                         -0.001   0.9995    
lcLive                                    6.615 3.71e-11 ***
rvcWhite victim                           1.194   0.2324    
rdcWhite defend.:lcLive                   0.001   0.9996    
rdcWhite defend.:rvcWhite victim          0.001   0.9995    
lcLive:rvcWhite victim                   -2.295   0.0217 *  
rdcWhite defend.:lcLive:rvcWhite victim  -0.001   0.9996    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for poisson family taken to be 1)

    Null deviance: 3.9592e+02  on 7  degrees of freedom
Residual deviance: 4.1224e-10  on 0  degrees of freedom
AIC: 51.682

Number of Fisher Scoring iterations: 20

> #Do we need to let the odds ratio vary by race of victim?
> (anovaout<-anova(d2,d3))
Analysis of Deviance Table

Model 1: we ~ rdc * lc + rdc * rvc + lc * rvc
Model 2: we ~ rdc * lc * rvc
  Resid. Df Resid. Dev Df Deviance
1         1    0.70074            
2         0    0.00000  1  0.70074
> # Calculate p value
> 1-pchisq(anovaout$Deviance[2],anovaout$Df[2])
[1] 0.4025351
> # Do we need to let the odds ratio vary by race of victim?
> # Compare to deviance test between d2 and d3 above.
> library(DescTools)#For BreslowDayTest
> BreslowDayTest(mm)

	Breslow-Day test on Homogeneity of Odds Ratios

data:  mm
X-squared = 0.38056, df = 1, p-value = 0.5373

> BreslowDayTest(mm,correct=TRUE)

	Breslow-Day Test on Homogeneity of Odds Ratios (with
	Tarone correction)

data:  mm
X-squared = 0.37944, df = 1, p-value = 0.5379

> # Compare to the Wald test for equality of OR
> coef(summary(d3))["rdcWhite defend.:lcLive:rvcWhite victim",]
     Estimate    Std. Error       z value      Pr(>|z|) 
-2.133184e+01  4.224717e+04 -5.049296e-04  9.995971e-01 
> #*******************************************************/
> #* Proportional Mortality: Fit two Poisson regression  */
> #* models, with offset.                                */
> #*******************************************************/
> #**************************************************************/
> #*          Montana Nickel Smelter Example                    */
> #* Causes of death for nickel smelters, from Breslow and Day  */
> #* (1987) Statistical Models in Cancer Research V. II, Inter- */
> #* national Agency for Research on Cancer.  The study is desc-*/
> #* ribled on pp. 349ff.  Variables are age group, calendar    */
> #* period, hire period, arsenic exposure, person-years at     */
> #* risk, deaths from all causes, deaths from respiratory can- */
> #* cer, deaths from all cancers, and deaths from cirulatory   */
> #* disease.                                                   */
> #* http://faculty.washington.edu/norm/BDdata/montana4.dat     */
> #**************************************************************/
> montana<-as.data.frame(scan("montana4.dat",
+    what=list(ageg=0,period=0,hire=0,arsenic=0,pyar=0,alld=0,
+       rcd=0,acd=0,cd=0)))
> montana$lpyar<-log(montana$pyar)
> summary(glm(rcd~factor(period)+factor(hire)+factor(ageg)+
+    arsenic+offset(lpyar),family=poisson,data=montana))$coef
                  Estimate Std. Error    z value      Pr(>|z|)
(Intercept)     -8.3354313 0.30066091 -27.723694 3.617907e-169
factor(period)2  0.5415218 0.21510528   2.517473  1.181999e-02
factor(period)3  0.7124689 0.21414554   3.327031  8.777650e-04
factor(period)4  0.6925791 0.23017823   3.008882  2.622112e-03
factor(hire)2   -0.4828380 0.15266151  -3.162801  1.562590e-03
factor(ageg)2    1.3816528 0.24675608   5.599265  2.152619e-08
factor(ageg)3    2.1692010 0.24439792   8.875693  6.949909e-19
factor(ageg)4    2.3088425 0.26987482   8.555235  1.176290e-17
arsenic          0.3174323 0.05371368   5.909710  3.427098e-09
> summary(glm(cd~factor(period)+factor(hire)+factor(ageg)+arsenic+
+    +offset(lpyar),family=poisson,data=montana))$coef
                   Estimate Std. Error     z value
(Intercept)     -6.07404150 0.12414604 -48.9265836
factor(period)2  0.18866370 0.09564960   1.9724463
factor(period)3  0.22527782 0.09462166   2.3808271
factor(period)4  0.13500748 0.10074529   1.3400873
factor(hire)2    0.02073941 0.07135351   0.2906572
factor(ageg)2    1.05337430 0.09959307  10.5767835
factor(ageg)3    1.80339197 0.09899671  18.2166864
factor(ageg)4    2.67187802 0.10394313  25.7051921
arsenic          0.04452792 0.02952771   1.5080043
                     Pr(>|z|)
(Intercept)      0.000000e+00
factor(period)2  4.855868e-02
factor(period)3  1.727382e-02
factor(period)4  1.802170e-01
factor(hire)2    7.713135e-01
factor(ageg)2    3.818421e-26
factor(ageg)3    3.805112e-74
factor(ageg)4   1.022650e-145
arsenic          1.315534e-01
> #**************************************************************/
> #* Now add the two causes of death, and do prop. mortality.   */
> #**************************************************************/
> summary(glm(cbind(rcd,cd)~factor(period)+factor(hire)+
+    factor(ageg)+arsenic,family=binomial,data=montana))$coef
                  Estimate Std. Error   z value     Pr(>|z|)
(Intercept)     -2.3070781 0.33692841 -6.847384 7.521281e-12
factor(period)2  0.3379935 0.23939707  1.411853 1.579932e-01
factor(period)3  0.4523566 0.23630087  1.914325 5.557867e-02
factor(period)4  0.5553701 0.25439285  2.183120 2.902696e-02
factor(hire)2   -0.4479115 0.16940110 -2.644088 8.191124e-03
factor(ageg)2    0.3242341 0.26783869  1.210557 2.260651e-01
factor(ageg)3    0.3881245 0.26576181  1.460422 1.441740e-01
factor(ageg)4   -0.2910999 0.29026972 -1.002860 3.159284e-01
arsenic          0.2736544 0.06292414  4.348957 1.367868e-05
> #*******************************************************/
> #* Logistic regregression for 2xK table.               */
> #*******************************************************/
> stage4$frx<-factor(stage4$rx)
> summary(glm(alive~frx,family=binomial,data=stage4))$coef
              Estimate Std. Error    z value     Pr(>|z|)
(Intercept) -1.5869651  0.3658393 -4.3378749 0.0000143867
frx2         0.2712883  0.4991369  0.5435148 0.5867754225
frx3         1.0273493  0.4608794  2.2291066 0.0258068132
frx4         0.4383423  0.4849243  0.9039398 0.3660273124
> #* Demonstraton of logistic regression with data grouped.*/
> stage4a<-as.data.frame(xtabs(~rx+alive,data=stage4))
> stage4a$frx<-factor(stage4a$rx)
> summary(glm(alive~frx,family=binomial,weight=Freq,
+    data=stage4a))$coef
              Estimate Std. Error    z value     Pr(>|z|)
(Intercept) -1.5869651  0.3658392 -4.3378762 1.438661e-05
frx2         0.2712883  0.4991369  0.5435148 5.867754e-01
frx3         1.0273493  0.4608793  2.2291070 2.580679e-02
frx4         0.4383423  0.4849243  0.9039397 3.660274e-01
> #**************************************************************/
> #* Stratified 2x2 tables using logistic regression            */
> #* Treat death penalty data as two 2x2 tables for lc vs rvc   */
> #* stratified on rdc, or as two 2x2 tables for lc vs rdc      */
> #* stratified on rvc.  The model without interactions presumes*/
> #* the same OR between lc and each of the stratifiers/cov-    */
> #* variates in both tables, and the coefficient of covariate  */
> #* is the log of the odds ratio.  So the Wald test for each   */
> #* covariate is an alternative to the Mantel-Hanszel test.    */
> #* The test of the interaction is an alternative to the Bres- */
> #* low-Day test for homogeity of odds ratios.  Note infinite  */
> #* MLE for model with interactions.                           */
> #**************************************************************/
> nonprim$surv<-(nonprim$lc=="Live")+0
> nonprim$rvf<-factor(nonprim$rvc,labels=c("White","Black"))
> nonprim$rdf<-factor(nonprim$rdc,labels=c("White","Black"))
> summary(glm(surv~rvf+rdf,
+    family=binomial,data=nonprim,weight=we))$coef
              Estimate Std. Error   z value     Pr(>|z|)
(Intercept)  2.8421051  0.4203344  6.761533 1.365391e-11
rvfBlack    -1.3242128  0.5193427 -2.549786 1.077891e-02
rdfBlack     0.4402222  0.4008880  1.098118 2.721531e-01
> summary(glmout<-glm(surv~rvf*rdf,
+    family=binomial,data=nonprim,weight=we))$coef
                    Estimate   Std. Error     z value
(Intercept)         2.782952    0.4206851  6.61528439
rvfBlack           -1.229603    0.5358319 -2.29475533
rdfBlack           14.242955 1006.6074794  0.01414946
rvfBlack:rdfBlack -13.857940 1006.6075640 -0.01376697
                      Pr(>|z|)
(Intercept)       3.708381e-11
rvfBlack          2.174715e-02
rdfBlack          9.887107e-01
rvfBlack:rdfBlack 9.890159e-01
> aov(glmout)
Call:
   aov(formula = glmout)

Terms:
                      rvf       rdf   rvf:rdf Residuals
Sum of Squares   0.551576  0.133215  0.000627 31.339122
Deg. of Freedom         1         1         1         3

Residual standard error: 3.232085
Estimated effects may be unbalanced
>