R/ranfor.R

In inecolr: Modeling and Plotting for Ecologist

#BOSQUE DE ARBOLES----
ranfor <- function(data, repetitions=1000, prop=0.85, n.digits = 5, plot=TRUE){
  userpars <-par(no.readonly = TRUE) #collect users parameters
  on.exit(par(userpars)) #reset users parameters on exit of the
  if(is.numeric(data[,1]) | is.factor(data[,1])){
  RESULTADOS <- matrix(0, 1000, ncol(data))
  colnames(RESULTADOS) <- c(names(data)[2 : ncol(data)], "<leaf>")
  
  NO.NODO <- integer()
  NO.NODO.VAR <- character()
  NO.NODO.SPLIT <- character()
  NO.NODO.Y <- numeric()
  NO.NODO.DEV <- numeric()
  NO.NODO.N <- numeric()
  
  for(k in 1:repetitions){
    M <- sample(1:nrow(data), size=round(nrow(data)*prop))
    data.2 <- data[M, ]
    modelo <- tree::tree(data.2)
    NODOS <- as.numeric((rownames(modelo$frame)))
    POS <- which(NODOS <= 1000)
    NODOS <- NODOS[POS]
    NO.NODO <- c(NO.NODO, NODOS)
    NO.NODO.VAR <- c(NO.NODO.VAR, as.character(modelo$frame[POS, 1]))
    NO.NODO.SPLIT <- c(NO.NODO.SPLIT, as.character(modelo$frame[, 5][POS,1]))
    NO.NODO.Y <- c(NO.NODO.Y, modelo$frame[POS, 4])
    NO.NODO.DEV <- c(NO.NODO.DEV, modelo$frame[POS, 3])
    NO.NODO.N <- c(NO.NODO.N, modelo$frame[POS, 2])
    
    POS <- cbind(as.numeric(rownames(modelo$frame)), match(as.character(modelo$frame[,1]), colnames(RESULTADOS)))
    POS<- POS[POS[,1]<= 1000,]
    RESULTADOS[POS] <-   RESULTADOS[POS]+1
  }
  
  #SUMARIZING TREE MODELS
  ANALISIS <- data.frame(NO.NODO, NO.NODO.VAR, NO.NODO.SPLIT, NO.NODO.Y, NO.NODO.DEV, NO.NODO.N)
  #attach(ANALISIS)
  rm(NO.NODO, NO.NODO.VAR, NO.NODO.SPLIT, NO.NODO.Y, NO.NODO.DEV, NO.NODO.N)
  REF <- c(2,4,8,16,32,64,128,256,512,1024,2048,4096)
  
  TABLA <- table(ANALISIS$NO.NODO, ANALISIS$NO.NODO.VAR)
  LEAVES <- rep(0, nrow(TABLA))
  for(i in 1:nrow(TABLA)) {
    if(sum((TABLA[i, 1] < TABLA[i, -1]))==0) LEAVES[i] <- i
  }
  
  TOREMOVE <- numeric()
  for(i in nrow(TABLA):1){
    KKK<- as.numeric(rownames(TABLA)[i])
    ASCENDENTES <- numeric()
    while(KKK>0){
      if(gtools::odd(KKK)){ASCENDENTES <- c(ASCENDENTES, (KKK-1)/2)
      KKK <- (KKK-1)/2
      }
      if(gtools::even(KKK) & KKK!=0){ASCENDENTES <- c(ASCENDENTES, (KKK)/2)
      KKK <- (KKK)/2
      }
    }
    if(sum(LEAVES[match(ASCENDENTES, as.numeric(rownames(TABLA)))], na.rm=T)>0) TOREMOVE<- c(TOREMOVE,i)
  }
  if(length(TOREMOVE) > 0) {
    LEAVES <- LEAVES[-TOREMOVE]
  TABLA.VAR <- TABLA[-TOREMOVE, ]
  TABLA.Y <- tapply(ANALISIS$NO.NODO.Y, list(ANALISIS$NO.NODO, ANALISIS$NO.NODO.VAR), mean, na.rm=T)[-TOREMOVE, ]
  TABLA.Yse <- tapply(ANALISIS$NO.NODO.Y, list(ANALISIS$NO.NODO, ANALISIS$NO.NODO.VAR), sd, na.rm=T)[-TOREMOVE, ]
  TABLA.DEV <- tapply(ANALISIS$NO.NODO.DEV, list(ANALISIS$NO.NODO, ANALISIS$NO.NODO.VAR), mean, na.rm=T)[-TOREMOVE, ]
  TABLA.N <- tapply(ANALISIS$NO.NODO.N, list(ANALISIS$NO.NODO, ANALISIS$NO.NODO.VAR), mean, na.rm=T)[-TOREMOVE, ]
  TABLA.SPLIT <- tapply(as.numeric(substr(ANALISIS$NO.NODO.SPLIT, 2, stringr::str_length(ANALISIS$NO.NODO.SPLIT))), list(ANALISIS$NO.NODO, ANALISIS$NO.NODO.VAR), mean, na.rm=T)[-TOREMOVE, ]
  } else  {TABLA.VAR <- TABLA
  
  TABLA.Y <- tapply(ANALISIS$NO.NODO.Y, list(ANALISIS$NO.NODO, ANALISIS$NO.NODO.VAR), mean, na.rm=T)
  TABLA.Yse <- tapply(ANALISIS$NO.NODO.Y, list(ANALISIS$NO.NODO, ANALISIS$NO.NODO.VAR), se, na.rm=T)
  TABLA.DEV <- tapply(ANALISIS$NO.NODO.DEV, list(ANALISIS$NO.NODO, ANALISIS$NO.NODO.VAR), mean, na.rm=T)
  TABLA.N <- tapply(ANALISIS$NO.NODO.N, list(ANALISIS$NO.NODO, ANALISIS$NO.NODO.VAR), mean, na.rm=T)
  TABLA.SPLIT <- tapply(as.numeric(substr(ANALISIS$NO.NODO.SPLIT, 2, str_length(ANALISIS$NO.NODO.SPLIT))), list(ANALISIS$NO.NODO, ANALISIS$NO.NODO.VAR), mean, na.rm=T)
  }
  
  RESULTADOS2 <- data.frame(matrix(0, nrow(TABLA.VAR), 6))
  colnames(RESULTADOS2) <- c("VAR.DIV", "PUNTO.CORTE", "DEVIANZA", "N", "Y", "SE")
  rownames(RESULTADOS2) <- rownames(TABLA.VAR)
  RESULTADOS2[, 1] <- as.numeric(rownames(RESULTADOS))
  for(i in 1:nrow(TABLA.VAR)){
    ESTEMAX <- which(TABLA.VAR[i,] == max(TABLA.VAR[i,]))[1]
    RESULTADOS2[i, 1] <- colnames(TABLA.VAR)[ESTEMAX]
    RESULTADOS2[i, 2] <-  TABLA.SPLIT[i, ESTEMAX]
    RESULTADOS2[i, 3] <-  TABLA.DEV[i, ESTEMAX]
    RESULTADOS2[i, 4] <-  round(TABLA.N[i, ESTEMAX],0)
    RESULTADOS2[i, 5] <-  TABLA.Y[i, ESTEMAX]
    RESULTADOS2[i, 6] <-  TABLA.Yse[i, ESTEMAX]
  }
  
  RESULTADOS2$Prop.Red.Dev <- 0
  RESULTADOS2$Prop.Red.Y <- 0
  k<-0
 for(i in rownames(RESULTADOS2)){
   k<-k+1
   if(i == 1) RESULTADOS2 [1, 7] <- 0
     nodos_anindados <- (as.numeric(i))*c(2, 2)+c(0,1)
     aqui_nodos_aninados <- match(nodos_anindados, rownames(RESULTADOS2))
     if(any(!is.na(aqui_nodos_aninados))) {
       RESULTADOS2[aqui_nodos_aninados, 7] <- (RESULTADOS2[k, 3] - RESULTADOS2[aqui_nodos_aninados, 3])/RESULTADOS2[k, 3]
     
       RESULTADOS2[aqui_nodos_aninados, 8] <- (RESULTADOS2[k, 5] - RESULTADOS2[aqui_nodos_aninados, 5])/(RESULTADOS2[k, 5] + RESULTADOS2[aqui_nodos_aninados, 5])
     }
 }
  RESULTADOS2
  
  MRD <- mean(RESULTADOS2[which(RESULTADOS2[,1] == "<leaf>"), 3], na.rm = TRUE)
  MED <-  1 - MRD/RESULTADOS2[1,3]
  RESULTADOS3 <- matrix(c(MRD, MED), 2, 1)
  rownames(RESULTADOS3) <- c("Mean_residual_deviance = ","Mean explained deviance")
  colnames(RESULTADOS3) <- c("_")
  
  RESULTADOS2[which(RESULTADOS2[,1] == "<leaf>"), 2] <- NA
  
  M2 <- list(RESULTADOS, RESULTADOS2, RESULTADOS3)

  
  
  if(plot=="TRUE"){ 
    ###PLOT CONSENSUS TREE
    NODOS <- as.numeric(rownames(M2[[2]]))
    N2 <- 1
    for(i in 1:1000){
      LN2 <- length(N2)
      ESTE <- which(N2==i)
      if(ESTE==1){
        N2 <- c(i*2,N2[1],(i*2)+1,N2[-1])
      }
      if(ESTE==LN2 & LN2!=1){
        N2 <- c(N2[-LN2], i*2, N2[LN2], (i*2)+1)
      }
      if(ESTE!=1 & ESTE!=LN2){
        N2 <- c(N2[1:(ESTE-1)], i*2, N2[ESTE], (i*2)+1, N2[(ESTE+1):LN2])
      }
    }
    
    
    N3 <- (gtools::even(N2))
    
    POS <- 1:length(N2)
    NODOS2 <- POS[match(NODOS, N2)]
    M2kk <-M2[[2]][order(NODOS2),]
    NODOS2 <- seq(0, 1, 1/(nrow(M2kk)+1)) ##PARA HACER LAS POSOCIONES DE LOS NODOS RELATIVAS
    par(mfrow=c(2, 1), mar =c(3, 4, 3, 1))
    plot(c(0,1), c(0, 1.01), type = "n", axes = FALSE, ylab = "Remanent deviance (%)", xlab = "")
    axis(2, las = 1, cex.axis = 1)
    
    for(i in 1:nrow(M2kk)){
      DEAQUI <- as.numeric(rownames(M2kk))[i]
      if(DEAQUI != 1){
        AACA <- floor(DEAQUI/2)
        ESTE <- which(as.numeric(rownames(M2kk))==AACA)
        lines(c(NODOS2[i+1], NODOS2[i+1]), M2kk[c(i, ESTE),3]/max(M2kk[,3]), lwd=2)
        if(gtools::even(DEAQUI)){
          AACA2 <- floor(DEAQUI/2)
          AACA <- DEAQUI+1
          ESTE <- which(as.numeric(rownames(M2kk))==AACA)
          ESTE2 <- which(as.numeric(rownames(M2kk))==AACA2)
          lines(c(NODOS2[i+1], NODOS2[ESTE+1]), M2kk[rep(ESTE2,2), 3]/max(M2kk[,3]), lwd=2)
          
          AACA <- floor(DEAQUI/2)
          ESTE <- which(as.numeric(rownames(M2kk))==AACA)
          text(mean(c(NODOS2[i+1], NODOS2[ESTE+1])), (M2kk[ESTE2, 3]/max(M2kk[,3]))+0.02, paste( M2kk[ESTE,1], "/", round(M2kk[ESTE,2], 2)))
        }
      }
      
    }
    axis(1, at=NODOS2[-c(1, length(NODOS2))][M2kk[,1]=="<leaf>"], labels = rownames(M2kk)[M2kk[,1]=="<leaf>"], tick = FALSE)
    #End of forest
    
    MEDIAS <- M2kk[,5][M2kk[,1]=="<leaf>"]
    ERROR <- M2kk[,6][M2kk[,1]=="<leaf>"]
    par(mar = c(4, 4,3, 1))
    if(is.numeric(data[,1])) meanerrba(MEDIAS, ERROR, names.arg =  rownames(M2kk)[M2kk[,1]=="<leaf>"], ylab = colnames(data)[1], xlab = "Node", ylim = c((min(MEDIAS-ERROR)), max((MEDIAS+ERROR))), cex.axis=1, cex.lab = 1) else {
      meanerrba(MEDIAS, ERROR, names.arg =  rownames(M2kk)[M2kk[,1]=="<leaf>"], ylab = "", xlab = "Node", ylim = c((min(MEDIAS-ERROR)), max((MEDIAS+ERROR))), cex.axis=1, cex.lab = 1, axes=FALSE)
      axis(2, labels=levels(data[,1]), at = 1:length(levels(data[,1])), las=1)
      abline(h=seq(0.5, length(levels(data[,1]))+0.5, 1), lty=2)
    }
  } 
  M2} else stop("First column must be numeric or a factor with two or more levels")
}#END FUNCTION