GLM_fit = function(pmm, X, family) {
  
  if (family == "gamma") {
    links = c("log", "inverse", "identity")   
#    links = c("log", "inverse")   
    
    # clean data and remove zeros
    precip_dt = na.omit(precip_dt)
    rownames(precip_dt) = NULL
    precip_dt$pmm = ifelse(precip_dt$pmm <=0, runif(1, 0.0001, 0.001), precip_dt$pmm)
    
  } else if (family == "binomial"){
    links = c("logit", "probit", "cauchit")
#    links = c("logit", "probit")

  }
  
  glm_aic = vector(length = length(links))
  N = length(pmm)
  
  combs = leaps(X,pmm, nbest=25)     #  GEt upto 25 combinations for each
  # number of predictors
  combos = combs$which
  
  ncombos = length(combos[,1])
  
  xpress=1:ncombos
  xmse = 1:ncombos
  
  for(j in 1:length(links)) {
    
    for(i in 1:ncombos) {
      
      xx = X[,combos[i,]]
      xx=as.data.frame(xx)
      if (family == "gamma") {
         zz=glm(pmm ~ ., data=xx, family = Gamma(link=links[j]), maxit=500)
      } else if (family == "binomial"){
         zz=glm(pmm ~ ., data=xx, family = binomial(link=links[j]), maxit=500)
      }
      
      xpress[i]=PRESS(zz)
      
      xmse[i] = sum((zz$res)^2) / (N - length(zz$coef))
    }
    
    # Test using AIC objective function
    zms=stepAIC(zz, scope=list(upper = ~., lower = ~1), trace=FALSE)
    glm_aic[j] = zms$aic
  } 
  
  aic_df = data.frame(glm_aic)
  rownames(aic_df) = links[1:3]
  
  print("Results of AIC for bestfit GLM")
  print(aic_df)
  
  sprintf("Choosing the GLM which minimizes AIC: %s family and %s link function.", family, links[which.min(glm_aic)])
  # bestmod = glm(pmm ~ ., data = X, family = Gamma(link=links[which.min(glm_aic)]))
  
  if (family == "gamma") {
    bestmod = glm(pmm ~ ., data = X, family = Gamma(link=links[which.min(glm_aic)]))
  } else if (family == "binomial") {
    bestmod = glm(pmm ~ ., data = X, family = binomial(link=links[which.min(glm_aic)]))
  } else { 
    print("Error!")
  }
  
  return(bestmod)
}