1 Instalação do pacote

Optimal Bayesian experimental design using the approximate coordinate exchange (ACE) algorithm.

O pacote [acebayes] tem funções para obtenção de delineamentos ótimos bayesianos usando o algorítmo de troca de coordenadas aproximado (approximate coordinate exchange algorithm). As funções do pacote são implementadas em baixo nível usando C++ por meio do pacote Rcpp.

Para instalar o pacote [acebayes] execute em usa sessão R as seguintes instruções.

install.packages("acebayes",
                 dependencies = TRUE,
                 repos = "https://cran-r.c3sl.ufpr.br/")

As instruções a carregam o pacote, exibem sua versão e conteúdo.

rm(list = objects())

library(acebayes)
packageVersion("acebayes")

## [1] '1.5'

ls("package:acebayes")

##  [1] "ace"                "aceglm"             "acenlm"             "acephase1"         
##  [5] "acephase2"          "assess"             "inideshlrnsel"      "inideshlrsig"      
##  [9] "inideslrnsel"       "inideslrsig"        "optdesbeetle"       "optdescomp15bad"   
## [13] "optdescomp15sig"    "optdescomp15sigDRS" "optdescomp18bad"    "optdeshlrbaa"      
## [17] "optdeshlrbad"       "optdeshlrnsel"      "optdeshlrsig"       "optdeslinmod"      
## [21] "optdeslrbaa"        "optdeslrbad"        "optdeslrnsel"       "optdeslrsig"       
## [25] "pace"               "paceglm"            "pacenlm"            "utilbeetle"        
## [29] "utilcomp15bad"      "utilcomp15sig"      "utilcomp15sigDRS"   "utilcomp18bad"     
## [33] "utilhlrbaa"         "utilhlrbad"         "utilhlrnsel"        "utilhlrsig"        
## [37] "utilityglm"         "utilitynlm"         "utillinmod"         "utillrbaa"         
## [41] "utillrbad"          "utillrnsel"         "utillrsig"

2 Importação e análise descritiva dos dados

Os dados estão arquivados em arquivo com campos separador pelo delimitador tabulação. O arquivo contém informação de 6 variáveis:

peso: é o peso ??peso médio de cada animal?? (kg).
dias: é a idade dos animais no dia da pesagem.
trea: corresponde ao tratamento aplicado.
box: é a unidade experimental e o valor corresponde ao identificador único das caixas onde são mantidos os animais.
bloc: fator de blocagem de condições ambientais e de operação dentro da unidade de produção.
repl: são as repetições. Os bloc agrupam essas repetições. Para efeito de modelagem, é melhor considerar repl como fator de blocagem.

library(tidyverse)

da <- read_tsv("dados_tratamentos.txt")
da

## # A tibble: 441 x 6
##      box trea  repl  bloc   dias   peso
##    <int> <chr> <chr> <chr> <int>  <dbl>
##  1     1 T6    R1    2         1 0.0430
##  2     2 T3    R1    2         1 0.0434
##  3     3 T4    R1    2         1 0.0433
##  4     4 T1    R2    3         1 0.0458
##  5     5 T6    R2    3         1 0.0460
##  6     6 T4    R3    4         1 0.0486
##  7     7 T5    R3    4         1 0.0489
##  8     8 T6    R4    1         1 0.0396
##  9     9 T2    R4    1         1 0.0398
## 10    10 T3    R4    1         1 0.0397
## # ... with 431 more rows

# da %>% xtabs(formula = ~trea + dias)
# da %>% xtabs(formula = ~trea + bloc)
da %>%
    xtabs(formula = peso ~ trea + repl + dias) %>%
    ftable() %>%
    round(digits = 4)

##           dias      1      7     14     21     28     35     42
## trea repl                                                      
## T1   R1        0.0431 0.1682 0.4340 0.9230 1.6200 2.2625 3.0450
##      R2        0.0458 0.1732 0.4376 0.8625 1.6050 2.2200 2.9725
##      R3        0.0487 0.1778 0.4318 0.9527 1.6365 2.3938 3.1338
##      R4        0.0396 0.1622 0.4165 0.9085 1.5782 2.2885 3.0141
##      R5        0.0454 0.1684 0.4411 0.9592 1.6368 2.3882 3.1776
##      R6        0.0430 0.1596 0.4090 0.9026 1.5667 2.3231 3.1179
##      R7        0.0419 0.1589 0.4150 0.9046 1.5895 2.3461 3.1263
##      R8        0.0461 0.1726 0.4312 0.9403 1.6077 2.3921 3.1947
##      R9        0.0504 0.1766 0.4345 0.8582 1.5897 2.2936 3.1385
## T2   R1        0.0428 0.1704 0.4452 0.9546 1.6625 2.3350 3.1275
##      R2        0.0460 0.1624 0.4356 0.9750 1.6306 2.3986 3.1057
##      R3        0.0490 0.1885 0.4622 1.0121 1.7115 2.4372 3.2000
##      R4        0.0398 0.1691 0.4283 0.9547 1.6382 2.3270 2.9757
##      R5        0.0459 0.1909 0.4938 0.9905 1.7711 2.4770 3.2905
##      R6        0.0427 0.1748 0.4424 0.9010 1.6231 2.3744 3.1551
##      R7        0.0412 0.1653 0.4387 0.9471 1.6447 2.4276 3.1276
##      R8        0.0457 0.1809 0.4604 0.9762 1.7205 2.3679 3.1859
##      R9        0.0505 0.1713 0.4510 0.9342 1.6947 2.3000 3.0934
## T3   R1        0.0434 0.1775 0.4456 0.8256 1.5385 2.2705 3.1013
##      R2        0.0468 0.1718 0.4356 0.9310 1.6026 2.3013 3.1167
##      R3        0.0490 0.1867 0.4427 0.9022 1.5324 2.3027 3.0946
##      R4        0.0397 0.1639 0.4399 0.8673 1.5865 2.3311 3.1689
##      R5        0.0462 0.1756 0.4319 0.9228 1.6213 2.2788 3.1218
##      R6        0.0428 0.1650 0.3964 0.8826 1.5342 2.3079 3.0434
##      R7        0.0426 0.1609 0.4018 0.8726 1.5474 2.3000 3.1092
##      R8        0.0459 0.1686 0.4224 0.9097 1.5934 2.3842 3.2184
##      R9        0.0501 0.1852 0.4722 0.8903 1.6333 2.3692 3.0829
## T4   R1        0.0433 0.1676 0.4454 0.9490 1.6474 2.4263 3.2041
##      R2        0.0461 0.1722 0.4497 0.9600 1.6667 2.3692 3.1821
##      R3        0.0486 0.1909 0.4778 1.0051 1.7054 2.4554 3.2139
##      R4        0.0403 0.1739 0.4410 0.9664 1.6500 2.4179 3.1447
##      R5        0.0460 0.1839 0.4631 1.0074 1.7333 2.4321 3.2897
##      R6        0.0434 0.1779 0.4650 0.9879 1.6947 2.4338 3.2122
##      R7        0.0425 0.1772 0.4550 0.9782 1.6684 2.4276 3.2378
##      R8        0.0461 0.1744 0.4663 0.9916 1.7250 2.4447 3.3378
##      R9        0.0503 0.1918 0.4735 0.9295 1.6575 2.3262 3.1463
## T5   R1        0.0434 0.1795 0.4539 1.0246 1.7388 2.4821 3.2513
##      R2        0.0460 0.1799 0.4684 1.0335 1.7946 2.6203 3.4171
##      R3        0.0489 0.1928 0.4968 1.0448 1.7974 2.5932 3.3417
##      R4        0.0398 0.1678 0.4221 0.9944 1.7077 2.4692 3.2855
##      R5        0.0460 0.1894 0.4909 1.0285 1.7372 2.4821 3.2872
##      R6        0.0440 0.1824 0.4768 0.9670 1.7475 2.4238 3.1550
##      R7        0.0416 0.1721 0.4612 1.0246 1.7405 2.5959 3.4000
##      R8        0.0461 0.1823 0.4779 1.0316 1.7908 2.5697 3.3066
##      R9        0.0506 0.1946 0.4952 1.0235 1.7862 2.4600 3.2875
## T6   R1        0.0430 0.1686 0.4542 0.9913 1.6868 2.4829 3.2553
##      R2        0.0460 0.1718 0.4479 0.9647 1.6566 2.4284 3.1446
##      R3        0.0485 0.1806 0.4579 0.9485 1.6679 2.4158 3.1797
##      R4        0.0396 0.1650 0.4414 0.9846 1.6838 2.4973 3.2378
##      R5        0.0458 0.1714 0.4460 0.9651 1.6692 2.4231 3.2038
##      R6        0.0429 0.1695 0.4474 0.9890 1.6846 2.4103 3.2077
##      R7        0.0418 0.1701 0.4609 0.9905 1.7205 2.4218 3.2654
##      R8        0.0463 0.1695 0.4536 1.0069 1.7397 2.4295 3.2308
##      R9        0.0507 0.1754 0.4669 1.0230 1.7514 2.4865 3.3472
## T7   R1        0.0429 0.1698 0.4392 0.9963 1.7141 2.5051 3.2987
##      R2        0.0458 0.1841 0.4595 1.0197 1.7079 2.5176 3.2122
##      R3        0.0493 0.1806 0.4624 1.0121 1.7105 2.5079 3.2171
##      R4        0.0397 0.1659 0.4258 0.9964 1.6974 2.5013 3.2397
##      R5        0.0462 0.1766 0.4532 0.7795 1.5679 2.4538 3.3269
##      R6        0.0433 0.1825 0.4681 1.0223 1.7808 2.5513 3.3692
##      R7        0.0428 0.1748 0.4755 1.0028 1.7563 2.5064 3.4244
##      R8        0.0458 0.1824 0.4754 1.0435 1.7962 2.5295 3.3013
##      R9        0.0507 0.1804 0.5066 1.0919 1.8656 2.7391 3.4891

da <- da %>%
    select(trea, repl, dias, peso) %>%
    mutate_at(.vars = 1:2, .fun = as_factor)

ggplot(data = da,
       mapping = aes(x = dias, y = peso, color = repl)) +
    geom_point() +
    geom_line() +
    facet_wrap(facets = ~trea, nrow = 2) +
    scale_color_discrete(name = "Blocos") +
    xlab("Idade dos animais (dias)") +
    ylab("Peso das unidades experimentais (kg)")

ggplot(data = da,
       mapping = aes(x = trea, y = peso, color = repl)) +
    geom_point() +
    geom_line(mapping = aes(group = repl)) +
    facet_wrap(facets = ~dias,
               scales = "free_y",
               nrow = 2) +
    scale_color_discrete(name = "Blocos") +
    xlab("Tratamentos") +
    ylab("Peso das unidades experimentais (kg)")

3 Ajuste da curva de crescimento

Para o ajuste da curva de crescimento considerou-se o modelo logístico que é amplamente usado na literatura

\[ f(x, \theta_{\text{ain}}, \theta_{\text{asp}}, \theta_{\text{xmd}}, \theta_{\text{scl}}) = \theta_{\text{ain}} + \frac{ \theta_{\text{asp}} - \theta_{\text{ain}} }{ 1 + \exp\left\{ -\frac{-(x - \theta_{\text{xmd}})}{\theta_{\text{scl}}} \right\} } \]

em que

\(x\) é a variável idade dos animais (dias),
\(\theta_{\text{ain}}\) é a assíntota inferior,
\(\theta_{\text{asp}}\) é a assintota superior e representa o peso máximo dos animais,
\(\theta_{\text{xmd}}\) é o tempo corresponde a troca de concavidade da curva de crescimento o onde está localizada a maior taxa de crescimento,
\(\theta_{\text{scl}}\) é o parâmetro inversamente proporcional à taxa de crescimento máxima.

O modelo foi inicialmente ajustado aos dados negligenciando a estrutura experimental de tratatmentos e repetições. Ou seja, foi ajustada uma única curva aos dados. Em seguida, foram efeito ajustes individualizados por unidade experimental, resultante da combinação de 7 tratamentos com 9 repetições que resulta em 63 unidades experimentais, portanto, 63 curvas ajustadas.

É feita análise do conjunto de curvas ajustadas para, juntamente com informações da literatura, permitir a elucidação de prioris para obtenção de delineamentos ótimos.

#-----------------------------------------------------------------------
# Ajuste de curvas de crescimento.

library(nlme)

plot(peso ~ dias, data = da)

# Modelo logístico.
n0 <- nls(peso ~ ain + (asp - ain)/(1 + exp(-(dias - xmd)/scl)),
          data = da,
          start = list(ain = 0.1, asp = 4, xmd = 30, scl = 8.5))
summary(n0)

## 
## Formula: peso ~ ain + (asp - ain)/(1 + exp(-(dias - xmd)/scl))
## 
## Parameters:
##     Estimate Std. Error t value Pr(>|t|)    
## ain -0.16272    0.01797  -9.056   <2e-16 ***
## asp  4.53927    0.08759  51.824   <2e-16 ***
## xmd 32.74626    0.36576  89.530   <2e-16 ***
## scl 10.11048    0.24069  42.007   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.06904 on 437 degrees of freedom
## 
## Number of iterations to convergence: 4 
## Achieved convergence tolerance: 5.112e-06

# Verifica se o modelo ficou bem ajustado.
plot(peso ~ dias, data = da)
with(as.list(coef(n0)), {
    curve(ain + (asp - ain)/(1 + exp(-(dias - xmd)/scl)),
          xname = "dias",
          add = TRUE)
})

#-----------------------------------------------------------------------
# Ajustando curvas por unidade experimental.

da <- da %>%
    mutate(ue = interaction(trea, repl, sep = "_"))
nlevels(da$ue)

## [1] 63

# Cria estrutura de dados agrupados.
db <- groupedData(peso ~ dias | ue,
                  data = da,
                  order.groups = FALSE)

nn0 <- nlsList(peso ~ ain + (asp - ain)/(1 + exp(-(dias - xmd)/scl)) | ue,
               data = db,
               start = as.list(coef(n0)))
# summary(nn0)
coef(nn0)

##               ain      asp      xmd       scl
## T6_R1 -0.16333610 4.554392 32.43074  9.954122
## T3_R1 -0.13063157 4.811604 35.40671 10.437553
## T4_R1 -0.15091560 4.529505 32.82732  9.944704
## T1_R1 -0.19261719 4.477241 33.41983 10.733849
## T5_R1 -0.17830938 4.465631 31.75738  9.999528
## T2_R1 -0.19640872 4.553903 33.16472 10.637944
## T7_R1 -0.15795452 4.562483 32.26068  9.776262
## T6_R2 -0.13423217 4.234861 31.54815  9.543785
## T3_R2 -0.19331188 4.805393 34.71475 10.995737
## T4_R2 -0.19272514 4.706803 33.65056 10.681716
## T1_R2 -0.14322841 4.167806 32.54291 10.066864
## T5_R2 -0.14841552 4.611977 31.68129  9.516348
## T2_R2 -0.14809032 4.215254 31.61772  9.705690
## T7_R2 -0.13366301 4.230716 30.88111  9.370584
## T6_R3 -0.14149272 4.425614 32.40658  9.871479
## T3_R3 -0.15195621 4.824545 35.24312 10.762038
## T4_R3 -0.17307061 4.523304 32.33693 10.257585
## T1_R3 -0.13277822 4.296586 32.07881  9.682581
## T5_R3 -0.15106194 4.462281 31.14844  9.610806
## T2_R3 -0.17778792 4.472030 32.13132 10.237776
## T7_R3 -0.13478147 4.263430 31.06218  9.434976
## T6_R4 -0.14181794 4.366174 31.65062  9.474751
## T3_R4 -0.15178904 4.814537 34.78013 10.382358
## T4_R4 -0.14674046 4.293270 31.82786  9.720289
## T1_R4 -0.14649302 4.187099 32.30733  9.879069
## T5_R4 -0.17490085 4.626501 32.64491 10.011881
## T2_R4 -0.13454293 3.833613 30.17812  9.285282
## T7_R4 -0.13997930 4.306504 31.30847  9.337489
## T6_R5 -0.15989281 4.519177 32.66882 10.040999
## T3_R5 -0.20526217 4.904282 35.15118 11.223116
## T4_R5 -0.22677765 4.999484 34.15551 11.104652
## T1_R5 -0.17332250 4.624103 33.37436 10.371286
## T5_R5 -0.21422512 4.835501 33.25984 10.837941
## T2_R5 -0.18503425 4.674221 32.60530 10.363592
## T7_R5 -0.04083950 4.694827 34.04613  8.814270
## T6_R6 -0.19612309 4.639687 33.03353 10.507609
## T3_R6 -0.11051401 4.223222 32.70462  9.492299
## T4_R6 -0.18281991 4.575191 32.67042 10.337570
## T1_R6 -0.14826872 4.567957 33.83981 10.138441
## T5_R6 -0.14809260 4.166163 30.84005  9.582782
## T2_R6 -0.13151402 4.461571 33.00409  9.828956
## T7_R6 -0.17528508 4.693083 32.28041  9.998157
## T6_R7 -0.22153272 4.885866 33.79881 10.909743
## T3_R7 -0.13570897 4.588863 34.17550 10.056190
## T4_R7 -0.18858087 4.791740 33.70736 10.594277
## T1_R7 -0.14051428 4.450743 33.16001  9.877093
## T5_R7 -0.16095615 4.739470 32.45892  9.799530
## T2_R7 -0.12212962 4.139530 31.23521  9.259975
## T7_R7 -0.23345342 5.414090 35.28904 11.227524
## T6_R8 -0.19757175 4.542767 32.27162 10.334930
## T3_R8 -0.14216417 4.798749 34.38883 10.161733
## T4_R8 -0.23115397 5.260526 35.15908 11.285110
## T1_R8 -0.15379328 4.698546 33.85374 10.244798
## T5_R8 -0.14417251 4.332603 30.74093  9.369203
## T2_R8 -0.20798336 4.642604 33.10396 10.762673
## T7_R8 -0.17837088 4.455072 31.26362  9.877058
## T6_R9 -0.21437963 4.993552 33.79924 10.820348
## T3_R9 -0.09994577 4.147492 31.77101  9.384206
## T4_R9 -0.19074202 4.812797 34.46329 11.042940
## T1_R9 -0.13889685 4.765778 34.83393 10.385254
## T5_R9 -0.21765623 4.778885 32.92761 10.840717
## T2_R9 -0.18132788 4.369300 32.45687 10.391279
## T7_R9 -0.14196828 4.541478 30.63144  9.207198

# Intervalos de confiança para os parâmetros em cada ajuste.
plot(intervals(nn0))

# Gráfico de pares das estimativas dos parâmetros.
pairs(nn0)

# Obteção das médias e desvios-padrões (útil para proposição das
# prioris).
apply(coef(nn0), MARGIN = 2, mean)

##        ain        asp        xmd        scl 
## -0.1635557  4.5611425 32.8275049 10.1235952

apply(coef(nn0), MARGIN = 2, sd)

##        ain        asp        xmd        scl 
## 0.03491849 0.27963732 1.31456514 0.58523027

theta <- coef(nn0) %>%
    as.data.frame() %>%
    rownames_to_column() %>%
    as_tibble() %>%
    separate(col = "rowname", into = c("treat", "repl"))
theta

## # A tibble: 63 x 6
##    treat repl     ain   asp   xmd   scl
##    <chr> <chr>  <dbl> <dbl> <dbl> <dbl>
##  1 T6    R1    -0.163  4.55  32.4  9.95
##  2 T3    R1    -0.131  4.81  35.4 10.4 
##  3 T4    R1    -0.151  4.53  32.8  9.94
##  4 T1    R1    -0.193  4.48  33.4 10.7 
##  5 T5    R1    -0.178  4.47  31.8 10.00
##  6 T2    R1    -0.196  4.55  33.2 10.6 
##  7 T7    R1    -0.158  4.56  32.3  9.78
##  8 T6    R2    -0.134  4.23  31.5  9.54
##  9 T3    R2    -0.193  4.81  34.7 11.0 
## 10 T4    R2    -0.193  4.71  33.7 10.7 
## # ... with 53 more rows

theta %>%
    group_by(treat) %>%
    summarise_if(is.numeric, sd)

## # A tibble: 7 x 5
##   treat    ain   asp   xmd   scl
##   <chr>  <dbl> <dbl> <dbl> <dbl>
## 1 T1    0.0190 0.216 0.868 0.322
## 2 T2    0.0313 0.270 0.999 0.557
## 3 T3    0.0346 0.281 1.23  0.627
## 4 T4    0.0290 0.288 1.08  0.536
## 5 T5    0.0282 0.220 0.926 0.546
## 6 T6    0.0334 0.240 0.810 0.517
## 7 T7    0.0513 0.360 1.58  0.688

4 Obtenção dos delineamentos ótimos

Os delineamentos serão avaliados considerando-se 7 pontos de suporte pois, como o cliclo de vida das aves é de 42 dias, com 7 pontos de suporte tem-se observações semanalmente espaçadas.

Serão comparados 6 delineamentos

Delineamento regularmente espaçado de 7 dias: {1, 7, 14, 21, 28, 35, 42}. Para ser uma progressão geométrica perfeita, o primeiro ponto de suporte deveria ser o 0. No entanto, esse é conjunto de pontos de suporte usado na prática.
Delineamento obtido pela abordagem de Santos & Santos (2008).
Delineamento Bayesiano com prioris uniformes para os parâmetros.
Delineamento Bayesiano com prioris Gaussianas para os parâmetros.
Delineamento 3 mas com restrição de intervalo mínimo de 5 dias entre observações.
Delineamento 4 mas com restrição de intervalo mínimo de 5 dias entre observações.

Todos os delineamentos serão avaliados dentro do domínimo real limitado em \([1, 42]\). A função utilidade ou critério de otimabilidade será o D-ótimo.

#-----------------------------------------------------------------------
# Configurações iniciais.

# Número de pontos de suporte.
sup_points <- 7

# Número de fatores (apenas o tempo).
n_fat <- 1

# Extremos do domínio.
dom <- range(da$dias)

# Delineamento regularmente espaçado.
del_re <- cbind(t = c(1, seq(7, 42, by = 7)))
del_re

##       t
## [1,]  1
## [2,]  7
## [3,] 14
## [4,] 21
## [5,] 28
## [6,] 35
## [7,] 42

4.1 Delimeanto ótimo Bayesiano com prioris uniformes

As prioris usadas foram elucidades via análise exploratória dos ajustes feitos por unidade experimental. Com isso delimitou-se os extremos das distribuições uniformes para cada parâmetro do modelo logístico. As prioris são independentes entre os parâmetros.

\[ \begin{align*} \theta_{\text{ain}} &\sim U(-0.5, 0.1)\\ \theta_{\text{asp}} &\sim U(3.5, 6)\\ \theta_{\text{xmd}} &\sim U(28, 38)\\ \theta_{\text{scl}} &\sim U(7, 14) \end{align*} \]

#-----------------------------------------------------------------------
# Priores uniformes para os parâmtros.

# Desenho inicial (v.a. uniforme).
start.d <- del_re

# Priores uniformes para os parâmetro, baseadas nos ajustes por unidade
# experimental.
prior_u <- list(support = cbind(ain = c(-0.5, 0.1),
                                asp = c(3.5, 6),
                                xmd = c(28, 38),
                                scl = c(7, 14)))

set.seed(1)
del_uni <- acenlm(formula = ~ain + (asp - ain)/(1 + exp(-(t - xmd)/scl)),
                  start.d = start.d,
                  prior = prior_u,
                  criterion = "D",
                  N1 = 20,
                  N2 = 0,
                  lower = min(dom),
                  upper = max(dom),
                  progress = TRUE)

## Phase I iteration 1 out of 20 (Current value = -11.45013) 
## Phase I iteration 2 out of 20 (Current value = -11.41656) 
## Phase I iteration 3 out of 20 (Current value = -11.40905) 
## Phase I iteration 4 out of 20 (Current value = -11.40807) 
## Phase I iteration 5 out of 20 (Current value = -11.40601) 
## Phase I iteration 6 out of 20 (Current value = -11.40563) 
## Phase I iteration 7 out of 20 (Current value = -11.40562) 
## Phase I iteration 8 out of 20 (Current value = -11.40562) 
## Phase I iteration 9 out of 20 (Current value = -11.40562) 
## Phase I iteration 10 out of 20 (Current value = -11.40562) 
## Phase I iteration 11 out of 20 (Current value = -11.40562) 
## Phase I iteration 12 out of 20 (Current value = -11.40562) 
## Phase I iteration 13 out of 20 (Current value = -11.40562) 
## Phase I iteration 14 out of 20 (Current value = -11.40562) 
## Phase I iteration 15 out of 20 (Current value = -11.40562) 
## Phase I iteration 16 out of 20 (Current value = -11.40562) 
## Phase I iteration 17 out of 20 (Current value = -11.40562) 
## Phase I iteration 18 out of 20 (Current value = -11.40561) 
## Phase I iteration 19 out of 20 (Current value = -11.40561) 
## Phase I iteration 20 out of 20 (Current value = -11.40561)

del_uni

## Non Linear Model 
## Criterion = Bayesian D-optimality 
## Formula: ~ain + (asp - ain)/(1 + exp(-(t - xmd)/scl))
## Method:  quadrature 
## 
## nr =  2 , nq =  8 
## Prior: uniform
## 
## Number of runs = 7
## 
## Number of factors = 1
## 
## Number of Phase I iterations = 20
## 
## Number of Phase II iterations = 0
## 
## Computer time = 00:00:09

# Mostra os pontos de suporte obtidos.
cbind(t = sort(del_uni$phase1.d))

##             t
## [1,]  1.00000
## [2,]  1.00000
## [3,] 18.62386
## [4,] 18.62604
## [5,] 33.02671
## [6,] 33.03921
## [7,] 42.00000

plot(density(del_uni$phase1.d, bw = 1),
     col = "black",
     xlab = "Pontos de suporte",
     ylab = "Densidade",
     main = NA)
rug(del_uni$phase1.d)

4.2 Delimeanto ótimo Bayesiano com prioris uniformes e restrição

O deliamento a seguir preserva as configurações anteriores. No entanto, é incluída a restrição de que as obervações consecutivas devem estar espaçadas de pelo menos 5 dias.

# Fornece a restrição de espaçamento mínimo igual ou superior a 5 dias.
limits <- function(d, i, j) {
    grid <- seq(from = 1,
                to = 42,
                length.out = 300)
    for (s in as.vector(d)[-i]) {
        grid <- grid[(grid <= (s - 5)) | (grid >= (s + 5))]
    }
    grid
}

del_uni_lim <- acenlm(formula = ~ain + (asp - ain)/(1 + exp(-(t - xmd)/scl)),
                      start.d = start.d,
                      prior = prior_u,
                      criterion = "D",
                      N1 = 20,
                      N2 = 0,
                      limits = limits,
                      lower = min(dom),
                      upper = max(dom),
                      progress = TRUE)

## Phase I iteration 1 out of 20 (Current value = -11.70372) 
## Phase I iteration 2 out of 20 (Current value = -11.70372) 
## Phase I iteration 3 out of 20 (Current value = -11.70372) 
## Phase I iteration 4 out of 20 (Current value = -11.70372) 
## Phase I iteration 5 out of 20 (Current value = -11.70372) 
## Phase I iteration 6 out of 20 (Current value = -11.70372) 
## Phase I iteration 7 out of 20 (Current value = -11.70372) 
## Phase I iteration 8 out of 20 (Current value = -11.70372) 
## Phase I iteration 9 out of 20 (Current value = -11.70372) 
## Phase I iteration 10 out of 20 (Current value = -11.70372) 
## Phase I iteration 11 out of 20 (Current value = -11.70372) 
## Phase I iteration 12 out of 20 (Current value = -11.70372) 
## Phase I iteration 13 out of 20 (Current value = -11.70372) 
## Phase I iteration 14 out of 20 (Current value = -11.70372) 
## Phase I iteration 15 out of 20 (Current value = -11.70372) 
## Phase I iteration 16 out of 20 (Current value = -11.70372) 
## Phase I iteration 17 out of 20 (Current value = -11.70372) 
## Phase I iteration 18 out of 20 (Current value = -11.70372) 
## Phase I iteration 19 out of 20 (Current value = -11.70372) 
## Phase I iteration 20 out of 20 (Current value = -11.70372)

del_uni_lim

## Non Linear Model 
## Criterion = Bayesian D-optimality 
## Formula: ~ain + (asp - ain)/(1 + exp(-(t - xmd)/scl))
## Method:  quadrature 
## 
## nr =  2 , nq =  8 
## Prior: uniform
## 
## Number of runs = 7
## 
## Number of factors = 1
## 
## Number of Phase I iterations = 20
## 
## Number of Phase II iterations = 0
## 
## Computer time = 00:00:07

# Mostra os pontos de suporte obtidos.
cbind(t = sort(del_uni_lim$phase1.d))

##              t
## [1,]  1.000000
## [2,]  6.073579
## [3,] 15.946488
## [4,] 21.000000
## [5,] 29.933110
## [6,] 35.000000
## [7,] 42.000000

plot(density(del_uni_lim$phase1.d, bw = 1),
     col = "black",
     xlab = "Pontos de suporte",
     ylab = "Densidade",
     main = NA)
rug(del_uni_lim$phase1.d)

4.3 Delimeanto ótimo Bayesiano com prioris Gaussianas

As prioris usadas foram elucidades via análise exploratória dos ajustes feitos por unidade experimental. Com isso determinou-se a média e variância para as distribuições Gaussianas de cada parâmetro do modelo logístico. As prioris são independentes entre os parâmetros.

\[ \begin{align*} \theta_{\text{ain}} &\sim N(-0.15, 0.04^2)\\ \theta_{\text{asp}} &\sim N(4.2, 0.3^2)\\ \theta_{\text{xmd}} &\sim N(32, 1^2)\\ \theta_{\text{scl}} &\sim N(10, 0.6^2) \end{align*} \]

#-----------------------------------------------------------------------
# Priores gaussianas com método de Monte Carlo.

prior_g <- list(mu = c(ain = -0.15,
                       asp = 4.2,
                       xmd = 32,
                       scl = 10),
                sigma2 = c(ain = 0.04,
                           asp = 0.3,
                           xmd = 1,
                           scl = 0.6)^2)

del_gau <- acenlm(formula = ~ain + (asp - ain)/(1 + exp(-(t - xmd)/scl)),
                  start.d = start.d,
                  prior = prior_g,
                  criterion = "D",
                  N1 = 20,
                  N2 = 0,
                  lower = min(dom),
                  upper = max(dom),
                  progress = TRUE)

## Phase I iteration 1 out of 20 (Current value = -11.15599) 
## Phase I iteration 2 out of 20 (Current value = -11.13244) 
## Phase I iteration 3 out of 20 (Current value = -11.13204) 
## Phase I iteration 4 out of 20 (Current value = -11.13193) 
## Phase I iteration 5 out of 20 (Current value = -11.13192) 
## Phase I iteration 6 out of 20 (Current value = -11.13192) 
## Phase I iteration 7 out of 20 (Current value = -11.13192) 
## Phase I iteration 8 out of 20 (Current value = -11.13192) 
## Phase I iteration 9 out of 20 (Current value = -11.13192) 
## Phase I iteration 10 out of 20 (Current value = -11.13192) 
## Phase I iteration 11 out of 20 (Current value = -11.13192) 
## Phase I iteration 12 out of 20 (Current value = -11.13192) 
## Phase I iteration 13 out of 20 (Current value = -11.13192) 
## Phase I iteration 14 out of 20 (Current value = -11.13192) 
## Phase I iteration 15 out of 20 (Current value = -11.13192) 
## Phase I iteration 16 out of 20 (Current value = -11.13192) 
## Phase I iteration 17 out of 20 (Current value = -11.13192) 
## Phase I iteration 18 out of 20 (Current value = -11.13192) 
## Phase I iteration 19 out of 20 (Current value = -11.13192) 
## Phase I iteration 20 out of 20 (Current value = -11.13192)

del_gau

## Non Linear Model 
## Criterion = Bayesian D-optimality 
## Formula: ~ain + (asp - ain)/(1 + exp(-(t - xmd)/scl))
## Method:  quadrature 
## 
## nr =  2 , nq =  8 
## Prior: normal
## 
## Number of runs = 7
## 
## Number of factors = 1
## 
## Number of Phase I iterations = 20
## 
## Number of Phase II iterations = 0
## 
## Computer time = 00:00:09

# Mostra os pontos de suporte obtidos.
cbind(t = sort(del_gau$phase1.d))

##             t
## [1,]  1.00000
## [2,]  1.00000
## [3,] 18.36051
## [4,] 18.36064
## [5,] 32.74902
## [6,] 32.75817
## [7,] 42.00000

plot(density(del_gau$phase1.d, bw = 1),
     col = "black",
     xlab = "Pontos de suporte",
     ylab = "Densidade",
     main = NA)
rug(del_gau$phase1.d)

4.4 Delimeanto ótimo Bayesiano com prioris Gaussianas e restrição

O deliamento a seguir preserva as configurações anteriores. No entanto, é incluída a restrição de que as obervações consecutivas devem estar espaçadas de pelo menos 5 dias.

del_gau_lim <- acenlm(formula = ~ain + (asp - ain)/(1 + exp(-(t - xmd)/scl)),
                      start.d = start.d,
                      prior = prior_g,
                      criterion = "D",
                      N1 = 20,
                      N2 = 0,
                      limits = limits,
                      lower = min(dom),
                      upper = max(dom),
                      progress = TRUE)

## Phase I iteration 1 out of 20 (Current value = -11.45578) 
## Phase I iteration 2 out of 20 (Current value = -11.45578) 
## Phase I iteration 3 out of 20 (Current value = -11.45578) 
## Phase I iteration 4 out of 20 (Current value = -11.45578) 
## Phase I iteration 5 out of 20 (Current value = -11.45578) 
## Phase I iteration 6 out of 20 (Current value = -11.45578) 
## Phase I iteration 7 out of 20 (Current value = -11.45578) 
## Phase I iteration 8 out of 20 (Current value = -11.45578) 
## Phase I iteration 9 out of 20 (Current value = -11.45578) 
## Phase I iteration 10 out of 20 (Current value = -11.45578) 
## Phase I iteration 11 out of 20 (Current value = -11.45578) 
## Phase I iteration 12 out of 20 (Current value = -11.45578) 
## Phase I iteration 13 out of 20 (Current value = -11.45578) 
## Phase I iteration 14 out of 20 (Current value = -11.45578) 
## Phase I iteration 15 out of 20 (Current value = -11.45578) 
## Phase I iteration 16 out of 20 (Current value = -11.45578) 
## Phase I iteration 17 out of 20 (Current value = -11.45578) 
## Phase I iteration 18 out of 20 (Current value = -11.45578) 
## Phase I iteration 19 out of 20 (Current value = -11.45578) 
## Phase I iteration 20 out of 20 (Current value = -11.45578)

del_gau_lim

## Non Linear Model 
## Criterion = Bayesian D-optimality 
## Formula: ~ain + (asp - ain)/(1 + exp(-(t - xmd)/scl))
## Method:  quadrature 
## 
## nr =  2 , nq =  8 
## Prior: normal
## 
## Number of runs = 7
## 
## Number of factors = 1
## 
## Number of Phase I iterations = 20
## 
## Number of Phase II iterations = 0
## 
## Computer time = 00:00:06

# Mostra os pontos de suporte obtidos.
cbind(t = sort(del_gau_lim$phase1.d))

##              t
## [1,]  1.000000
## [2,]  6.073579
## [3,] 15.946488
## [4,] 21.000000
## [5,] 29.933110
## [6,] 35.000000
## [7,] 42.000000

plot(density(del_gau_lim$phase1.d, bw = 1),
     col = "black",
     xlab = "Pontos de suporte",
     ylab = "Densidade",
     main = NA)
rug(del_gau_lim$phase1.d)

4.5 Delineamento de Santos e Santos

flogis <- function(x, ain, asp, xmd, scl) {
    ain + (asp - ain)/(1 + exp(-(x - xmd)/scl))
}

santos_santos <- function(n, theta, rx, f) {
    #-------------------------------------------------------------------
    #' @param n números de pontos para posicionar dentro do domínio.
    #' @param theta lista nomeada com os parâmetros da curva f.
    #' @param rx amplitude do domínio em x, vetor de dois números.
    #' @param f função que depende o vetor x da lista nomeada theta.
    #' @return retorna os pontos 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)
}

del_ss <- santos_santos(n = sup_points,
                        theta = prior_g$mu,
                        rx = dom,
                        f = flogis)
del_ss <- cbind(t = del_ss)
del_ss

##           t
## [1,]  1.000
## [2,] 11.250
## [3,] 21.500
## [4,] 26.625
## [5,] 31.750
## [6,] 36.875
## [7,] 42.000

4.6 Comparação dos resultados

# Delineamentos bayesianos obtidos.
del_baye <- list(del_uni,
                 del_uni_lim,
                 del_gau,
                 del_gau_lim)
names(del_baye) <- c("U", "UR", "N", "NR")

del_sup <- c(list(del_re,
                  del_ss),
             lapply(del_baye, FUN = "[[", "phase1.d"))
names(del_sup) <- c(c("RE", "SS"), names(del_baye))
# del_sup

# Tabela com os pontos de suporte de cada delineamento.
tb_del <- lapply(del_sup, FUN = as.vector)
tb_del <- lapply(tb_del, FUN = sort)
tb_del <- do.call(cbind, tb_del)
knitr::kable(tb_del,
             digits = rep(2, 6),
             caption = "Tabela 1. Delineamentos obtidos com as metodologias consideradas. RE: regularmente espaçaço; SS: Santos & Santos (2008); U: Bayesiano com priores uniformes; UR: Bayesiano com priores uniformes com restrições; N: Bayesiano com priores Gaussianas; NR: Bayesiano com priores Gaussianas e restrição;")

Tabela 1. Delineamentos obtidos com as metodologias consideradas. RE: regularmente espaçaço; SS: Santos & Santos (2008); U: Bayesiano com priores uniformes; UR: Bayesiano com priores uniformes com restrições; N: Bayesiano com priores Gaussianas; NR: Bayesiano com priores Gaussianas e restrição;
RE	SS	U	UR	N	NR
1	1.00	1.00	1.00	1.00	1.00
7	11.25	1.00	6.07	1.00	6.07
14	21.50	18.62	15.95	18.36	15.95
21	26.62	18.63	21.00	18.36	21.00
28	31.75	33.03	29.93	32.75	29.93
35	36.88	33.04	35.00	32.76	35.00
42	42.00	42.00	42.00	42.00	42.00

# Obtenção das eficiências relativas com relação aos delineamento
# Bayesianos.

# Eficiencias absolutas.
util <- sapply(del_baye,
               FUN = function(baye) {
                   sapply(del_sup[c(3:6, 1:2)],
                          FUN = function(sup) {
                              baye$utility(d = sup)
                          })
               })

# assess(d1 = del_uni, d2 = del_uni_lim)

# Eficiências relativas.
util_rel <- util
for (i in seq_along(del_baye)) {
    util_rel[, i] <- 100 * util_rel[i, i]/util_rel[, i]
}

# Amplitude de variação em eficiência absoluta.
range(util)

## [1] -11.82661 -11.13192

# Amplitude de variação em eficiência relativa.
range(util_rel)

## [1]  95.86157 102.90926

knitr::kable(util,
             digits = 3,
             caption = "Tabela 2. Eficiência absoluta dos delineamentos (linhas) em relação a cada um dos delineamentos Bayesianos (colunas).")

Tabela 2. Eficiência absoluta dos delineamentos (linhas) em relação a cada um dos delineamentos Bayesianos (colunas).
	U	UR	N	NR
U	-11.406	-11.406	-11.134	-11.134
UR	-11.704	-11.704	-11.456	-11.456
N	-11.407	-11.407	-11.132	-11.132
NR	-11.704	-11.704	-11.456	-11.456
RE	-11.821	-11.821	-11.590	-11.590
SS	-11.827	-11.827	-11.612	-11.612

# Tabela com as efiências relativas (linhas) comparadas com cada
# delineamento bayesiano (colunas).
knitr::kable(util_rel,
             digits = 2,
             caption = "Tabela 3. Eficiência relativa dos delineamentos (linhas) em relação a cada um dos delineamentos Bayesianos (colunas). Os valores de 100% na diagonal são a referência. Valores maiores indicam maior eficiência.")

Tabela 3. Eficiência relativa dos delineamentos (linhas) em relação a cada um dos delineamentos Bayesianos (colunas). Os valores de 100% na diagonal são a referência. Valores maiores indicam maior eficiência.
	U	UR	N	NR
U	100.00	102.61	99.98	102.89
UR	97.45	100.00	97.17	100.00
N	99.98	102.60	100.00	102.91
NR	97.45	100.00	97.17	100.00
RE	96.48	99.00	96.04	98.84
SS	96.44	98.96	95.86	98.65

#-----------------------------------------------------------------------
# Construção do gráfico.

# Pontos para desenhar a curva.
x <- seq(min(dom), max(dom), length.out = 51)
f <- function(x, theta) {
    flogis(x,
           ain = theta[1],
           asp = theta[2],
           xmd = theta[3],
           scl = theta[4])
}

# Simulando da distribuição conjunta a priori dos parâmetros.
set.seed(4321)
pri <- apply(prior_u$support,
             MARGIN = 2,
             FUN = function(v) {
                 runif(50, min(v), max(v))
             })

# Curvas simuladas.
cur <- cbind(
    expand.grid(x = x,
                id = seq_len(nrow(pri)),
                KEEP.OUT.ATTRS = FALSE),
    y = c(apply(pri, MARGIN = 1, FUN = f, x = x)))

i <- 0
labs <- c("Expaçamento regular",
          "Santos & Santos (2008)",
          "Bayesiano Uniforme",
          "Bayesiano Uniforme Restrito",
          "Bayesiano Gaussiano",
          "Bayesiano Gaussiano Restrito")
gg <- lapply(del_sup,
             FUN = function(sup) {
                 i <<- i + 1
                 ggplot(data = cur,
                        mapping = aes(x = x, y = y)) +
                     # geom_point() +
                     geom_line(mapping = aes(group = id)) +
                     geom_vline(xintercept = c(sup), linetype = 2) +
                     xlab("Idade dos animais (dias)") +
                     ylab("Peso do animal (kg)") +
                     ggtitle(labs[i]) +
                     theme(plot.title = element_text(size = 11))
             })

cap <- "Figura 1. Curvas simuladas a partir das priores uniformes. As linhas verticais tracejadas indicam os pontos de suporte de cada delineamento avaliado."
do.call(gridExtra::grid.arrange, gg)

Figura 1. Curvas simuladas a partir das priores uniformes. As linhas verticais tracejadas indicam os pontos de suporte de cada delineamento avaliado.

5 Discussão

Foram avaliados no total 6 delineamentos. Destes, 4 são delineamentos Bayesianos que foram obtidos usando-se o algorítmo da troca de coordenada aproximada, implementado no pacote acebayes. Além destes, o algorítmo sequencial de Santos & Santos (2008) foi usando. O último delineamento é o regularmente espaçado que corresponde a prática usual de ensaios para obtenção de curva de crescimento de aves de corte.

Os delineamentos estão disponíveis na Tabela 1 e são apresentados na Figura 1. Verifica-se que:

O delineamento de Santos & Santos (2008) resultou em uma disposição de pontos de suporte mais concentrada do lado direito, região do domínio onde a curva de crescimento apresenta maior taxa.
Os delineamentos Bayesianos uniforme e gaussiano ficaram bem próximos e apresentaram, do ponto de vista prático, apenas 4 pontos de suporte. Os três primeiros pontos de suporte foram duplicados. Arredondando-se para o inteiro mais próximo e tomando apenas os valores únicos, dos delineamentos são {1, 19, 33, 42} e {1, 18, 33, 42} para o uniforme e o gaussiano, respectivamente.
Os delineamentos Bayesianos uniforme e gaussiano com restrição de espaçamento mínimo de 5 dias foram, para a precisão considerada, idênticos nos pontos de suporte. Dessa forma, conclui-se que a restrição eliminou o efeito das priores da determinação dos pontos de suporte.

Com relação à eficiência relativa, conforme a Tabela 3, tem-se que:

Os delineamentos Bayesianos uniforme e Gaussiano foram os mais eficientes (a maior eficência é a de 100% nas respectivas colunas).
Os delineamentos Bayesianos uniforme e Gaussiano com restrição são menos eficientes que suas contrapartes não restritas. Isso mostra que a restrição no domínio reduziu a eficiência considerando o critério D-ótimo. Por outro lado, a eficiência não passou de 3% e dessa forma, talvez seja mais aconselhável, assumir na prática os delineamentos com restrição que, por terem mais pontos de suporte, permitem que modelos mais complexos sejam também ajustados aos dados.
A eficiência do delineamento de Santos & Santos (2008) foi inferior aos delineamentos com restrição, porém, a diferença não ultrapassou 3%. Dessa forma, dada a facilidade de obtenção de delieamentos com o algorítmo sequencial, considera-se para casos similares ao deste trabalho que é uma abordagem interessante.
O delineamento regularmente espaçado ficou muito próximo aos delineamentos Bayesianos com restrição. Esse fato era esperado já que a restição de separação mínima de 5 dias força um delineamento aproximadamente regular dentro de um intervalo de 42 dias.

Delineamentos ótimo bayesianos com o pacote `acebayes`

Luciana Helena Kowalski e Walmes Marques Zeviani

12 de Agosto de 2018

1 Instalação do pacote

2 Importação e análise descritiva dos dados

3 Ajuste da curva de crescimento

4 Obtenção dos delineamentos ótimos

4.1 Delimeanto ótimo Bayesiano com prioris uniformes

4.2 Delimeanto ótimo Bayesiano com prioris uniformes e restrição

4.3 Delimeanto ótimo Bayesiano com prioris Gaussianas

4.4 Delimeanto ótimo Bayesiano com prioris Gaussianas e restrição

4.5 Delineamento de Santos e Santos

4.6 Comparação dos resultados

5 Discussão

Delineamentos ótimo bayesianos com o pacote acebayes

Luciana Helena Kowalski e Walmes Marques Zeviani

12 de Agosto de 2018

1 Instalação do pacote

2 Importação e análise descritiva dos dados

3 Ajuste da curva de crescimento

4 Obtenção dos delineamentos ótimos

4.1 Delimeanto ótimo Bayesiano com prioris uniformes

4.2 Delimeanto ótimo Bayesiano com prioris uniformes e restrição

4.3 Delimeanto ótimo Bayesiano com prioris Gaussianas

4.4 Delimeanto ótimo Bayesiano com prioris Gaussianas e restrição

4.5 Delineamento de Santos e Santos

4.6 Comparação dos resultados

5 Discussão

Delineamentos ótimo bayesianos com o pacote `acebayes`