vars=5:12
dependientes=datos[,vars]
independientes=datos[,-vars]

fits=
  lapply(                                                       #Aplico a cada nombre de las variables dependientes
    colnames(dependientes),
      function(x)                                               #Esta función donde:
        step(                                                   #Se usa la función step para seleccionar el mejor modelo AIC quitando variables
          glm(                                                  #El modelo glm de comienzo, donde se usa todas las variables d
            reformulate(                                        #Reformulate es para crear formula, donde la respuesta es el caracter x (la variable de la función, donde será reemplazado por cada nombre de colnames(dependientes)
              termlabels = colnames(independientes),            #La variables independientes
              response =  x,                                    #La variable dependientes
              intercept =T
            ),
            data = cbind(dependientes,independientes),        #Los datos
            family = binomial("logit")                        #El modelo Logit
          )
        )
      )

#Modelos Probits
fits2=
lapply(
  colnames(dependientes),
  function(x)
    step(
      glm(
        reformulate(
          termlabels = colnames(independientes),
          response =  x,
          intercept =T
        ),
        data = cbind(dependientes,independientes),
        family = binomial("probit")
      )
    )
  )
fits=c(fits,fits2) #Concatenar los modelos

fits3=
  do.call(
    rbind,
    lapply(
      c(1:22),
      function(i){
        if(dim(summary(fits[[i]])$coefficients)[1]!=0) {
          data.frame(
            dep=names(fits[[i]]$coefficients),
            summary(fits[[i]])$coefficients,
            respond=as.character(fits[[i]]$formula)[2],
            modelo=i,
            AIC=AIC(fits[[i]]),
            BIC=BIC(fits[[i]]),
            Desvianza=fits[[i]]$deviance,
            R=1-fits[[i]]$df.residual/fits[[i]]$df.null
          )
        } else {
         NULL
        }
      }
    )
  )