x<-c(189,104)
n<-c(11034,11037)
p<-prop.test(x,n)
p1<-p$estimate[1]
p2<-p$estimate[2]
D<-p1-p2
D
sigma<-sqrt(p1*(1-p1)/n[1]+p2*(1-p2)/n[2])
sigma
U<-(p1-p2)+qnorm(1-0.025)*sigma
L<-(p1-p2)-qnorm(1-0.025)*sigma
L
U
######Relative Risk#########
RR<-p1/p2
RR
L1<-exp((log(p1/p2))+qnorm(1-0.025)*sqrt((1-p1)/(n[1]*p1)+(1-p2)/(n[2]*p2)))
U1<-exp((log(p1/p2))-qnorm(1-0.025)*sqrt((1-p1)/(n[1]*p1)+(1-p2)/(n[2]*p2)))
L1
U1
#######Odds Ratio########
R1<-(x[1]*(n[2]-x[2]))/(x[2]*(n[1]-x[1]))
R1
ASE<-sqrt(1/x[1]+1/x[2]+1/(n[2]-x[2])+1/(n[1]-x[1]))
L2<-exp(log(R1)+qnorm(1-0.025)*ASE)
U2<-exp(log(R1)-qnorm(1-0.025)*ASE)
L2
U2
############ Eteghad be zendegi ##################
x<-factor(rep(c(1,2),c(625,502),lables=c("Female","Male")))
y<-factor(rep(c(0,1,0,1),c(116,509,104,398),lables=c("Yes","No"))) 
table(x,y) 
chisq.test(x,y) 
chisq.test(x,y)$observed
chisq.test(x,y)$expected
chisq.test(x,y)$residuals
###################### Asprin Example ###################
x<-factor(rep(c(1,2),c(676,648),lables=c("Asprin","placebo"))) 
y<-factor(rep(c(0,1,0,1),c(658,18,620,28),lables=c("Yes","No"))) 
table(x,y) 
chisq.test(x,y) 
chisq.test(x,y)$observed
chisq.test(x,y)$expected
chisq.test(x,y)$residuals
################## Gerayesh Siasi######################
x<-factor(rep(c(1,2),c(577,403),lables=c("Female","Male"))) 
y<-factor(rep(c(1,2,3,1,2,3),c(279,73,225,165,47,191),lables=c("Democrat","Independent","Republican"))) 
table(x,y) 
chisq.test(x,y) 
chisq.test(x,y)$observed
chisq.test(x,y)$expected
chisq.test(x,y)$residuals
#################### infant Example#########################
#A function for computing Pearson correlation for IxJ tables & Mantel-Haenszel, M2
#pearson correlation for IxJ tables
#table = IxJ table or a matrix
#rscore=vector of row scores
#cscore=vector of column scores
pears.cor=function(table, rscore, cscore)
{
    dim=dim(table)
    rbar=sum(margin.table(table,1)*rscore)/sum(table)
    rdif=rscore-rbar
    cbar=sum(margin.table(table,2)*cscore)/sum(table)
    cdif=cscore-cbar
    ssr=sum(margin.table(table,1)*(rdif^2))
    ssc=sum(margin.table(table,2)*(cdif^2))
    ssrc=sum(t(table*rdif)*cdif)
    pcor=ssrc/(sqrt(ssr*ssc))
    pcor
    M2=(sum(table)-1)*pcor^2
    M2
    result=c(pcor, M2)
    result
}


rscore=c(0,0.5,1.5,4,7)
cscore=c(0,1)
x<-factor(rep(c("0","1","2","3","4"), c(17114,14502,793,127,38), levels=c("0","<1","1-2","3-5",">=6")))
y<-factor(rep(c(0,1,0,1,0,1,0,1,0,1),c(17066,48,14464,38,788,5,126,1,37,1),lables=c("not Have","Have")))
table<-table(x,y)
pears.cor(table,rscore,cscore)



#################### Tea Example#########################
x<-factor(rep(c(1,2),c(4,4),lables=c("Milk","Tea"))) 
y<-factor(rep(c(0,1,0,1),c(3,1,1,3),lables=c("G.Milk","G.Tea"))) 
table(x,y) 
fisher.test(x,y,alternative="greater") 
##############clinic example#####################
x=factor(rep(c(1,2,1,2,1,2),c(10,10,7,6,8,9),labels=c("Drug","Placebo"))) 
y=factor(rep(c(0,1,0,1,0,1,0,1,0,1,0,1),c(4,6,8,2,3,4,5,1,3,5,6,3),labels=c("Yes","No")))  
z=factor(rep(c(1,2,3),c(20,13,17),labels=c("clinic1","clinic2","clinic3")))
table(x,y,z)
mantelhaen.test(x,y,z)
#################China Example#######################
cities<-c("Beijing","Shanghai","Shenyang","Nanjing","Harbin","Zhengzhou","Taiyuan","Nanchang")
City<-factor(rep(cities,rep(4,length(cities))),levels=cities)
Smoker<-factor(rep(rep(c("Yes","No"),c(2,2)),8),levels=c("Yes","No"))
Cancer<-factor(rep(c("Yes","No"),16),levels=c("Yes","No"))
count<-c(126,100,35,61,908,688,497,807,913,747,336,598,235,172,58,
121,402,308,121,215,182,156,72,98,60,99,11,43,104,89,21,36)
chismoke<-data.frame(City,Smoker,Cancer,count)
chismoke
x<-tapply(count,list(Smoker,Cancer,City),c)
names(dimnames(x))<-c("Smoker","Cancer","City")
ftable(x,row.vars=c("City","Smoker"),col.vars="Cancer")
mantelhaen.test(x)
#####################Snoring Example #####################
snoring <-matrix(c(24,1355,35,603,21,192,30,224), ncol=2, byrow=TRUE)
dimnames(snoring) <-list(snore=c("never","sometimes","often","always"),heartdisease=c("yes","No"))
scores.a <- c(0,2,4,5) 
scores.b <- c(0,2,4,6) 
scores.c <- 0:3 
scores.d <- 1:4 
snoring.lg.a <- glm( snoring ~ scores.a, family=binomial() ) 
snoring.lg.b <- glm( snoring ~ scores.b, family=binomial() ) 
snoring.lg.c <- glm( snoring ~ scores.c, family=binomial() ) 
snoring.lg.d <- glm( snoring ~ scores.d, family=binomial ()) 
coef(snoring.lg.a)
coef(snoring.lg.b)
coef(snoring.lg.c )
coef(snoring.lg.d) 
predict(snoring.lg.a, type="response" )
predict(snoring.lg.b, type="response" )
predict(snoring.lg.c, type="response" )
predict(snoring.lg.d, type="response" )
snoring.logit <-glm( snoring ~ scores.a, family=binomial(link="logit"))
summary(snoring.logit)
predict(snoring.logit, type="response")
######################Mice Exaple##################
data=matrix(scan(),byrow=T,ncol=4)
1.5 96 44 120
1.5 168 37 80
1.5 336 43 80
1.5 504 35 60
3.5 0.5 29 100
3.5 1 53 200
3.5 2 13 40 
3.5 3 75 200
3.5 5 23 40
3.5 7 152 280
3.5 14 55 80
3.5 24 98 140
3.5 48 121 160 
dim(data)
colnames(data)<-c("no2","time","y","n")
data<-data.frame(data)
y=data$y
n=data$n
no2=data$no2
time=data$time
logno2=log(no2) 
logtime=log(time) 
logno2=(logno2-mean(logno2))/sqrt(var(logno2))
logtime=(logtime-mean(logtime))/sqrt(var(logtime)) 
combdat=cbind(y,nminusy=n-y)
model1=glm(combdat~logno2+logtime,family=binomial(link=logit),data=data)
print(summary(model1))
residuals(model1) 
fitted.values(model1) 
predict.glm(model1) 
anova(model1)
model2=glm(combdat~logno2*logtime,family=binomial(link=logit),data=data) 
print(summary(model2))