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

ds=read.table("TxYears.txt") ;
ds$t = ds$Year - 0.5; ds$t[9:10]=c(4.25,4.25)  ; 

# observed rates

ds$obs.rate = round(ds$events/ds$P.Years,4) ; ds

# plots of these

plot(ds$t[ds$Tx==0],ds$obs.rate[ds$Tx==0],col="red",cex=2,
 ylim=c(-0.01,0.07),xlim=c(0,5),pch=19,
 ylab="Rate (events/PY)", xlab="Year")
segments(0,0,5,0)
points(ds$t[ds$Tx==1],ds$obs.rate[ds$Tx==1],col="blue",pch=19,cex=2)

# lines fitted to observed rates 

# choose 1 weighting scheme 1 by commenting out the other one ..
# come back and set the weights to the other, and re-fit and plot..

  ds$w=ds$P.Years  # weights = PT
  ds$w=rep(1,10);  # weights of 1 each

fit=lm(obs.rate ~ t + Tx, data=ds, weights=w)
str(fit) # objects
lines(ds$t[ds$Tx==0],
  fit$fitted.values[ds$Tx==0],pch=20,col="red",lwd=2)
lines(ds$t[ds$Tx==1],
  fit$fitted.values[ds$Tx==1],pch=20,col="blue",lwd=2)
  
# go back and set the other set of weights... refit and redraw
  
 # ============  Poisson Regression models (# Control arm only...) ======

ds.noTx = ds[ds$Tx==0,3:5] ; 

# for MULTIPLICATIVE MODELS...

ds.noTx$log.PY = log(ds.noTx$P.Years)

ds.noTx

mult.model.fit=glm(events ~ t,
                   offset=log.PY,
                   family=poisson(link="log"),
                   data=ds.noTx)

summary(mult.model.fit)

str(mult.model.fit)

exp(mult.model.fit$coefficients)

exp(mult.model.fit$linear.predictors)

mult.model.fit$fitted.values

# plot
Time=seq(0,5,0.1);
fitted.ID.mult = exp(
  mult.model.fit$coefficients[1] + mult.model.fit$coefficients[2]*Time
               )
plot(Time,fitted.ID.mult, ylim=c(0,0.04), type="l")

# etc..

# EXTEND to the 2-sample case, with tx=1 and tx=0 ...
# (left as part of exercise..)

# === for ADDITIVE RATE MODELS...

ds.noTx$P.Years.times.t = ds.noTx$P.Years * ds.noTx$t; ds.noTx

add.model.fit = glm(events ~ -1 + P.Years + P.Years.times.t,
                    family=poisson(link="identity"),
                    data=ds.noTx)
                     
summary(add.model.fit)

str(add.model.fit)

add.model.fit$fitted.values

add.model.fit$fitted.rates = add.model.fit$fitted.values / ds.noTx$P.Years
add.model.fit$fitted.rates

# plots

Time=seq(0,5,0.1);
fitted.ID.add = add.model.fit$coefficients[1] + add.model.fit$coefficients[2]*Time
plot(Time,fitted.ID.add, ylim=c(0,0.04), type="l",col="red")

lines(Time,fitted.ID.mult, ylim=c(0,0.04), type="l",col="blue")

text(4,0.021,"Additive",col="red")
text(4,0.032,"Multiplicative",col="blue")


# etc...

# Again, extension to model involving tx left as exercise