# --------- Ayas et al data on PI's in Residents --------

#### Incidence, Overall ####

cases=c(498); InternMonths=c(17003);

## regression through origin, Gaussian variation around mean

summary( lm(cases ~ -1+InternMonths) )



## regression through origin, Poisson variation around mean

summary( glm(cases ~ -1+InternMonths, family=poisson(link=identity)) )


## Model log[rate], and Poisson variation around expected no. of cases

summary( glm(cases ~ 1, family=poisson, offset=log(InternMonths) ) )


#### Incidence, extened vs non-extended ####

cases=c(35,46); PT=c(26667,60763); e=c(1,0); ePT = e*PT; log.o = log(PT)
ds <- data.frame(cases=cases, opps=PT, extended=e, extended.opps = ePT, log.opps=log.o); 
ds
attach(ds)

# regression to obtain rate ratio

summary(glm(cases ~ -1 + opps + extended.opps, family=poisson(link=identity)))

# regression to obtain rate difference

summary(glm(cases ~ extended, family=poisson(link=log),offset=log.opps))


 #     ---------- Quebec Death rates ------------
 
 ds=read.table("quedata.txt",header=T)
 
 # death rates '91 vs '71 
 
 sum(ds$deaths[ds$year==1991])/sum(ds$population[ds$year==1991])
 sum(ds$deaths[ds$year==1971])/sum(ds$population[ds$year==1971])
 
 
 ds7191=ds[(ds$age > 40) & (ds$age < 85 ) & ( (ds$year==1971) | (ds$year==1991)),]
 
 attach(ds7191)
 
 # age and sex specific death rates '91  
 
y91=ds[(ds$age > 40) & (ds$age < 85 ) & (ds$year==1991) ,]
y91$age=y91$age - 40

y91m=y91[y91$male==1,]; y91f=y91[y91$male==0,];
(y91m$deaths/y91m$population) / (y91f$deaths/y91f$population)
mean( (y91m$deaths/y91m$population) / (y91f$deaths/y91f$population) )

plot(y91$age,  log( y91$deaths / y91$population ) )

summary( glm(deaths ~ age + male, family=poisson, offset=log(population),data=y91) )

summary( glm(deaths ~ age + male + age*male, family=poisson, offset=log(population),data=y91) )

m7191=ds7191[male==1,]; f7191=ds7191[male==0,]

plot(f7191$age,  log( f7191$deaths / f7191$population ) )



# ------------ Uganda rct of circumcision to prevent hiv infection --


Male circumcision for HIV prevention in men in Rakai,
Uganda: a randomised trial The Lancet Feb 25 2007

t.from t.to intervention participants events   py 
   0    6        1            2263      14   1172.1
   0    6        0            2319      19   1206.7
   6   12        1            2235       5   1190.7
   6   12        0            2229      14   1176.3
  12   24        1             964       3    989.7 
  12   24        0             980      12   1008.7
  
ds=read.table("UgandaTrial.txt",header=T); ds

fit= glm(events ~ intervention, family=poisson, offset=log(py),data=ds) ;
summary(fit)
0.01*round(100*exp(fit$coefficients))

0.1*round(10*fit$fit)

attach(ds); after6 = (t.from >= 6)

fit= glm(events ~ intervention + after6, family=poisson, offset=log(py),data=ds) ;
summary(fit)

fit= glm(events ~ intervention + + after6 + intervention*after6, family=poisson, offset=log(py),data=ds) ;
summary(fit)


fit= glm(events ~ intervention + as.factor(t.from), family=poisson, offset=log(py),data=ds) ;
summary(fit)

fit= glm(events ~ intervention + as.factor(t.from) + intervention * as.factor(t.from),
  family=poisson, offset=log(py),data=ds) ;
summary(fit)
exp(fit$coefficients)
exp(fit$coefficients)

# ------------ Clayton Hills Ch 22 regression --


Table 22.6   -------- R analysis

   age  exposed events  py 
   45      0      4   607.9
   45      1      2   311.9
   55      0      5  1272.1
   55      1     12   878.1
   65      0      8   888.9 
   65      1     14   667.5
  
ds=read.table("chtable226.txt",header=T); ds

fit= glm(events ~ exposed+as.factor(age), family=poisson, offset=log(py),data=ds) ;
summary(fit)
0.01*round(100*exp(fit$coefficients))

0.1*round(10*fit$fit)
ds$events


Table 22.6 ---------- SAS analysis 

data a;
input   age  exposed events  py; log_py  = log(py);
lines;
        45      0      4   607.9
        45      1      2   311.9
        55      0      5  1272.1
        55      1     12   878.1
        65      0      8   888.9
        65      1     14   667.5
;

proc genmod data=a; * by default uses
 class age ; * early sas versions used highest value as default ref. category;
 model events = exposed age / dist=poisson offset=log_py ;
 run;
 
 * can use ods options to see the fitted values;
 
 
 
 Table 22.6 ---------- Stata analysis 
 
* Clayton Hills Table 22.6
 
clear
input   age  exposed events  py

        45      0      4   607.9
        45      1      2   311.9
        55      0      5  1272.1
        55      1     12   878.1
        65      0      8   888.9
        65      1     14   667.5
        end

generate logpy  = log(py)

xi i.age

* beta_hats 

glm events exposed _Iage_55 _Iage_65, family(poisson) link(log) offset(logpy)

* exp[beta_hats] 

glm events exposed _Iage_55 _Iage_65, eform family(poisson) link(log) offset(logpy)

* fitted no. of cases (added to dataset)

predict e_cases