Delineamento ótimo para curva de crescimento de aves
Luciana Helena Kowalski & Walmes Marques Zeviani
Descrição do ensaio
TODO
Objetivos da análise
- Ajustar o modelo de crescimento logístico para cada unidade experimental.
- Obter o delineamento ótimo pelo método proposto por Santos e Santos (2008).
DOS SANTOS, M. I. R.; DOS SANTOS, P. M. R. Sequential experimental designs for nonlinear regression metamodels in simulation. Simulation Modelling Practice and Theory, v. 16, n. 9, p. 1365–1378, 2008. Elsevier BV. Disponível em: http://dx.doi.org/10.1016/j.simpat.2008.07.001.
Análise
#-----------------------------------------------------------------------
# Carregando pacotes.
library(latticeExtra)## Loading required package: lattice## Loading required package: RColorBrewerlibrary(reshape2)
library(nlme)
#-----------------------------------------------------------------------
# Leitura e preparação dos dados.
# list.files()
#
# da <- gdata::read.xls(
# xls = "~/Dropbox/Bolsista/Luciana/Dados_tratamentos.xlsx",
# header = TRUE)
# str(da)
#
# v <- names(da)[1:4]
# r <- grep(x = names(da), "^Peso_\\d+d$", value = TRUE)
#
# # Seleciona variáveis.
# db <- subset(da, select = c(v, r))
# str(db)
#
# # Empilha, de wide para long.
# db <- melt(data = db,
# id.vars = v,
# variable.name = "dias",
# value.name = "peso")
#
# # Extrai os dias da string.
# db$dias <- as.integer(gsub(x = as.character(db$dias),
# pattern = "\\D",
# replacement = ""))
#
# # Variáveis com nome em caixa baixa e com 4 digitos.
# names(db) <- tolower(names(db))
# names(db) <- substr(names(db), 1, 4)
# str(db)
#
# write.csv2(db,
# file = "cres-aves.csv",
# row.names = FALSE)
# args(write.table)Análise exploratória
#-----------------------------------------------------------------------
# Visualização dos dados.
rm(list = ls())
db <- read.csv2("cres-aves.csv")
str(db)## 'data.frame': 441 obs. of 6 variables:
## $ box : int 1 2 3 4 5 6 7 8 9 10 ...
## $ trea: Factor w/ 7 levels "T1","T2","T3",..: 6 3 4 1 6 4 5 6 2 3 ...
## $ repl: Factor w/ 9 levels "R1","R2","R3",..: 1 1 1 2 2 3 3 4 4 4 ...
## $ bloc: Factor w/ 7 levels "1","2","3","4",..: 2 2 2 3 3 4 4 1 1 1 ...
## $ dias: int 1 1 1 1 1 1 1 1 1 1 ...
## $ peso: num 0.043 0.0434 0.0433 0.0458 0.0461 ...xyplot(peso ~ dias | trea,
xlab = "Idade das aves (dias)",
ylab = "Peso por ave (kg)",
data = db,
type = "o",
groups = box,
as.table = TRUE)
xyplot(peso ~ dias | bloc,
xlab = "Idade das aves (dias)",
ylab = "Peso por ave (kg)",
data = db,
groups = trea,
auto.key = list(columns = 4),
as.table = TRUE) +
glayer(panel.smoother(...))
xyplot(peso ~ reorder(trea, peso) | factor(dias),
data = db,
xlab = "Tratamentos",
ylab = "Peso por ave (kg)",
type = c("p", "a"),
scales = list(y = list(relation = "free")),
as.table = TRUE)
Ajuste das curvas por unidade experimental
O modelo logístico na parametrização centro-escala será usado \[
f(x) = \dfrac{\theta_a}{1+\exp\left\{-\dfrac{(x-\theta_m)}{\theta_s}\right\}},
\] em que \(\theta_a\) é a assintota superior do modelo (Asym), \(\theta_m\) é o ponto de inflexão da curva (xmid) e \(\theta_s\) é o parâmetro de escala (scal) que proporcional a taxa no ponto de inflexão.
#-----------------------------------------------------------------------
# Ajuste de curvas de crescimento.
plot(peso ~ dias,
data = db,
xlab = "Idade (dias)",
ylab = "Peso da ave (kg)")
# Ajuste de uma única curva usando o modelo logístico
n0 <- nls(peso ~ SSlogis(dias, Asym, xmid, scal),
data = db)
summary(n0)##
## Formula: peso ~ SSlogis(dias, Asym, xmid, scal)
##
## Parameters:
## Estimate Std. Error t value Pr(>|t|)
## Asym 4.03624 0.04248 95.02 <2e-16 ***
## xmid 31.10929 0.21992 141.46 <2e-16 ***
## scal 8.32030 0.10017 83.06 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.07816 on 438 degrees of freedom
##
## Number of iterations to convergence: 0
## Achieved convergence tolerance: 5.204e-07# Curva ajustada sobre os dados.
plot(peso ~ dias,
data = db,
xlab = "Idade (dias)",
ylab = "Peso da ave (kg)")
with(as.list(coef(n0)), {
curve(Asym/(1 + exp(-(x - xmid)/scal)), add = TRUE)
})
# Estimativas dos parâmetros.
summary(n0)##
## Formula: peso ~ SSlogis(dias, Asym, xmid, scal)
##
## Parameters:
## Estimate Std. Error t value Pr(>|t|)
## Asym 4.03624 0.04248 95.02 <2e-16 ***
## xmid 31.10929 0.21992 141.46 <2e-16 ***
## scal 8.32030 0.10017 83.06 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.07816 on 438 degrees of freedom
##
## Number of iterations to convergence: 0
## Achieved convergence tolerance: 5.204e-07#-----------------------------------------------------------------------
# Ajustando curvas por unidade experimental.
db <- groupedData(peso ~ dias | box,
data = db,
order.groups = FALSE)
nn0 <- nlsList(peso ~ SSlogis(dias, Asym, xmid, scal) | box,
data = db)
# summary(nn0)
# Estimavas dos parâmetros.
coef(nn0)## Asym xmid scal
## 1 4.074180 30.93322 8.223688
## 2 4.215229 33.23719 8.740823
## 3 4.055735 31.29113 8.280800
## 4 3.725411 30.97364 8.351081
## 5 3.878549 30.45274 8.117228
## 6 4.019443 30.72285 8.387978
## 7 4.082601 30.07314 8.122272
## 8 3.988902 30.53358 8.020265
## 9 3.540491 29.33596 7.876164
## 10 4.194146 32.61104 8.533987
## 11 4.155080 31.04425 8.512133
## 12 4.030185 31.08774 8.290863
## 13 4.185917 30.75151 8.208222
## 14 4.024709 30.88894 8.347558
## 15 3.829236 30.30808 7.979090
## 16 4.139182 31.31767 8.452893
## 17 4.477828 32.22067 8.619995
## 18 4.350558 32.09547 8.646756
## 19 4.242481 32.49709 8.485671
## 20 3.830979 30.63237 8.337133
## 21 4.182363 32.74966 8.650834
## 23 3.846003 31.18225 8.454008
## 24 4.106182 30.87169 8.138230
## 25 4.211915 30.54462 8.062130
## 26 4.064900 32.06287 8.647269
## 27 4.054495 31.42865 8.479457
## 28 3.914903 30.85261 8.224941
## 29 3.929676 30.09685 8.067882
## 30 4.100888 31.00701 8.169904
## 31 3.952377 30.29044 7.930843
## 32 4.083623 32.18001 8.695355
## 33 4.203142 31.48549 8.566117
## 34 4.051004 31.45181 8.409545
## 35 4.028197 31.55440 8.324931
## 36 4.036118 32.00361 8.387689
## 37 4.071263 32.36570 8.415653
## 38 4.147470 31.55825 8.474382
## 39 3.991151 29.80485 7.968504
## 40 3.981902 30.50118 8.243026
## 41 4.005884 29.96387 8.101521
## 42 4.109818 30.75293 8.491337
## 43 4.264007 31.43032 8.492016
## 45 3.931001 31.04219 8.402710
## 46 3.988584 30.31369 8.163357
## 47 3.909101 29.96832 8.027854
## 48 3.824782 30.40474 8.121887
## 49 4.173473 32.86862 8.825741
## 50 3.999483 31.02423 8.310454
## 51 3.969915 30.53400 8.331385
## 52 3.759298 30.86721 8.202826
## 53 3.891787 30.57271 8.146667
## 54 4.517963 33.47308 8.328610
## 55 4.114298 30.81932 8.381084
## 56 3.805704 29.77064 8.048159
## 57 4.029266 31.03147 8.353626
## 58 3.872845 31.52011 8.207400
## 59 3.986676 31.60098 8.266596
## 60 4.263288 31.04489 8.168680
## 61 3.988179 30.90293 8.435877
## 62 4.146654 31.98475 8.466267
## 63 4.213107 29.80086 7.916188
## 64 4.094210 31.89360 8.736762
## 65 3.861390 30.83758 8.257330# Intervalos de confiança para os parâmetros.
plot(intervals(nn0), layout = c(NA, 1))
#-----------------------------------------------------------------------
# Diagnóstico.
xyplot(residuals(nn0) ~ fitted(nn0)) +
layer(panel.abline(h = 0, lty = 2))
xyplot(residuals(nn0) ~ dias | box,
data = db) +
layer(panel.abline(h = 0, lty = 2))
# ATTENTION: o padrão dos resíduos mostra que atendimento dos
# pressupostos não está satisfatório. É preciso usar um modelo mais
# flexível, por exemplo. Essa preocupação será resolvida depois.Obtenção do delineamento ótimo
#-----------------------------------------------------------------------
# Delineamento ótimo baseado na proposta de Santos & Santos.
#' @name santos_santos
#' @author Walmes Marques Zeviani (walmes@@ufpr.br).
#' @description Função que obtem o delineamento ótimo de acordo com a
#' metodologia de Santos \& Santos (2008).
#' @param n números de pontos para posicionar dentro do domínio
#' \code{rx}.
#' @param theta lista nomeada com os parâmetros da função \code{f()}.
#' @param rx amplitude do domínio em \code{x}, vetor de dois números,
#' extremo inferior e superior.
#' @param f função que depende o vetor \code{x} da lista nomeada
#' \code{theta}.
#' @return retorna um vetor de \code{n} pontos que são os valores ótimos
#' para \code{x}.
santos_santos <- function(n, theta, rx, f) {
# Mínimo e amplitude de x.
mx <- min(rx)
dx <- diff(rx)
# Para transformar os valores de x para [0, 1] ou o inverso.
resizex <- function(x, mx = mx, dx = dx, to01 = TRUE) {
if (to01) {
(x - mx)/dx
} else {
mx + dx * x
}
}
x <- rx
p <- 1
while (p <= n - 2) {
# Lista com argumentos.
a <- c(x = list(x),
theta)
# Vetor de coordenadas y.
ry <- do.call(f, args = a)
x <- resizex(x, mx, dx, TRUE)
# Para padronizar as coordenadas y.
m <- min(ry)
d <- abs(diff(range(ry)))
y <- (ry - m)/d
# Calcula as distâncias.
k <- length(x)
i <- 1
dd <- numeric(k - 1)
for (i in 1:(k - 1)) {
xy <- cbind(x[i:(i + 1)], y[i:(i + 1)])
dd[i] <- dist(x = xy)[1]
}
M <- which.max(dd)
xnew <- approxfun(x = y[M:(M + 1)],
y = x[M:(M + 1)])(
v = median(y[M:(M + 1)]))
x <- sort(c(x, xnew))
x <- resizex(x, mx, dx, FALSE)
p <- p + 1
}
xx <- seq(rx[1], rx[2], length.out = 30)
a <- c(x = list(xx), theta)
y <- do.call(f, args = a)
plot(y ~ xx, type = "o")
abline(v = x, lty = 2, col = 2)
return(x)
}
#-----------------------------------------------------------------------
# Domínio.
rx <- range(db$dias)
rx## [1] 1 42# Pontos usados (dias).
xcurren <- unique(sort(db$dias))
xcurren## [1] 1 7 14 21 28 35 42# Número de pontos que usados.
n <- length(xcurren)
n## [1] 7# Função do modelo logístico.
flogis <- function(x, Asym, xmid, scal) {
Asym/(1 + exp((xmid - x)/scal))
}
# Parâmetros estimados para a primeira unidade experimental.
theta <- as.list(coef(nn0)[1, ])
# Sequência de valores de x.
x <- seq(rx[1], rx[2], length.out = 30)
# Lista com x e theta.
a <- c(x = list(x), theta)
# Chama a função para obter os valores de y correspondentes dos x.
y <- do.call(flogis, args = a)
# Faz o gráfico.
plot(y ~ x, type = "o")
# Aplica o algoritmo de Santos e Santos.
santos_santos(n = n,
theta = theta,
rx = rx,
f = flogis)
## [1] 1.000 11.250 21.500 26.625 31.750 36.875 42.000#-----------------------------------------------------------------------
# Aplica para todas as undiades experimentais.
# Lista com os parâmetros de cada unidade experimental.
id <- rownames(coef(nn0))
theta_list <- by(data = coef(nn0),
INDICES = id,
FUN = function(x) {
as.list(unlist(x))
})
length(theta_list)## [1] 63# Retorna o delineamento para cada unidade experimental.
del <- lapply(theta_list,
FUN = function(theta) {
santos_santos(n = n,
theta = theta,
rx = rx,
f = flogis)
})# Converte lista para matriz.
del <- do.call(rbind, del)
# Verifica se existe variação.
apply(del, MARGIN = 2, FUN = sd)## [1] 0 0 0 0 0 0 0# Obtém os valores médios.
xoptmal <- apply(del, MARGIN = 2, FUN = mean)
# Delineamento atual vs delineamento ótimo.
cbind(xcurren, xoptmal)## xcurren xoptmal
## [1,] 1 1.000
## [2,] 7 11.250
## [3,] 14 21.500
## [4,] 21 26.625
## [5,] 28 31.750
## [6,] 35 36.875
## [7,] 42 42.000# Os delineamentos vistos no gráfico.
par(mfrow = c(1, 2))
with(as.list(coef(n0)), {
curve(flogis(x, Asym, xmid, scal),
from = 1,
to = 42,
main = "Atual")
abline(v = xcurren, col = "blue", lty = 2)
curve(flogis(x, Asym, xmid, scal),
from = 1,
to = 42,
main = "Ótimo")
abline(v = xoptmal, col = "red", lty = 2)
})
#-----------------------------------------------------------------------
# Determinar a eficiência relativa do delineamento contra a referência.
# Para retornar a matriz gradiente avaliada em x e theta.
der <- deriv3(expr = ~Asym/(1 + exp((xmid - x)/scal)),
namevec = c("Asym", "xmid", "scal"),
function.arg = function(x, Asym, xmid, scal) { NULL })
dets <- with(as.list(coef(n0)), {
# Obtém o valor da derivada primeira do modelo em relação a theta.
Xo <- attr(der(x = xoptmal, Asym, xmid, scal), which = "gradient")
Xc <- attr(der(x = xcurren, Asym, xmid, scal), which = "gradient")
# maximizar o det(X'X) corresponde ao critério D-ótimo.
c(xoptmal = det(crossprod(Xo)),
xcurren = det(crossprod(Xc)))
})
# Quanto mais ótimo baseado no D-ótimo.
dets[1]/dets[2]## xoptmal
## 1.129097# Aplica para todas as unidades experimentais.
opt <- sapply(theta_list,
FUN = function(theta) {
dets <- with(theta, {
Xo <- attr(der(x = xoptmal, Asym, xmid, scal),
which = "gradient")
Xc <- attr(der(x = xcurren, Asym, xmid, scal),
which = "gradient")
c(xoptmal = det(crossprod(Xo)),
xcurren = det(crossprod(Xc)))
})
dets[1]/dets[2]
})
# Ganho médio em eficiência considerando D-ótimo de 13%. Ou seja, o
# determinante de X'X é 13% maior considerando o delineamento ótimo em
# relação ao delineamento atual.
mean(opt)## [1] 1.129555sd(opt)## [1] 0.00921698#-----------------------------------------------------------------------