#####################
#Regresion Logistica#
#####################

rm(list=ls())
library(ISLR)
attach(Default)
data=Default
data
head(data,n=4)

modelo=glm(default~balance,family=binomial,data)
summary(modelo)
anova(modelo,test='Chisq')

1-pchisq(modelo$deviance,9998)
#no rechazamos el modelo

anova(modelo,test='Chisq')
a=modelo$deviance

predict(modelo, data.frame(balance=c(1000,2000)), type='response')

modelo2=glm(default~student,family=binomial,data)
summary(modelo2)
predict(modelo2, data.frame(student=c("Yes","No")), type='response')


modelo3=glm(default~.,family=binomial,data)

summary(modelo3)
predict(modelo3, data.frame(student="Yes",balance=1500,income=40000), type='response')

1-pchisq(modelo3$deviance,9996)
#no rechazamos el modelo

anova(modelo,modelo3,test='Chisq')
#Hay evidencia que el modelo3 es mejor que modelo
AIC(modelo)-AIC(modelo3)
#Idem

pred=predict(modelo3,data.frame=data,type="response")
prediccion=rep("No",10000)
prediccion[pred>.5]="Yes"
tabla=table(prediccion,default)
mean(prediccion==default)

#Risk ratio RR (P(Y=1|x)/P(Y=1|x) siendo x y x'que difieren solo en la variable estudiante)

p1=predict(modelo3,data.frame(student="No",balance=1500,income=40000),type="response")
p2=predict(modelo3, data.frame(student="Yes",balance=1500,income=40000), type='response')

p1/p2
#el riesgo de entrar en default para un mismo valor de balance y de ingreso si uno no es estudiante
#de si uno es estudiante es casi el doble


##############
#OTRO EJEMPLO#
##############
#
rm(list=ls())
#Datos AGE_CHD: Presencia o ausencia de enfermedad coronaria con la edad
#Variable X (AGE): cuantitativa
#variable respuesta Y (CHD): 1 si presenta enfermedad, 0 si no presenta enfermedad
#levantamos datos
data=read.table("AGE_CHD.csv",dec=",",sep=";",header=T)
summary(data)
attach(data)
data$CHD=as.factor(data$CHD)
summary(data)


modelo=glm(CHD~AGE,family=binomial,data)
summary(modelo)

#beta_0 y beta_1 resultan ser significativos.
#El modelo es ln(p/(1-p))=-5.30945+0.11092X lo cual tambien 
#se puede escribir como p=1/(1+e^{-5.30945+0.11092X})

#Interpretacion coeficiente
#beta_1=0.11092=ln(OR) lo cual significa que OR=1.1173 e indica como aumenta
#el riesgo de padecer enfermedad coronaria al aumentar la edad de un año.
#Si la edad aumenta 20 años: e^{20*beta_1)=9,2...


#Una predicci?n:
pp=predict(modelo,data.frame(AGE=36),type="response")

#plot
plot(AGE, fitted(modelo), xlab="edad", ylab="ajustados", col="red")

#bondad de ajuste mediante la desvianza
anova(modelo, test='Chisq')

#H0: modelo = modelo nulo
#H1: rechazo H0
#p-valor sigue rechazo la hipótesis nula de que nuestro modelo 
#se asemeje a un modelo tan simple como el modelo nulo. 



#REGRESIÓN LOGÍSTICA MÚLTIPLE.
#Datos Tibet

##Esos datos correponden a dos tipos raciales diferentes en los que se
#practicaron diferentes medidas de longitud y ancho de craneo y cara.


rm(list=ls())
source("datos tibet.dat")
summary(Tibet)
attach(Tibet)
names(Tibet)

#modelo1

modelo1=glm(Type~Length+Breadth+Height,family=binomial,Tibet)
summary(modelo1)
anova(modelo1, test='Chisq')

modelo2=glm(Type~Length,family=binomial,Tibet)
summary(modelo2)

#Comparo modelos

anova(modelo1,modelo2,test='Chisq')

#NO TENEMOS EVIDENCIA PARA CONCLUIR QUE EL MODELO CON MÁS PARÁMETROS (MODELO1)
#SEA MEJOR QUE EL MODELO CON MENOS PARÁMETROS (MODELO2) ASÍ QUE NOS 
#QUEDAMOS CON EL MODELO MÁS SIMPLE.
#AL COMPARAR LOS AIC LLEGAMOS A LA MISMA CONCLUSIÓN:

AIC(modelo1)
AIC(modelo2)

rm( list=ls() )

source('datos tibet.dat')
str( Tibet )

## Ejercicio 2 y 3 (se usan datos de Tibet para el 2, los datos simulados quedan de ejercicio).

##1

## Esta es una funcion para fijar la funcion de de verosimilitud a maximizar, tiene como argumentos los datos de x e y,
#los cuales se fijan dentro de otra funcion de argumento beta que es la verosimilitud
la.L <- function(y, x) function(beta)
  t( y ) %*% ( x %*% beta ) - t( rep(1, length(y)) ) %*% log( 1 + exp( x %*% beta ) )

## Se determinan la verosimilitud con los datos de Type como respuesta y Length como varaible independiente.
L <- la.L(as.numeric(Tibet$Type)-1, cbind(1, Tibet$Length) )
L(c(1, 0))

## Para maximizar se usa el argumento 'control'. Otra forma es no usarlo y en la definicion de la funcion 
#a maximizar poner el resultado por -1, esto hace que el problema sea el de hallar un minimo que es el objetivo por defecto de la funcion optim.
optim( c(1, -50), L, method='BFGS', control = list(fnscale= -1) )$par
optim( c(1050, -56440), L, method='BFGS', control = list(fnscale= -1) )$par

## Podemos poner lo anterior dentro de una funcion que tome como argumento la L.
estimacion <- function(L, valor.inicial=c(1, 0), metodo='BFGS')
  optim( valor.inicial, L, method=metodo, control=list( fnscale = -1 ) )$par

## Ahora el trabajo se reduce a crear la verosimilitud con la.L() y pasarla como argumento a estimacion()
L <- la.L(as.numeric(Tibet$Type)-1, cbind(1, Tibet$Length) )
estimacion( L )

## Se puede poner directo la funcion la.L()
estimacion( la.L(as.numeric(Tibet$Type)-1, cbind(1, Tibet$Length) ) )

## Tambien se puede definir dentro de estimacion() la funcion que crea la verosimilitud.

## Resultados con glm.
glm( Tibet$Type ~ Tibet$Length, family=binomial)

modelo=glm( Tibet$Type ~ Tibet$Length, family=binomial)
attributes( modelo)
y.techo <- modelo$fitted.values
plot( exp( y.techo ) / (1 + exp( y.techo ) ), ylab='p estimado', ylim=c(0, 1) )

predict(modelo,data=Tibet)

## Ejercicio 4

rm( list=ls() )

if ( !require( nnet ) )
{
  chooseCRANmirror( graphics=F )
  install.packages( 'nnet' )
  library( nnet )
}

( m <- multinom( Species ~ ., data = iris ) )
summary( m )
predict(m, data=iris, type='class')
predict(m, data=iris, type='probs')

## Se extrae una muestra de 25 elementos por cada especie, se predice y se compara con los valores reales.
samp <- c(sample(1:50, 25), sample(51:100, 25), sample(101:150, 25))
table( iris$Species[ -samp], predict(m, iris[ -samp,], type='class'))
mean( iris$Species[ -samp ] == predict(m, iris[ -samp, ], type='class') )

## Ahora comparamos con todos.
table( iris$Species, predict(m, iris, type='class') )
mean( iris$Species == predict(m, iris, type='class') )