#EXERCISES 1: Solution to Q3

#Numerical studies

#First some calculation functions
#Function pi.est computes the estimate of Pi from a vector 
#x where x[1] = "Number of cases", x[2] = "Number of non-cases"

pi.est<-function(x){x[1]/(x[1]+x[2])}

#Function relative.risk computes the relative risk (or risk ratio) estimate pi1/pi0
#that is P(F|E)/P(F|E') 
relative.risk<-function(x){(x[1,1]/(x[1,1]+x[2,1]))/(x[1,2]/(x[1,2]+x[2,2]))}

#Function odds.ratio.cases computes the odds on exposure in cases
#that is P(E|F)/P(E'|F) 
odds.ratio.cases<-function(x){x[1,1]/x[1,2]}

#Function odds.ratio.controls computes the odds on exposure in controls
#that is P(E|F')/P(E'|F') 
odds.ratio.controls<-function(x){x[2,1]/x[2,2]}

#Function odds.ratio computes the estimate 
#pi1/(1-pi1)
#-------------
#pi0/(1-pi0)
odds.ratio<-function(x){(x[1,1]*x[2,2]/(x[1,2]*x[2,1]))}


##############################################################################

#Q3

table3i.cohort<-cbind(c(210,90),c(120,180))
table3ii.cohort<-cbind(c(230,70),c(120,80))
table3iii.cohort<-cbind(c(170,130),c(80,2200))
table3iv.cohort<-cbind(c(120,180),c(130,170))

#(i)
rd.i_pi.est(table3i.cohort[,1])-pi.est(table3i.cohort[,2])
rr.i_relative.risk(table3i.cohort)
or.i_odds.ratio(table3i.cohort)

print(paste("Risk Difference Table (i) =",round(rd.i,3)))
print(paste("Relative Risk Table   (i) =",round(rr.i,3)))
print(paste("Odds Ratio Table      (i) =",round(or.i,3)))

#(ii)
rd.ii_pi.est(table3ii.cohort[,1])-pi.est(table3ii.cohort[,2])
rr.ii_relative.risk(table3ii.cohort)
or.ii_odds.ratio(table3ii.cohort)

print(paste("Risk Difference Table (ii) =",round(rd.ii,3)))
print(paste("Relative Risk Table   (ii) =",round(rr.ii,3)))
print(paste("Odds Ratio Table      (ii) =",round(or.ii,3)))

#(iii)
rd.iii_pi.est(table3iii.cohort[,1])-pi.est(table3iii.cohort[,2])
rr.iii_relative.risk(table3iii.cohort)
or.iii_odds.ratio(table3iii.cohort)

print(paste("Risk Difference Table (iii) =",round(rd.iii,3)))
print(paste("Relative Risk Table   (iii) =",round(rr.iii,3)))
print(paste("Odds Ratio Table      (iii) =",round(or.iii,3)))

#(iv)
rd.iv_pi.est(table3iv.cohort[,1])-pi.est(table3iv.cohort[,2])
rr.iv_relative.risk(table3iv.cohort)
or.iv_odds.ratio(table3iv.cohort)

print(paste("Risk Difference Table (iv) =",round(rd.iv,3)))
print(paste("Relative Risk Table   (iv) =",round(rr.iv,3)))
print(paste("Odds Ratio Table      (iv) =",round(or.iv,3)))


#######################################################################################################

#Q4(i)  Assuming a cohort study: look at relative risks

table4i0<-cbind(c(70,30),c(80,120))
table4i1<-cbind(c(160,40),c(40,60))

rr.4i0_relative.risk(table4i0)
rr.4i1_relative.risk(table4i1)

print(paste("Relative Risk (X=0) :",round(rr.4i0,3)))
print(paste("Relative Risk (X=1) :",round(rr.4i1,3)))

#Is X related to exposure ?  Need to check the E vs X margin (combining over F and F')
#Doing this we get the table
#
#     E   E'
#X=0 100 200
#X=1 200 100
#
#and then look at, for example, P(X=0|E) and P(X=0|E'), or the odds ratio

table4iEX<-cbind(c(sum(table4i0[,1]),sum(table4i0[,2])),c(sum(table4i1[,1]),sum(table4i1[,2])))

table4iEX.OR<-(table4iEX[1,1]*table4iEX[2,2])/(table4iEX[1,2]*table4iEX[2,1])

#If the OR is not equal to 1, then X and E are related.

#Actually, we can check this formally using the STANDARD ERROR for the OR
#and construct a 95 % CI (on the log scale)

table4iEX.ORSE<-sqrt(sum(1/table4iEX))
table4iEX.95CI<-exp(c(log(table4iEX.OR)-1.96*table4iEX.ORSE,log(table4iEX.OR)+1.96*table4iEX.ORSE))

print(table4iEX.95CI)
if(table4iEX.95CI[1] > 0 & table4iEX.95CI[2] > 0){
	print("X is related to E")
}
else if(table4iEX.95CI[1] < 0 & table4iEX.95CI[2] < 0){
	print("X is related to E")
}
else{
	print("X is NOT related to E")
}



#(ii)
table4ii0<-cbind(c(90,10),c(60,40))
table4ii1<-cbind(c(80,120),c(20,180))
rr.4ii0_relative.risk(table4ii0)
rr.4ii1_relative.risk(table4ii1)
print(paste("Relative Risk (X=0) :",round(rr.4ii0,3)))
print(paste("Relative Risk (X=1) :",round(rr.4ii1,3)))
table4iiEX<-cbind(c(sum(table4ii0[,1]),sum(table4ii0[,2])),c(sum(table4ii1[,1]),sum(table4ii1[,2])))
table4iiEX.OR<-(table4iiEX[1,1]*table4iiEX[2,2])/(table4iiEX[1,2]*table4iiEX[2,1])
table4iiEX.ORSE<-sqrt(sum(1/table4iiEX))
table4iiEX.95CI<-exp(c(log(table4iiEX.OR)-1.96*table4iiEX.ORSE,log(table4iiEX.OR)+1.96*table4iiEX.ORSE))
print(table4iiEX.95CI)
if(table4iiEX.95CI[1] > 0 & table4iiEX.95CI[2] > 0){
	print("X is related to E")
}
else if(table4iiEX.95CI[1] < 0 & table4iiEX.95CI[2] < 0){
	print("X is related to E")
}
else{
	print("X is NOT related to E")
}

#############################################################################################
#Q5
table5<-cbind(c(96,109),c(104,666))

table5.OR_odds.ratio(table5)

table5.ORSE<-sqrt(sum(1/table5))
table5.95CI<-exp(c(log(table5.OR)-1.96*table5.ORSE,log(table5.OR)+1.96*table5.ORSE))

table5.OR
table5.ORSE
table5.95CI