"read.f"<-
function(file = "gdata.txt")
{  #  use "File Save As" to download ascii file  gdata.txt
   #  use "File Save As" to download this ascii file ( sp-tabAIV.txt )
      # run Splus and
   # > source("sp-tabAIV.txt")
   # > read.f()
   # > tab.aiv<- tabaiv.f()
   #     note: at this point
   # "dt" is data frame with all cancer data from read.f
   # "tab.aiv" is output object from function tabaiv.f
   # type the next line to wirte out Table AIV into file "TableAIV.txt" 
   # > sink("TableAIV.txt");Info; tab.aiv; sink()
# read.f 
# Input: ADS for All Cancer data as Described in Appendix ( gdata.txt)
# Output:  data frame dt for Poisson Regresion
#
	varb <- c("age", "f", "XG", "IG", "B", "S", "L", "y1", "e1", "Dmrem", 
		"PY")
	d <- read.table(file, col.names = varb)	# read all Ca data
#
#             define variables and factors for Poisson regression
#               see Statistical Methods   Eqs. 1 and 3 
#
	ra <- d$y1/d$e1
	ageS <- (d$age - 52.5)/100
	D <- d$Dmrem/100	# dose in Sv
#
# define factors F XG IG B S and L
#
	d$f <- factor(d$f, levels = c(1, 2, 3), labels = c("X-10", "Y-12", 
		"MULT"))	#
	d$XG <- factor(d$XG, levels = c(1, 2, 3, 4, 5, 6, 7, 8, 9, 10), labels
		 = c("zero", ".01-5", "5-10", "10-", "20-", "40-", "80-", 
		"160-", "320-", "640+"))	#
	d$IG <- factor(d$IG, levels = c(1, 2, 3), labels = c("EN", "EM", "NE"))
	d$S <- factor(d$S, levels = c(1, 2), labels = c("NonMth", "Mth"))
	d$B <- factor(d$B, levels = c(1, 2, 3, 4, 5), c("< 1900", "1900-09", 
		"1910-19", "1920-29", "1930+"))
	d$L <- factor(d$L, levels = c(1, 2), labels = c("1 yr+", "<1yr"))	
	# output data frame dt
	dt <<- data.frame(d, ra, ageS, D)	#
#  Person Years panel-1 of  Table AIV  
	p1aiv <- round(tapply(d$PY, list(d$f, d$XG), sum))
	p1aiv
}
"tabaiv.f"<-
function(d = dt)
{
#  tabaiv.f is Splus Function for Oak Ridge Workers Mortality Study ORNL-6785
#      Input:  data frame dt from Splus function read.f()
#       Output: Table AIV of ORNL-6785 ( Radiation Research Table AI)
#
	p1 <- round(tapply(d$PY, list(d$f, d$XG), sum))
	TotPY <- round(tapply(d$PY, list(d$XG), sum))	#
# panel 2 Observed and Expected ( U.S. W males) Facility by Dose Group
#
	po <- tapply(d$y1, list(d$f, d$XG), sum)
	pe <- round(tapply(d$e1, list(d$f, d$XG), sum), 2)
	ty <- tapply(d$y1, list(d$XG), sum)
	te <- round(tapply(d$e1, list(d$XG), sum), 2)
	smr <- round((100 * ty)/te, 2)
	p2 <- rbind(po[1,  ], pe[1,  ], po[2,  ], pe[2,  ], po[3,  ], pe[3,  ])
	Total <- apply(p2, 1, sum, na.rm = T)
	p2 <- cbind(p2, Total)
	dimnames(p2)[[1]] <- c("X-10_Obs", "X-10_Exp", "Y-12_Obs", "Y-12_Exp", 
		"MULT_Obs", "MULT_Exp")
	p2t <- rbind(ty, te, smr)
	dimnames(p2t)[[1]] <- c("ToT_Ob", "ToT_Ex", "SMR")
	aiv.p2 <- p2
	pt <- sum(d$PY)
	out <- list(Pyrson.Years = p1, Total.PY = TotPY, Obs.Exp = p2, Total.OE
		 = p2t)	#
#  panel 3 Adjust Expected Deaths
#  fit main effects model with all terms in Table AV except Dose (XG)
#  save glm output in meff 
	options(contrasts = "contr.treatment")
	meff <<- glm(formula = ra ~ ageS + B + S + L + IG + f - 1, weights = e1,
		family = poisson, weight = e1, data = d)
	adjy1 <- d$e1 * fitted(meff)
	p3 <- round(tapply(adjy1, list(d$f, d$XG), sum), 2)
	Total <- round(tapply(adjy1, list(d$XG), sum), 2)
	Adj.smr <- round((100 * ty)/Total, 2)
	p3t <- rbind(Total, Adj.smr)
	out <- list(Pyrson.Years = p1, Total.PY = TotPY, Obs.Exp = p2, Total.OE
		 = p2t, Adj.Exp = p3, Total.Adj = p3t)
	out	# output most of Table AIV
}
Info<-c("ORNL-6785  Table AIV Summary For X-10/Y-12 White Males Ten Year Lag") 
# end of sp-tabAIB.txt