diff --git a/scripts/ce089-06.R b/scripts/ce089-06.R index 03e668c6c532cf0406ff98ce11e688529ea90552..bbd05609dd913191cee4bad5556531773cf0524d 100644 --- a/scripts/ce089-06.R +++ b/scripts/ce089-06.R @@ -63,8 +63,8 @@ cbind(n0 = rowMeans(b0) - coef(n0), n1 = rowMeans(b1) - coef(n1)) # Variância. -cbind(n0 = apply(b0, MARGIN = 1, var), - n1 = apply(b1, MARGIN = 1, var)) +round(cbind(n0 = apply(b0, MARGIN = 1, var), + n1 = apply(b1, MARGIN = 1, var)), digits = 5) # Função para calcular o erro quadrático médio. eqm <- function(x) { @@ -72,8 +72,8 @@ eqm <- function(x) { } # Erro quadrático médio. -cbind(n0 = apply(b0, 1, FUN = eqm), - n1 = apply(b1, 1, FUN = eqm)) +round(cbind(n0 = apply(b0, 1, FUN = eqm), + n1 = apply(b1, 1, FUN = eqm)), digits = 5) #----------------------------------------------------------------------- # Curvas ajustadas de onde é possÃvel determinar uma banda de confiança. @@ -205,6 +205,7 @@ str(st) # Gráfico de pares. pairs(b0) +pairs(b0$t) vcov(b0) cov2cor(vcov(b0)) @@ -381,3 +382,65 @@ rowSums(result)/N save.image(file = "my_results.RData") #----------------------------------------------------------------------- + +PaulaTb3.12 <- + structure(list( + vol = c(3.7, 3.5, 1.25, 0.75, 0.8, 0.7, 0.6, 1.1, 0.9, 0.9, 0.8, + 0.55, 0.6, 1.4, 0.75, 2.3, 3.2, 0.85, 1.7, 1.8, 0.4, + 0.95, 1.35, 1.5, 1.6, 0.6, 1.8, 0.95, 1.9, 1.6, 2.7, + 2.35, 1.1, 1.1, 1.2, 0.8, 0.95, 0.75, 1.3), + razao = c(0.825, 1.09, 2.5, 1.5, 3.2, 3.5, 0.75, 1.7, 0.75, + 0.45, 0.57, 2.75, 3, 2.33, 3.75, 1.64, 1.6, 1.415, + 1.06, 1.8, 2, 1.36, 1.35, 1.36, 1.78, 1.5, 1.5, 1.9, + 0.95, 0.4, 0.75, 0.03, 1.83, 2.2, 2, 3.33, 1.9, 1.9, + 1.625), + resp = c(1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, + 0, 1, 0, 0, 0, 0, 1, 0, 1, 0, 1, 0, 1, 0, 0, 1, 1, 1, + 0, 0, 1)), + .Names = c("vol", "razao", "resp"), + row.names = c(NA, 39L), + class = "data.frame") + +layout(1) + +# Visualização dos dados. +xyplot(resp ~ vol, + data = PaulaTb3.12, + groups = resp) + +# Ajuste do modelo considerando resposta binária. +m0 <- glm(resp ~ vol, + data = PaulaTb3.12, + family = binomial) +summary(m0) + +# Estimativa da DL_50 a partir do modelo ajustado. +dl <- -coef(m0)[1]/coef(m0)[2] + +# Curva ajustada. +plot(resp ~ vol, data = PaulaTb3.12) +curve(m0$family$linkinv(coef(m0)[1] + coef(m0)[2] * x), + add = TRUE) +abline(v = dl, h = 0.5, lty = 2) + +dl_50 <- function(dataset, index) { + m0 <- glm(resp ~ vol, + data = dataset[index, ], + family = binomial) + dl <- -coef(m0)[1]/coef(m0)[2] + return(dl) +} + +dl50 <- boot(data = PaulaTb3.12, + statistic = dl_50, + R = 2999) + +summary(dl50) + +hist(dl50) + +boot.ci(dl50, + type = "all", + index = c(1, 1)) + +#-----------------------------------------------------------------------