#######################################################
#
# Regressao logistica
#
# Nesse exemplo o interesse esta' em modelar a variavel
# dicotomica result (1=better,0=not better) atraves de 3
# variaveis explicativas: tipo de tratamento, sexo do
# individuo e idade do individuo.
#
# column 1: treatment (0=placebo)
# column 2: sex (1=male)
# column 3: age
# column 4: result(1=better)
#######################################################

data = matrix(scan("arthritis.txt"),byrow=T,ncol=4)

x1 = data[,1]
x2 = data[,2]
x3 = data[,3]
y   = data[,4]


#######################
# Modelo 1: Por tratamento
#######################
mylogit = glm(y~x1, family = "binomial")
summary(mylogit)

# P(better | placebo)
num=sum(exp(mylogit$coef*c(1,0)))
num/(1+num)

# P(better | tratamento)
num=sum(exp(mylogit$coef*c(1,1)))
num/(1+num)

#######################
# Modelo 2: Por sexo
#######################
mylogit = glm(y~x2, family = "binomial")
summary(mylogit)

# P(better | mulher)
num=sum(exp(mylogit$coef*c(1,0)))
num/(1+num)

# P(better | homen)
num=sum(exp(mylogit$coef*c(1,1)))
num/(1+num)

#######################
# Modelo 3: Por idade
#######################
mylogit = glm(y~x3, family = "binomial")
summary(mylogit)

# P(better | idade)
plot(x3,y,xlab="Idade (x3)",ylab="Prob(Better|x3)",pch=16)
points(x3,mylogit$fitted.values,col=2,pch=16)

#################################
# Modelo 4: Por tratamento e sexo
#################################
mylogit = glm(y~x1+x2, family = "binomial")
summary(mylogit)

# P(better | placebo,mulher)
num=sum(exp(mylogit$coef*c(1,0,0)))
num/(1+num)

# P(better | placebo,homem)
num=sum(exp(mylogit$coef*c(1,0,1)))
num/(1+num)

# P(better | tratamento,mulher)
num=sum(exp(mylogit$coef*c(1,1,0)))
num/(1+num)

# P(better | tratamento,homen)
num=sum(exp(mylogit$coef*c(1,1,1)))
num/(1+num)

####################################
# Modelo 5: Por tratamento, sexo e idade
####################################
mylogit = glm(y~x1+x2+x3, family = "binomial")
summary(mylogit)

xx3 = seq(min(x3),max(x3),length=100)
# P(better | placebo,mulher,idade)
fit1 = cbind(1,0,0,xx3)%*%mylogit$coef
fit1 = exp(fit1)/(1+exp(fit1))
# P(better | placebo,homem,idade)
fit2 = cbind(1,0,1,xx3)%*%mylogit$coef
fit2 = exp(fit2)/(1+exp(fit2))
# P(better | tratamento,mulher,idade)
fit3 = cbind(1,1,0,xx3)%*%mylogit$coef
fit3 = exp(fit3)/(1+exp(fit3))
# P(better | tratamento,homem,idade)
fit4 = cbind(1,1,1,xx3)%*%mylogit$coef
fit4 = exp(fit4)/(1+exp(fit4))

plot(x3,y,xlab="Idade (x3)",ylab="Prob(Better|x1,x2,x3)",pch=16)
lines(xx3,fit1,col=2,lwd=3,lty=2)
lines(xx3,fit2,col=3,lwd=3,lty=2)
lines(xx3,fit3,col=2,lwd=3)
lines(xx3,fit4,col=3,lwd=3)
legend(25,0.95,legend=c("Tratamento,Mulher","Tratamento,Homen","Placebo,Mulher","Placebo,Homen"),col=c(2,3,2,3),lty=c(1,1,2,2),lwd=3,bty="n")