rm(list=ls(all=TRUE)) cols=c("city","sex","doseg","agexg","calg","kerma","PY","adjPY","num.entering", "age","agex","tsx","cal","sv","gam","neut","lymphoma","NHL","leukemia","AML","ALL","CML","ATL","MML") d<-read.table("/data/abomb/HEMA87.DAT", header=F,col.names=cols); d=d[d$adjPY>0,] #remove any recs with zero py d=d[d$kerma==1,] # take only kerma < 4 Gy attach(d); #summary(d) svc=cut(sv,c(-1,.01,.4,10),labels=c("low","med","high")) agec=cut(age,seq(0,80,10)) adjPY=adjPY*1e4 BKPY=tapply(adjPY,list(sv=svc,age=agec),sum,na.rm=TRUE) # pool across city and sex to get background expected values BKPY sum(BKPY) PY=tapply(adjPY,list(sv=svc,age=agec,sex=sex,city=city),sum,na.rm=TRUE) # 1 = Hiroshima; 2 = Nagasaki; 1 = Male; 2 = Female dimnames(PY)$sex<-c("males","females") dimnames(PY)$city<-c("Hiroshima","Nagasaki") PY sum(PY) Mat=list() for (h in c("NHL","AML","ALL","CML","MML") ) { # h="CML" print(h) canc=d[,h] bkCounts=tapply(canc,list(sv=svc,age=agec),sum,na.rm=TRUE) # pool across city and sex to get background expected values bkCounts counts=tapply(canc,list(sv=svc,age=agec,sex=sex,city=city),sum,na.rm=TRUE) dimnames(counts)$sex<-c("males","females") dimnames(counts)$city<-c("Hiroshima","Nagasaki") counts ntsx=tsx[canc>0] # eliminate cells without cancer so as to NOT incorporate empty cell tsx's into averages nsvc=svc[canc>0] nagec=agec[canc>0] nsex=sex[canc>0] ncity=city[canc>0] TSX=tapply(ntsx,list(sv=nsvc,age=nagec,sex=nsex,city=ncity),mean,na.rm=TRUE) dimnames(TSX)$sex<-c("males","females") dimnames(TSX)$city<-c("Hiroshima","Nagasaki") bkRate=bkCounts[1,]/BKPY[1,] # first row is low dose group which is used here as the background/control group TSX #note the short times to cancer in the high dose group E=array(dim=c(3,8,2,2),dimnames=list(sv=c("low","med","high"),age=c(levels(agec)),sex=c("males","females"),city=c("Hiroshima","Nagasaki"))) for (i in c("low","med","high")) for (k in c("males","females")) for (l in c("Hiroshima","Nagasaki")) E[i,,k,l]=bkRate*PY[i,,k,l] E TOE=array(dim=c(9,8,2,2),dimnames=list(sv=c("lowTSX","lowO","lowE","medTSX","medO","medE","highTSX","highO","highE"),age=c(levels(agec)),sex=c("males","females"),city=c("Hiroshima","Nagasaki"))) #print(TOE) for (i in c("low","med","high")) for (k in c("males","females")) for (l in c("Hiroshima","Nagasaki")) for (m in c("TSX","O","E")) { if (m=="TSX") TOE[paste(i,m,sep=""),,k,l]=TSX[i,,k,l] if (m=="O") TOE[paste(i,m,sep=""),,k,l]=counts[i,,k,l] if (m=="E") TOE[paste(i,m,sep=""),,k,l]=E[i,,k,l] } Mat[[h]]=TOE } Mat detach(d)