EJERCICIOS DE TOXICOLOGIA
Victor Tello Mercado
2023-09-12
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.
Fijar directorio y 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_identityp1
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.
Fijar directorio y 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 neem
neem1 <- glm((response / total) ~ log10(dose),
data = df2[df2$nominal_dose != 0, ],
subset = c(application == "neem"),
weights = total,
family = binomial(link = "probit"))Extracto huacatay
huacatay1 <- glm((response / total) ~ log10(dose),
data = df2[df2$nominal_dose != 0, ],
subset = c(application == "huacatay"),
weights = total,
family = binomial(link = "probit"))Una vez que los modelos están ejecutados, se aplica el 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)