logistic.R

datapp <- read.csv("data_pib.csv",header=T)
datap2 <- datapp[datapp$sex<99,]
datap <- datap2[datap2$skin<99,]
attach(datap)
plot(age,skin)
plot(length,skin)

tumor.table <- table(skin,age)
tumor.table
ages <- as.real(colnames(tumor.table))
p.tumor <- tumor.table[2,]/(tumor.table[1,]+tumor.table[2,])
p.tumor
plot(ages,p.tumor,ylim=c(0,1),pch=20,xlab="Age",ylab="Probability of Skin Tumor")
no.age <- tumor.table[1,]+tumor.table[2,]
plot(ages,p.tumor,ylim=c(0,1),pch=20,xlab="Age",ylab="Probability of Skin Tumor",cex=no.age)

library(ggplot2)
p <- qplot(ages,p.tumor,size=no.age,ylim=c(0,1),xlim=c(1,16))
p <- p+geom_point(color="red") + scale_area(to=c(1,15))
p + xlab("Age of Bullhead") + ylab("Probability of Skin Tumor")

res <- glm(skin~age,family=binomial)
summary(res)
plot(res,1)

nd <- data.frame(age=seq(1,16))
predict(res,newdata=nd)
predict(res,newdata=nd,type="response")

plot(ages,p.tumor,ylim=c(0,1),pch=20,xlab="Age",ylab="Probability of Skin Tumor",cex=no.age/4)
matlines(nd,predict(res,newdata=nd,type="response"),lwd=2)

pred.res <- predict(res,newdata=nd,se.fit=TRUE)
lower.bound <- pred.res$fit - qt(.975,14)*pred.res$se.fit
upper.bound <- pred.res$fit + qt(.975,14)*pred.res$se.fit
p.lower.bound <- exp(lower.bound)/(1+exp(lower.bound))
p.upper.bound <- exp(upper.bound)/(1+exp(upper.bound))
 
plot(ages,p.tumor,ylim=c(0,1),pch=20,xlab="Age",ylab="Probability of Skin Tumor",cex=no.age/4)
matlines(nd,predict(res,newdata=nd,type="response"),lwd=2)
matlines(nd,p.lower.bound,lwd=2,lty=2,col="red")
matlines(nd,p.upper.bound,lwd=2,lty=2,col="red")

res <- glm(skin~age+tsex,family=binomial)
summary(res)

library(MASS)
res <- stepAIC(glm(skin~(age+length+weight+tsex)^2,family=binomial))
summary(res)

res <- stepAIC(glm(skin~(age+length+weight+tsex)^2+pibloc,family=binomial))
summary(res)

res1 <- glm(skin ~ age + length + tsex + age:length + age:tsex,family=binomial)
res2 <- glm(skin ~ age + length + tsex + age:length + age:tsex + pibloc,family=binomial)

summary(res1)
summary(res2)

anova(res1,res2)
AIC(res1)
AIC(res2)