########primer ejemplo
logit<-function(x){return(1/(1+exp(-x)))}


x1<-rbinom(300,1,0.5)
y<-rep(NA,300);
for (i in 1:300)
  y[i]<-rbinom(1,1,logit(2*x1[i]))
fit <- glm(y ~ x1, family=binomial(link=logit))
table(x1,y)
x<-c(0,1)
a<-cbind(c(82,129),c(73,16))
fit <- glm(a ~ x-1, family=binomial(link=logit))
summary(fit)





#simulamos datos
x1<-(1:60)/60 -0.5; x2<-rnorm(60,0,1); x3<-rnorm(60,0,1);
y<-rep(0,60);
logit<-function(x){return(1/(1+exp(-x)))}

for (i in 1:60)
  y[i]<-rbinom(1,1,logit(10*x1[i]))

#ajustamos algunos modelos

fit <- glm(y ~ x1+x2+x3, family=binomial(link=logit))

fit <- glm(y ~ x1+x3-1, family=binomial(link=logit))

fit <- glm(y ~ x1-1, family=binomial(link=logit))

plot(x1, fitted(fit),type="l",xlim=c(-1,1))

rug(x1[y==1],col="green",side=3);rug(x1[y==0],col="red")



#caracteristicas importantes en el output son "residual deviance",
#las estimaciones y intervalos de confianza de los parametros 
#y -si no hay demasiado incertidumbre en el modelo - el numero 
#de casos bien clasificados (lo mejor ser basarlo en un conjunto de  observaciones 
#no usadas en la estimacion de los parametros)

########segundo ejemplo

d<-read.table("c:\\datamat.txt"); d1<-as.data.frame(d)

a<-glm(cbind(V1,V2)~V3+V4+V5+V6+V7,family=binomial(logit),data=d1);


a<-glm(cbind(V1,V2)~V3*V4+V3*V5+V3*V6+V3*V7+V4*V5+V4*V6+V4*V7+V5*V6+V5*V7+V6*V7,family=binomial(logit),data=d1);

#aqui se leen las observaciones como tabla con frecuencias; 
#atraves de la funcion "table()" se puede armar esta tabla a partir de las 
#observaciones individuales


# como los predictores toman un conjunto finito de valores, en casoque 
# en las celdas hay suficiente observaciones, se pueden comparar las frecuencias 
# observadas vs estimadas y la normalidad de los residuos 
# a traves del siguiente comando

qqnorm(a$residuals)

######### 
# el siguiente codigo es para enumerar todos los posibles modelos con terminos 
# que aparecen en "alfabeto"; es MUY lenta si alfabeto es grande!!!
# se arranca con itera() y el resultado esta en "kans2" 
#########

alfabet2<-c("V3","V4","V5","V6", "V7","V3*V4",
"V3*V5","V3*V6","V3*V7","V4*V5","V4*V6","V4*V7","V5*V6","V5*V7","V6*V7")

alfabet<-c("V3","V4","V5","V6", "V7")

ultimo<-length(alfabet) 
n<-ultimo+1 
model<<-rep(0,n)

model[1]<--999999999999 ;kans2<<-NULL

#la siguiente funcion genera todas las combinaciones

itera <- function() { work<-TRUE;
 while(work)
 {if (model[n]< ultimo)
   {model[n]<<-model[n]+1;
    procesa(model[-1]);
   }
  else
  {j<-n;
   while(model[j]>=ultimo-(n-j))
       {j <-j-1; }
   if (j >1)
    {model[j]<<-model[j]+1;
     for (k in (j+1):n) model[k]<<-model[j]+(k-j)
     procesa(model[-1]) ;
    }
    else work<-FALSE

   }

 }
}

procesa <-function(m) {print(m); s<-" ";
 i<-1;
 while (m[i]==0)
   {i<-i+1}
 for (k in i:(n-1))
  {
   if (m[k]>0) {s<-paste(s, alfabet[m[k]])
                 if (k<(n-1)) s<-paste(s, "+");
                }
  }

 print(s)
#observa como en el siguiente formula, se convierte un string en una expresion
a<-glm(as.formula(paste("cbind(V1,V2)~",s)),family=binomial(logit),data=d1);
kans2<<-rbind(kans2,c(a$df.residual,a$deviance)); }

###########  START ###################
itera()
plot(kans2,ylim=c(0,60),xlab="d.f",ylab="residual deviance")
points(16:30,qchisq(0.95,16:30),col="red",type="l")

########### 
########  tercer ejemplo 

crabs <- read.table("crab.dat",header=T) ;
attach(crabs);
 y <-ifelse(satellites>0, 1, 0); 
 plot(width,jitter(y,0.2))