source("ghq.dat") attach(ghq) # p<-case/(case+nocase) ghq<-data.frame(ghq,p) attach(ghq) # fit.lin<-glm(p~ghqscore) fit.log<-glm(p~ghqscore,family=binomial,weights=case+nocase) summary(fit.lin) summary(fit.log) # # # lin.pred<-predict(fit.lin) log.pred<-exp(predict(fit.log))/(1+exp(predict(fit.log))) cbind(lin.pred,log.pred) # Case<-case[1:11]+case[12:22] Nocase<-(nocase[1:11]+nocase[12:22]) P<-Case/(Case+Nocase) # ylim<-range(log.pred,lin.pred,P) win.graph() par(pty="s") plot(ghqscore[1:11],lin.pred[1:11],xlab="GHQ",ylab="Probability of case",type="l",ylim=ylim) lines(ghqscore[1:11],log.pred[1:11],lty=2) points(ghqscore[1:11],P) legend(0.3,1,c("Linear","Logistic"),lty=1:2) # # fit1<-glm(p~ghqscore,family=binomial,weights=case+nocase,data=ghq) fit2<-glm(p~ghqscore+sex,family=binomial,weights=case+nocase,data=ghq) fit3<-glm(p~ghqscore*sex,family=binomial,weights=case+nocase,data=ghq) anova(fit1,fit2,fit3) # summary(fit2) # resd<-residuals(fit2,type="deviance") resp<-residuals(fit2,type="pearson") # resds<-resd/sqrt(1-lm.influence(fit2)$hat) resps<-resp/sqrt(1-lm.influence(fit2)$hat) ghq.pred<-predict(fit2) # par(mfrow=c(1,3)) plot(1:22,resds,xlab="Observation number",ylab="Standardized deviance residual") plot(1:22,resps,xlab="Observation number",ylab="Standardized Pearson residual") # plot(ghq.pred,resds,xlab="Predicted value",ylab="Standardized deviance residual") plot(ghq.pred,resps,xlab="Predicted value",ylab="Standardized Pearson residual") # plot(ghqscore,resds,xlab="GHQ",ylab="Standardized deviance residual")