R/simpmod.r

In inecolr: Modeling and Plotting for Ecologist

simpmod <- function(M, alpha=0.05){
  NAMES <- character()
  k<-0
  History <- list()
  Msimp <- M
  BETAS_NAMES <- rownames(anova(Msimp))
  BETAS_NAMES <- BETAS_NAMES[BETAS_NAMES != "NULL"]
  BETAS_NAMES <- BETAS_NAMES[BETAS_NAMES != "(Intercept)"]
  BETAS_NAMES <- BETAS_NAMES[BETAS_NAMES != "Residuals"]
  BETAS_NAMES_LIST <- str_split(BETAS_NAMES, ":")
  NIVEL_INTER <- unlist(lapply(BETAS_NAMES_LIST, length))
  for(i in sort(unique(NIVEL_INTER), decreasing = T)){
    INT_NIV_X <- which(NIVEL_INTER == i)
    TO_DROP <- numeric()
    for(j in INT_NIV_X[length(INT_NIV_X):1]){
      BETAS_NAMES_LISTA <-  str_split(BETAS_NAMES[-j], ":")
      FOCAL_BETA_NAME <- str_split(BETAS_NAMES[j], ":")
      IN_OTHER_INT <- numeric(length(BETAS_NAMES_LISTA))
      for(w in unlist(FOCAL_BETA_NAME)) IN_OTHER_INT <- IN_OTHER_INT + unlist(lapply(lapply(BETAS_NAMES_LISTA, str_detect, w), sum))
      
      if(sum(IN_OTHER_INT == i) < 1) {
        Msimp2 <- eval(parse(text = paste("update(Msimp,", "~. -", BETAS_NAMES[j], ")",sep="")))
        if(class(Msimp)[1]=="lm") {
          TEST <- anova(Msimp2, Msimp)
          if(TEST[2,6] > alpha ) {
            k <- k + 1
            History[[k]] <- TEST
            NAMES <- c(NAMES, BETAS_NAMES[j])
            TO_DROP <- paste(TO_DROP, paste("-", BETAS_NAMES[j],sep=""))
            BETAS_NAMES <- BETAS_NAMES[-j]
          }
        }
        if(class(Msimp)[1]=="glm") {
          TEST <- anova(Msimp2, Msimp, test="Chi")
          if(TEST[2,5] > alpha ) {
            k <- k + 1
            History[[k]] <- TEST
            NAMES <- c(NAMES, BETAS_NAMES[j])
            TO_DROP <- paste(TO_DROP, paste("-", BETAS_NAMES[j],sep=""))
            BETAS_NAMES <- BETAS_NAMES[-j]
          }
        }
        
        if(class(Msimp)[1]=="negbin") {
          TEST <- anova(Msimp2, Msimp)
          if(TEST[2,7] > alpha ) {
            k <- k + 1
            History[[k]] <- TEST
            NAMES <- c(NAMES, BETAS_NAMES[j])
            TO_DROP <- paste(TO_DROP, paste("-", BETAS_NAMES[j],sep=""))
            BETAS_NAMES <- BETAS_NAMES[-j]
          }
        }
        
      }
    }
    if(length(TO_DROP) != 0) Msimp <- eval(parse(text = paste("update(Msimp, ~.", TO_DROP, ")", sep=""))) 
  }
  if(k>0) {History[[k+1]] <- Msimp
  names(History) <- c(paste("-",NAMES,sep=""), "Simplified_Model")
  History
  } else History <- list(Simplified_Model =Msimp )
  History
}