R/funcion_vb_total.R

In PabloMBooster/CAMPtools:

# funcion para calcular los parametros de crecimiento por cada lote

# falta incluir el vector de t_inicial, po ahor asume t_inicial = 5 dias

funcion_vb_total <- function(data = DSB_AB, sector = DSB_AB$Lote,
                             t_inicial = NULL){
  if(is.null(t_inicial)){
    t_inicial = 5
  }
  calculo_vb <- lapply(split(data, sector, drop = TRUE),function(x){
    print(x$Lote[1])
    Semana_inicio <- x$Semana[which.min(x$Fecha)]
    Semana_fin    <- x$Semana[which.max(x$Fecha)]
    Fecha_inicio  <- min(x$Fecha)
    Fecha_fin     <- max(x$Fecha)
    Paltos_muestreados  <- length(x$Fecha)
    Variedad     <- x$Variedad[1]
    Cantidad     <- x$Cantidad[1]
    Fundo        <- as.character(x$Fundo[1])
    Parcela      <- as.character(x$Parcela[1])
    Lote         <- as.character(x$Lote[1])
    LotesInpseccion <- x$LotesInpseccion[1]
    numero_dias = as.numeric(Fecha_fin - Fecha_inicio)
    nlen = length(unique(x$Fecha))

    # para calcula Winf y k es necesario 4 muestreos
    if(numero_dias > 20 & nlen > 3){
      Peso_medio <- tapply(x$Peso, x$Fecha, mean, na.rm = TRUE)
      Fecha      <- as.Date(names(Peso_medio), format = "%Y-%m-%d")
      t          <- as.numeric(diff(Fecha))
      t          <- c(t_inicial, t)
      t          <- cumsum(t)
      Peso_medio <- as.numeric(Peso_medio)

      datos = data.frame(t = t, Peso_medio = Peso_medio)

      model   = nls(Peso_medio ~ Winf*(1-exp(-k*(t)))^3, data = datos, control = list(maxiter = 1000),
                    start = c(Winf = 227, k = 0.02))

      Winf <- as.numeric(coefficients(model)[1])
      k    <- as.numeric(coefficients(model)[2])

    }else{
      Winf <- NA
      k    <- NA
    }

    cbind.data.frame(Semana_inicio, Semana_fin, Fecha_inicio, Fecha_fin,
                     Paltos_muestreados, Variedad, Cantidad, Fundo,
                     Parcela, Lote, LotesInpseccion, numero_dias, Winf, k)
  })

  require(dplyr)
  resultado <- suppressWarnings(calculo_vb %>% lapply(as.data.frame) %>% bind_rows())

  resultado = control_vb(resultado)
  write.csv(resultado, "parametros_vb.csv")
  return(resultado)
}

control_vb <- function(vb){
  k_outliers = remove_outliers(vb$k)
  vb$revisar = "bien"
  vb$revisar[is.na(k_outliers)] = "revisa"
  return(vb)
}