########## full dataset


ds=read.table("fruitflyDataR.txt",header=T)

# 1 partner

ds.1=ds[ds$partners==1,]

lm(longevity~active,data=ds.1)

require(survival)

coxph(Surv(longevity)~active,data=ds.1)



# added, march 16

attach(ds.1);
thorax.z = (thorax-mean(thorax))/sd(thorax); summary(thorax.z)

fit = coxph( Surv(longevity) ~ active + thorax.z )


# estimate curve for active=0, thorax.z = 0 .. *** from scratch ***

  ## create no. of zero-equivalents
  
   hr.hat=exp(fit$coefficients[1]*active + fit$coefficients[2]*thorax.z)

  ## get sizes of risksets, using each weight in portfolio
  
   n.risk=summary(survfit(Surv(longevity),weights=hr.hat))$n.risk; n.risk ;
   
  ## get numbers of deaths at each time
   
   n.event=summary(survfit(Surv(longevity)))$n.event ; n.event ;
   
  ## get times of the events
    
   event.times=summary(survfit(Surv(longevity)))$time ; event.times;

  ## nelson-aalen survival estimate
  
    nelson.aalen.S = exp(-cumsum(n.event/n.risk)); nelson.aalen.S;
    plot(event.times,nelson.aalen.S)

# estimate S curve for active=0, thorax.z = 0 .. using built in routines...

ph.fit <- coxph( Surv(longevity) ~ active + thorax.z)

 ## fitted curves for 2 profiles, active = 0, thorax.z = 0 and active = 0, thorax.z = 1 
 
curves = survfit( ph.fit, newdata=data.frame(active = c(0,0),thorax.z = c(0,1) ) ) 
 
plot( c(0,curves$time),c(1,curves$surv[,1]),type="l")
lines(c(0,curves$time),c(1,curves$surv[,2]),col="red")

# for interest .. km-curve for those active = 0

  give weight of 0 to (ie exclude) active ones by using 1-active as weight
  
  km.fit.inactive = survfit(Surv(longevity),weights=(1-active))
  summary(km.fit.inactive)
  plot(km.fit.inactive)
  
  
s2=survfit( ph.fit, newdata=data.frame(active = c(0,0),thorax.z = c(0,1))  )
  
  
Stratified analyses

thorax.stratum = 1 + ( thorax.z > -.67) + ( thorax.z > 0) + ( thorax.z > 0.67 );


ph.fit.stratified <- coxph( Surv(longevity) ~ active + strata(thorax.stratum));
summary(ph.fit.stratified)




######## subset of 10 (Fig1 and Fig3)


# usual way

require(survival)
longevity = c(30,35,40,45,50,65,70,75,80,90);
dead      = c( 1, 1, 1, 1, 1, 1, 1, 1, 1, 1);
active    = c( 1, 0, 1, 0, 1, 1, 0, 1, 0, 0);
coxph(Surv(longevity,dead) ~ active)


# directly  n's in each riskset

n0          = c(5,5,4,4,3,3,3,2);
n1          = c(5,4,4,3,3,2,1,1);

# was case in active group?

case.active = c(1,0,1,0,1,1,0,1); 

# mh

mh=sum(case.active*(n0-(1-case.active))/(n0+n1))/
   sum((1-case.active)*(n1-case.active)/(n0+n1))

### likelihood and log.likelihood functions 'from scratch'   
   
likelihood  =function(B){prod( (exp(B*case.active))/(n0 + exp(B)*n1) ) }
log.likelihood=function(B){sum( log((exp(B*case.active))/(n0 + exp(B)*n1) ) ) }

# range of (trial) coefficient over which to to calculate likelihood

# modify seq(uence) as needed...  

#  syntax is seq(starting values, ending value, increment)

b = seq(-5,5,1); 

lik = lapply(b,likelihood); plot(b,lik);

log.lik = lapply(b,log.likelihood); plot(b,log.lik);

mh=sum(case.active*(n0-(1-case.active))/(n0+n1))/sum((1-case.active)*(n1-case.active)/(n0+n1))