Configurar salida y cargar librerías

knitr::opts_chunk$set(echo = T, comment=" ", message=F, warning=F,options(scipen=999, root.dir = "D:/MANEJO FITOSANITARIO 2023/R_para_SEMINARIOS"))
#install.packages("ecotox")
library(ecotox)
#install.packages("openxlsx")
library(openxlsx)
#install.packages("ggplot2")
library(ggplot2)
#install.packages("dplyr")
library(dplyr)

## 
## Attaching package: 'dplyr'

## The following objects are masked from 'package:stats':
## 
##     filter, lag

## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union

#install.packages("knitr")
library(knitr)
#install.packages("DT")
library(DT)
#install.packages("gridExtra")
library(gridExtra)

## 
## Attaching package: 'gridExtra'

## The following object is masked from 'package:dplyr':
## 
##     combine

EJERCICIO 1

• Se aplicó un extracto botánico etanólico en 6 concentraciones, sobre una especie de insecto. Se solicito calcular las concentraciones 50 y 90.

Cargar datos

df<-as.data.frame(read.xlsx("D:/MANEJO FITOSANITARIO 2023/R_para_SEMINARIOS/CL50.xlsx", sheet=2))

Cálculo CL50 y CL90 para extracto (ext)

ext <- LC_probit((response / total) ~ log10(dose),
               p = c(50, 90),
               weights = total,
               data = df[df$nominal_dose != 0, ],
               subset = c(application == "extracto"))

Mostrar los resultados de las CL50 y CL90

DT::datatable(ext)

Graficar la curva dosis-respuesta

lc_ext <- subset(df, application %in% c("extracto"))
p1 <- ggplot(data = lc_ext[lc_ext$nominal_dose != 0, ],
              aes(x = log10(dose), y = (response / total))) +
   geom_point() + ggtitle("Gráfico Dosis-Probit (Extracto Botánico)")
   geom_smooth(method = "glm",
               method.args = list(family = binomial(link = "probit")),
               aes(weight = total), colour = "#FF0000", se = TRUE)

  mapping: weight = ~total 
  geom_smooth: na.rm = FALSE, orientation = NA, se = TRUE
  stat_smooth: na.rm = FALSE, orientation = NA, se = TRUE, method.args = list(family = list(family = "binomial", link = "probit", linkfun = function (mu) 
  qnorm(mu), linkinv = function (eta) 
  {
      thresh <- -qnorm(.Machine$double.eps)
      eta <- pmin(pmax(eta, -thresh), thresh)
      pnorm(eta)
  }, variance = function (mu) 
  mu * (1 - mu), dev.resids = function (y, mu, wt) 
  .Call(C_binomial_dev_resids, y, mu, wt), aic = function (y, n, mu, wt, dev) 
  {
      m <- if (any(n > 1)) 
          n
      else wt
      -2 * sum(ifelse(m > 0, (wt/m), 0) * dbinom(round(m * y), round(m), mu, log = TRUE))
  }, mu.eta = function (eta) 
  pmax(dnorm(eta), .Machine$double.eps), initialize = {
      if (NCOL(y) == 1) {
          if (is.factor(y)) 
              y <- y != levels(y)[1]
          n <- rep.int(1, nobs)
          y[weights == 0] <- 0
          if (any(y < 0 | y > 1)) 
              stop("y values must be 0 <= y <= 1")
          mustart <- (weights * y + 0.5)/(weights + 1)
          m <- weights * y
          if ("binomial" == "binomial" && any(abs(m - round(m)) > 0.001)) 
              warning(gettextf("non-integer #successes in a %s glm!", "binomial"), domain = NA)
      }
      else if (NCOL(y) == 2) {
          if ("binomial" == "binomial" && any(abs(y - round(y)) > 0.001)) 
              warning(gettextf("non-integer counts in a %s glm!", "binomial"), domain = NA)
          n <- (y1 <- y[, 1]) + y[, 2]
          y <- y1/n
          if (any(n0 <- n == 0)) 
              y[n0] <- 0
          weights <- weights * n
          mustart <- (n * y + 0.5)/(n + 1)
      }
      else stop(gettextf("for the '%s' family, y must be a vector of 0 and 1's\nor a 2 column matrix where col 1 is no. successes and col 2 is no. failures", "binomial"), domain = NA)
  }, validmu = function (mu) 
  all(is.finite(mu)) && all(mu > 0 & mu < 1), valideta = function (eta) 
  TRUE, simulate = function (object, nsim) 
  {
      ftd <- fitted(object)
      n <- length(ftd)
      ntot <- n * nsim
      wts <- object$prior.weights
      if (any(wts%%1 != 0)) 
          stop("cannot simulate from non-integer prior.weights")
      if (!is.null(m <- object$model)) {
          y <- model.response(m)
          if (is.factor(y)) {
              yy <- factor(1 + rbinom(ntot, size = 1, prob = ftd), labels = levels(y))
              split(yy, rep(seq_len(nsim), each = n))
          }
          else if (is.matrix(y) && ncol(y) == 2) {
              yy <- vector("list", nsim)
              for (i in seq_len(nsim)) {
                  Y <- rbinom(n, size = wts, prob = ftd)
                  YY <- cbind(Y, wts - Y)
                  colnames(YY) <- colnames(y)
                  yy[[i]] <- YY
              }
              yy
          }
          else rbinom(ntot, size = wts, prob = ftd)/wts
      }
      else rbinom(ntot, size = wts, prob = ftd)/wts
  }, dispersion = 1)), method = glm
  position_identity

p1

EJERCICIO 2

• Se realizaron dos ensayos con dos extractos botánicos (neem y huacatay). Se pide comparar las concentraciones letales 50 y 90 entre ambos extractos.

cargar datos

df2<-as.data.frame(read.xlsx("D:/MANEJO FITOSANITARIO 2023/R_para_SEMINARIOS/CL50.xlsx", sheet=3))

Cálculo CL50 y CL90 para el extracto de neem

neem <- LC_probit((response / total) ~ log10(dose),
               p = c(50, 90),
               weights = total,
               data = df2[df2$nominal_dose != 0, ],
               subset = c(application == "neem"))

CL50 y CL90 calculadas para el extracto de neem

DT::datatable(neem)

Cálculo CL50 y CL90 para el extracto de huacatay

huacatay <- LC_probit((response / total) ~ log10(dose),
               p = c(50, 90),
               weights = total,
               data = df2[df2$nominal_dose != 0, ],
               subset = c(application == "huacatay"))

CL50 y CL90 para el extracto de huacatay

DT::datatable(huacatay)

Comparación de los efectos de los dos extractos

-Uso de glm() para determinar los valores de CL usando el modelo probit para “neem” y “huacatay”

Extracto de neem

neem1 <- glm((response / total) ~ log10(dose),
         data = df2[df2$nominal_dose != 0, ],
         subset = c(application == "neem"),
         weights = total,
         family = binomial(link = "probit"))

Extracto de huacatay

huacatay1 <- glm((response / total) ~ log10(dose),
         data = df2[df2$nominal_dose != 0, ],
         subset = c(application == "huacatay"),
         weights = total,
         family = binomial(link = "probit"))

Test de razón (ratio test) para comparar los valores de CL50 y CL90.

ratios1 <- ratio_test(model_1 = neem1, model_2 = huacatay1, 
                     percentage = 50, 
                     compare = "neem - huacatay")
ratios2 <- ratio_test(model_1 = neem1, model_2 = huacatay1, 
                     percentage = 90, 
                     compare = "neem - huacatay")

Resultados de la comparación para la CL50

DT::datatable(ratios1)  

Resultados de la comparación para la CL90

DT::datatable(ratios2)

Graficar la curva dosis-respuesta

lc_neem <- subset(df2, application %in% c("neem"))
lc_huacatay <- subset(df2, application %in% c("huacatay"))

Gráfico dosis-respuesta extracto de neem

par(mfrow = c(1, 2))
n <- ggplot(data = lc_neem[lc_neem$nominal_dose != 0, ],
              aes(x = log10(dose), y = (response / total))) +
   geom_point() + ggtitle("Gráfico Dosis-Probit (Neem)") +
   geom_smooth(method = "glm",
               method.args = list(family = binomial(link = "probit")),
               aes(weight = total), colour = "#FF0000", se = TRUE)

Gráfico dosis-respuesta extracto de huacatay

h <- ggplot(data = lc_huacatay[lc_huacatay$nominal_dose != 0, ],
              aes(x = log10(dose), y = (response / total))) +
     geom_point() + ggtitle("Gráfico Dosis-Probit (Huacatay)") +
   geom_smooth(method = "glm",
               method.args = list(family = binomial(link = "probit")),
               aes(weight = total), colour = "#FF0000", se = TRUE)

Mostrar gráfico

grid.arrange(n, h, ncol=2)