# us death rates respca, by age (rows) and period (cols) ,and age weights

rm(all)

us=read.table("usratesRespCa.txt",header=T)
names(us)
rate=us[,2:5]

wts=c(37.4, 30.1, 21.5, 11.0) /100

# all observations

all=read.table("montanadata.txt",header=T)

names(all)

# limit to those hired before 1925

ds = all[all$hire==1,]

# add column of weights

w = wts[ds$a]

ds$w = w

# add column of us death rates

us.rates=c();
for (i in 1:length(ds$a)) {
r=ds[i,1] ; c=ds[i,2];
us.rates=c(us.rates,rate[r,c])
} 

ds$us.rates = us.rates

# ben

tapply(ds$d, factor(ds$exp), sum)


# table 3.3

observed.deaths = rep(0,4);
expected.deaths=rep(0,4);
person.time = rep(0,4);
standardized.rate=rep(0,4);
standardized.us.rate=rep(0,4);
for (k in 1:4){
observed.deaths[k] = sum(ds[ds$exp==k,"d"]);
person.time[k] = sum(ds[ds$exp==k,"py"]);
standardized.rate[k] = 1000*sum(ds[ds$exp==k,"w"]*ds[ds$exp==k,"d"]
                                              /ds[ds$exp==k,"py"]
                                )/
                            sum(ds[ds$exp==k,"w"]);
standardized.us.rate[k] = sum(ds[ds$exp==k,"w"] * ds[ds$exp==k,"us.rates"])/
                          sum(ds[ds$exp==k,"w"]);
expected.deaths[k] = (1/1000)*sum(ds[ds$exp==k,"py"] * ds[ds$exp==k,"us.rates"]);
 };
observed.deaths
person.time
crude.rate = 1000*(observed.deaths / person.time); crude.rate
standardized.rate
standardized.us.rate
cmf.percent = 100*standardized.rate/standardized.us.rate;cmf.percent
expected.deaths
smr=100*observed.deaths/expected.deaths;smr
ratio.smrs=smr/smr[1]; ratio.smrs
adjusted.expected = sum(observed.deaths)*(expected.deaths/sum(expected.deaths));
adjusted.expected
chisq.test.homogeneous.smrs = sum((observed.deaths-adjusted.expected)^2/adjusted.expected);
chisq.test.homogeneous.smrs
exposure.x =1:4
o.minus.e.tilde = observed.deaths-adjusted.expected
chisq.test.trend.num = ( sum(exposure.x*o.minus.e.tilde) )^2;
chisq.test.trend.den = sum( (exposure.x^2)*adjusted.expected );
chisq.test.trend.den = chisq.test.trend.den - 
   (sum(exposure.x*adjusted.expected))^2/sum(observed.deaths)
chisq.test.trend.smrs = chisq.test.trend.num/chisq.test.trend.den;
chisq.test.trend.smrs