# nov 6, 2011

setwd("/Users/jameshanley/Documents/Courses/634/VasectomyMI")

ds=read.table("vas-mi-data.txt",header=T) ; 

ds$non.v.mi = 1- ds$v.mi ; ds[1:10,]

# outcome .. among all pairs... 
with(ds, table(non.v.mi, v.mi) ) ; 

# among pairs where neither smoked
with(ds[ds$non.v.smoker==0 & ds$v.smoker==0,], 
  table(non.v.mi, v.mi) ) ;
  
# among pairs where both smoked
with(ds[ds$non.v.smoker==1 & ds$v.smoker==1,], 
  table(non.v.mi, v.mi) ) ;
  
# among smoking-matched pairs
with(ds[ds$non.v.smoker == ds$v.smoker,], 
  table(non.v.mi, v.mi) ) ;
  
# where non.vasectomized man smoked, vasectomized man did not
with(ds[ds$non.v.smoker==1 & ds$v.smoker==0,], 
  table(non.v.mi, v.mi) ) ;
  
# where non.vasectomized man didn't smoke, vasectomized man did
with(ds[ds$non.v.smoker==0 & ds$v.smoker==1,], 
  table(non.v.mi, v.mi) ) ;
  
# fully  matched
with(
ds[ds$non.v.smoker == ds$v.smoker & ds$non.v.obese == ds$v.obese,],
  table(non.v.mi, v.mi) )

# M-H HR, matched on smoking...

length(ds[ds$non.v.smoker == ds$v.smoker,1])
with(
ds[ds$non.v.smoker == ds$v.smoker,],
  sum( v.mi*(1-non.v.mi) ) / sum( (1-v.mi)*non.v.mi) ) 

# M-H HR, matched on smoking and obesity...

length(ds[ds$non.v.smoker == ds$v.smoker & ds$non.v.obese == ds$v.obese,1])
with(
ds[ds$non.v.smoker == ds$v.smoker & ds$non.v.obese == ds$v.obese,],
  sum( v.mi*(1-non.v.mi) ) / sum( (1-v.mi)*non.v.mi) ) 


# for conditional logistic regression .. 1 record per man

# from scratch

lik=function(V,S,O) { 
v.tickets =   V * (O^ds$v.obese)     * (S^ds$v.smoker);
non.v.tickets =   (O^ds$non.v.obese) * (S^ds$non.v.smoker);
p = v.tickets/(v.tickets+non.v.tickets);
l = (p^ds$v.mi) * ( (1-p)^ds$non.v.mi );
return(sum(log(l))) }
lik(1,1,1) ; lik(1.1,4.0,3.1)

V.range=2^seq(-2,2,0.25)
O.range=2^seq(-2,2,0.25)
S.range=2^seq(-2,2,0.25)
for (V in V.range){
	for (S in S.range){
		for ( S in S.range ) {
			    print(round(c(V,S,O,lik(V,S,O)),2)) }
			}
}

library(survival)

no=1:36
dv = cbind(no, rep(1,36),ds[,c(3,4,2)]); 
  names(dv)=c("pair.no","v","obese","smoker","mi") ; dv
dn = cbind(no, rep(0,36),ds[,c(5,6,7)]); 
  names(dn)=c("pair.no","v","obese","smoker","mi") ; dn
d=rbind(dv,dn) ; 
ord=c(); for (i in 1:36) ord=c(ord,i,36+i); d=d[ord,] ; 
d[1:6,] ; d[67:72,] 

summary(clogit(mi ~ v                  + strata(pair.no),data=d))

summary(clogit(mi ~ v + obese + smoker + strata(pair.no),data=d))

d$age = 58 # arbitrary

summary(coxph(Surv(age,mi) ~ v + obese + smoker + strata(pair.no),data=d))



with(ds, table(non.v.mi, v.mi,non.v.smoker,v.smoker,non.v.obese,v.obese) )