sink('c:/MLG/ss.sai')
date()
diamond.dat<-scan(what=list(res=0,cos=0,tot=0,y=0))
1 1  958 27
1 2  160  1
1 3   65  3
1 4  293 12
2 1 2670 67
2 2  338 11
2 3  237 11
2 4  654 23
3 1  543  7
3 2   70  4
3 3  110  3
3 4  260  7

attach(diamond.dat)
res<-diamond.dat$res
res<-factor(res)
res<-C(res,treatment)
cos<-diamond.dat$cos
cos<-factor(cos)
cos<-C(cos,treatment)
y<-diamond.dat$y
n<-diamond.dat$tot
Xmat<-cbind(y,n-y)
print("Modelo 1")
diamond1.fit<-glm(Xmat~cos,family=binomial())
summary(diamond1.fit)
diamond2.fit<-glm(Xmat~res,family=binomial())
summary(diamond2.fit)
diamond3.fit<-glm(Xmat~cos+res,family=binomial())
summary(diamond3.fit)
diamond4.fit<-glm(Xmat~cos+res+cos*res,family=binomial())
summary(diamond4.fit)
anova(diamond1.fit,diamond3.fit,diamond4.fit,test="Chi")
anova(diamond2.fit,diamond3.fit,diamond4.fit,test="Chi")
##########teste da RV
 anova(diamond1.fit,diamond3.fit,test="Chi")
 RV<-deviance(diamond1.fit)-deviance(diamond3.fit)
pv<-1-pchisq(RV,2)
#teste F
 anova(diamond1.fit,diamond3.fit,test="F")
#############teste escore
##modelo aditivo
X<-model.matrix(diamond3.fit)
X1<-X[,c(5,6)]
X2<-X[,c(1,2,3,4)]
rp<-resid(diamond1.fit,type="pearson")
w<-diamond1.fit$weights
W<-diag(w)
A<-solve(t(X2)%*%W%*%X2)
Ci<-A%*%t(X2)%*%W%*%X1
Ri<-X1-X2%*%Ci
varo<-solve(t(Ri)%*%W%*%Ri)
SR<-t(rp)%*%sqrt(W)%*%X1%*%varo%*%t(X1)%*%sqrt(W)%*%rp
pv<-1-pchisq(SR,ncol(X1))
#######################################################
########Teste wald
phi<-1
betah<-coef(diamond3.fit)[c(5,6)]
betah<-matrix(betah,2,1)
w<-diamond3.fit$weights
W<-diag(w)
A<-solve(t(X2)%*%W%*%X2)
Ci<-A%*%t(X2)%*%W%*%X1
Ri<-X1-X2%*%Ci
var1<-t(Ri)%*%W%*%Ri
Wald<-phi*t(betah)%*%var1%*%betah
pv<-1-pchisq(Wald,2)
######teste  da razao de verossimilhança

 #########dados agrupados

 diamondg.dat<-scan(what=list(res=0,cos=0,tot=0,y=0))
1 1  958 27
1 1  160  1
1 3   65  3
1 4  293 12
2 1 2670 67
2 1  338 11
2 3  237 11
2 4  654 23
3 1  543  7
3 1   70  4
3 3  110  3
3 4  260  7

attach(diamondg.dat)
res<-diamondg.dat$res
res<-factor(res)
res<-C(res,treatment)
cos<-diamondg.dat$cos
cos<-factor(cos)
cos<-C(cos,treatment)
y<-diamondg.dat$y
n<-diamondg.dat$tot
Xmat<-cbind(y,n-y)
print("Modelo 1")
diamond1.fit<-glm(Xmat~cos,family=binomial())
summary(diamond1.fit)
diamond2.fit<-glm(Xmat~res,family=binomial())
summary(diamond2.fit)
diamond3.fit<-glm(Xmat~cos+res,family=binomial())
summary(diamond3.fit)
diamond4.fit<-glm(Xmat~cos+res+cos*res,family=binomial())
summary(diamond4.fit)
anova(diamond1.fit,diamond3.fit,diamond4.fit,test="Chi")
anova(diamond2.fit,diamond3.fit,diamond4.fit,test="Chi")
#####################
###################Modelo final so com a variavel grau de parentesco
summary(diamond1.fit)