#               APPENDIX 2 COMPUTATIONAL DETAILS
#
# POISSON REGRESSION ANALYSIS OF ILLNESS AND INJURY SURVEILLANCE DATA
# Technical Report ORISE-09-OEWH-0176 (2012)   short name (PRA_IISP)
#    Available at URL  http://www.csm.ornl.gov/~frome/priisp
#    or mirror site http://home.comcast.net/~fromestat/priisp
#                File Name readdata.txt at Step 7 (this file)
#                To obtain results at Step 7 of website
# Line 9         RUN R and open a txt file using sink(), i.e.
#                sink("PRA_IISPdemo.txt") 
#                then  SOURCE THIS FILE WITH echo=TRUE, i.e.
# line 12   source("readdata.txt",echo=TRUE,max.deparse.length =3000) 
# 
#  Input data frame drtw with IISP RTW data from binary file drtw.rda
           load("./RFIISP/drtw.rda")
#    with column names  
#    "yr","fa","sx","og","ag","pyr","all","rsp","inj","hrt","msc")
#
#  see technical report PRA_IISP in Item 1 for additional Information
 
#     "yr","fa","sx","og","ag"       are factor covariates
#     "all","rsp","inj","hrt","msc"  are counts of RTW health events 
#
#####################################################################
#
#  The folder RFIISP contains utility functions that were used to 
#  obtain some of the numerical results in PRA_IISP--- see Example 1
#  NOTE:
#  THESE FUNCTIONS ARE PROVIDED FOR THE READERS CONVENIENCE AND
#  ARE SPECIFIC TO THE IISP DATA STRUCTURES AND METHODS
#--------------------------------------------------------------------
  source("./RFIISP/read.RF.R")
  ogn7<-read.RF()# input R functions and Occupational Group long names
  library(MASS)   # attach MASS libraay
#####################################################################
#
#     calculate marginal table of rsp events, person years,
#     and rsp event rates for Occupational and  Age Groups
#     for SNL--- see Table 1 in PRA_IISP report
#---------------------------------------------------------------------
   evcount.oa(d=drtw,rfa="SNL",ryr="ALL",event="rsp",pc=F,ognames=ogn7)

######################################################################
#
#     Fit the Main Effects Model(MEM) to RTW rsp event data
#     --- see Table 2 in PRA_IISP Technical Report
#---------------------------------------------------------------------
  h <- make.h(drtw,rfa="SNL",event="rsp")
  make.h  #  this function extracts an event at facility rfa

######################################################################
#
#  fit the MEM using R function glm()--- see glm help file for details
#
#--------------------------------------------------------------------
     mf<- glm(y ~ ag+sx+og+yr,fa=poisson,data=h,offset= log(nn) )
     summary(mf)

#####################################################################
#
#  calculate analysis of deviance table for facility SNL rsp events
#            see Table 3 in PRA_IISP report 
#---------------------------------------------------------------------
     aodme4(drtw,"SNL")

#####################################################################
#
#      Use utility function compare.estimates() --- see listing below
#      TO COMPARE PARAMETER ESTIMATES AND STANDARD ERRORS
#      FOR FOUR METHODS BASED ON IWLS ALGORITHM
# 
#      1) Poisson ML    2) QL/M1     3) NB ML and    4) QL/M2
#
#      for IISP data for facility= "PTX"  event= "rsp"
#
#      THIS IS EXAMPLE 2 IN TECHNICAL REPORT PRA_IISP.
#-------------------------------------------------------------------
    compare.estimates(drtw,rfa="PTX",event="rsp",mf=NA)

#####################################################################
#
#      List R utility function compare.estimates()
#-------------------------------------------------------------------
   compare.estimates
#
#####################################################################
#
#    Close file PRA_IISPdemo.txt and exit R Dec 1 2012
#----------------------------------------------------------------------
     sink()

######################################################################
#
#THE OUTPUT WILL APPEAR IN FILE   PRA_IISPdemo.txt in R working folder 
#  
#  End of File readdata.txt
#####################################################################