4 Modelo beta-binomial

Assumindo o modelo binomial ou de Poisson, assume-se que a probabilidade de incidência no talhão é constante em todo o talhão, porém no caso de haver agregação essa suposição não é razoável e o índice de dispersão é maior que um. Supondo que a probabilidade de incidência é variável, pode-se modelar os dados de quadrat counts por uma distribuição beta-binomial.

A distribuição beta-binomial é uma composição da distribuição binomial com a densidade beta, (Skellam 1948). A distribuição beta-binomial tem dois parâmetros: π, que representa a incidência média no talhão, e θ que mede a variação da incidência. O parâmetro θ, no contexto dos dados de incidência de doenças, pode ser interpretado como um índice de agregação espacial da doença. Se θ > 0, o padrão da doença é agregado (Griffiths 1979).

As estimativas de máxima verossimilhança para a incidência média p e o parâmetro de agregação θ do modelo beta-binomial não tem expressões fechadas. Para obter essas estimativas, utiliza-se um procedimento iterativo de minimização (maximização) numérica, que pode não convergir.


  > disp.quadrats(v303.geo, dx = 2, dy = 5, death = 1:3, model = "beta-binomial")


  $2x5
         N    n  nN   prob  theta p.value  pattern
  Av1  100 9.45 945 0.0238 0.0986 0.01229 Agregate
  Av2  100 9.45 945 0.0283 0.0692 0.04202 Agregate
  Av3  100 9.45 945 0.0778 0.0695 0.01948 Agregate
  Av4  100 9.45 945 0.0948 0.0993 0.00174 Agregate
  Av5  100 9.45 945 0.1041 0.0781 0.00829 Agregate
  Av6  100 9.45 945 0.1570 0.0984 0.00147 Agregate
  Av7  100 9.45 945 0.2520 0.0829 0.00392 Agregate
  Av8  100 9.45 945 0.2568 0.0700 0.01075 Agregate
  Av9  100 9.45 945 0.3175 0.0617 0.01882 Agregate
  Av10 100 9.45 945 0.3357 0.0722 0.00781 Agregate
  Av11 100 9.45 945 0.4925 0.0936 0.00124 Agregate
  Av12 100 9.45 945 0.5522 0.1153 0.00016 Agregate
  Av13 100 9.45 945 0.6086 0.1602 0.00000 Agregate
  Av14 100 9.45 945 0.6210 0.1633 0.00000 Agregate
  Av15 100 9.45 945 0.6530 0.1822 0.00000 Agregate
  Av16 100 9.45 945 0.6714 0.1892 0.00000 Agregate
  Av17 100 9.45 945 0.7616 0.1175 0.00013 Agregate
  Av18 100 9.45 945 0.8219 0.1190 0.00016 Agregate
  Av19 100 9.45 945 0.8262 0.1202 0.00015 Agregate
  Av20 100 9.45 945 0.8730 0.1045 0.00097 Agregate
  Av21 100 9.45 945 0.9272 0.1406 0.00018 Agregate
  Av22 100 9.45 945 0.9272 0.1406 0.00018 Agregate
  Av23 100 9.45 945 0.9272 0.1406 0.00018 Agregate
  Av24 100 9.45 945 0.9711 0.1058 0.00730 Agregate
  Av25 100 9.45 945 0.9711 0.1058 0.00730 Agregate

4.1 Detalhes da estimação

O procedimento de estimação dos parâmetro da distribuição beta-binomial, esta implementado em detalhes na função betabinom.citrus().


  > args(betabinom.citrus)


  function (data, dx, dy = dx, check.args = NULL, ini.prob = NULL,
      ini.theta = NULL, usage = c("fitdistr", "mle"), ...)
  NULL

Nos argumentos dessa função precisamos entrar com os dados no argumento data, as dimensões dos quadrats nos argumentos dx e dy e as informações dos códigos identificadores do status das plantas em forma de lista nomeada.

Duas funções podem ser utilizadas na minimização numérica: a fitdistr() do pacote MASS, (Venables & Ripley 2002) e a função mle() do pacote stats4. A primeira, retorna simplesmente as estimativas dos parâmetros e dos erros-padrão obtidos do hessiano numérico A segunda, retorna um resultado que possibilita explorar as verossimilhanças perfilhadas. A escolha dessa função é informada no argumento usage.

Nos dados de ppc podemos utilizar a função fitdistr, fazendo:


  > (bb1ita <- betabinom.citrus(gItajobi, dx = 10))


  valid evaluations: 1
      prob      theta
    0.00638   0.01934
   (0.00236) (0.01330)

4.2 Verossimilhanças perfilhadas

Utilizando usage=’mle’, retorna-se um objeto da classe mle. Com esse resultado, podemos estudar as verossimilhanças perfilhadas, para aproximar intervalos de confiança para as estimativas e visualizar.


  > bb.ita <- betabinom.citrus(gItajobi, dx = 10, usage = "mle")


  valid evaluations: 1

Obtendo o sumário:


  > summary(bb.ita)


  Maximum likelihood estimation
  
  Call:
  mle(minuslogl = function(prob, theta) -sum(dbetabinom(x = y,
      size = size, prob = prob, theta = theta, log = TRUE)), start = list(prob = ini.prob,
      theta = ini.theta))
  
  Coefficients:
        Estimate Std. Error
  prob   0.00638    0.00236
  theta  0.01934    0.01330
  
  -2 log L: 62.6

Para obter o intervalo de confiança para as estimativas a partir das verossimilhanças perfilhadas, faz-se:


  > prof <- profile(bb.ita)
  > confint(prof)


          2.5 % 97.5 %
  prob  0.00288 0.0146
  theta 0.00463 0.0890

As verossimilhanças perfilhadas podem ser plotadas, como na Figura 2.


  > par(mfrow = c(1, 2), mar = c(3, 3, 3, 0.5), mgp = c(2, 1,
  +     0))
  > plot(prof, absVal = FALSE)


pict

Figura 2: Visualização da verossimilhança perfilhada para os parâmetros da distribuição beta-binomial aplicada aos dados de Itajobi