Setup

Instalando o pacote

#install.packages("INLA", repos=c(getOption("repos"), 
#                 INLA="https://inla.r-inla-download.org/R/stable"),
#                 dep=TRUE)

A instalação do pacote pode demorar algum tempo. arquivos de até 200Mb são baixados

Exemplo 1 - Sementes

A base em questão, “Seeds”, contém dados de um experimento fatorial 2 por 2. Cada linha representa uma planta. De cada planta foram retiradas n sementes e anotadas quantas germinaram .

Cada linha contém informação sobre as sementes vindas de um planta específica. As variáveis são referentes a informação de cada planta e são as seguintes:

#Carregando o pacote
library(INLA)
data(Seeds)  #Carregando dados no ambiente

dado <- data.frame(y = Seeds$r,
                   N = Seeds$n,
                   x1 = Seeds$x1,
                   x2 = Seeds$x2,
                   planta = Seeds$plate)

Análise Descritiva

summary(dado)
##        y               N               x1               x2             planta  
##  Min.   : 0.00   Min.   : 4.00   Min.   :0.0000   Min.   :0.0000   Min.   : 1  
##  1st Qu.: 8.00   1st Qu.:16.00   1st Qu.:0.0000   1st Qu.:0.0000   1st Qu.: 6  
##  Median :17.00   Median :39.00   Median :0.0000   Median :1.0000   Median :11  
##  Mean   :20.19   Mean   :39.57   Mean   :0.4762   Mean   :0.5238   Mean   :11  
##  3rd Qu.:26.00   3rd Qu.:51.00   3rd Qu.:1.0000   3rd Qu.:1.0000   3rd Qu.:16  
##  Max.   :55.00   Max.   :81.00   Max.   :1.0000   Max.   :1.0000   Max.   :21
plot(dado)

Modelo Proposto

\[ y_i \sim \text{Binomial}(N_i,p_i) \\ \eta_i = \text{logit}(p_i) = \beta_0 +\beta_1x_1 + \beta_2 x_2 + \epsilon_i \\ \epsilon_i \sim N(0, \sigma^2_\epsilon) \\ \beta_i \sim N(0,1000) \\ \pi(\sigma_\epsilon) = \lambda e^{-\lambda\sigma_\epsilon} \]

Foram escolhidas prioris normais para os \(\beta_i\) e exponencial para \(\sigma_\epsilon\) com \(\lambda = 4,6\) de forma que \(\pi(\sigma_\epsilon> 1) = 0,01\)

Especificando Modelo no pacote INLA

#Familia da verossimilhanca proposta
familia = "binomial"

#Função de ligação
fam.contr = list(control.link = list(model = "logit")) 


#Especificacao da priori como definido anteriormente
hyper = list(theta = list(prior = "pc.prec", param = c(1, 0.01)))

#mais informações sobre a priori pc para precisão:
#inla.doc("gev2")
              
#Estipulando equação do modelo
formula = y ~ x1 + x2 + f(planta, model = "iid", hyper = hyper)

O modelo é então ajustado passando-se os parâmetros especificados com a função inla()

res1 = inla(formula = formula, data = dado, 
            family = familia, Ntrials = N, 
            control.family = fam.contr)

Resultados do ajuste

summary(res1)
## 
## Call:
##    c("inla(formula = formula, family = familia, data = dado, Ntrials = N, 
##    ", " control.family = fam.contr)") 
## Time used:
##     Pre = 1.61, Running = 0.306, Post = 0.0623, Total = 1.98 
## Fixed effects:
##               mean    sd 0.025quant 0.5quant 0.975quant   mode kld
## (Intercept) -0.388 0.177     -0.733   -0.391     -0.027 -0.397   0
## x1          -0.354 0.226     -0.824   -0.347      0.073 -0.333   0
## x2           1.033 0.217      0.599    1.034      1.460  1.037   0
## 
## Random effects:
##   Name     Model
##     planta IID model
## 
## Model hyperparameters:
##                        mean      sd 0.025quant 0.5quant 0.975quant mode
## Precision for planta 111.07 1231.38       2.98    10.91     304.93 6.06
## 
## Expected number of effective parameters(stdev): 9.91(2.94)
## Number of equivalent replicates : 2.12 
## 
## Marginal log-Likelihood:  -68.44

Efeitos aleatórios:

res1$summary.random$planta
##    ID        mean        sd  0.025quant    0.5quant 0.975quant         mode
## 1   1 -0.28720994 0.2737350 -0.89550438 -0.25826247 0.16905303 -0.186445632
## 2   2 -0.08036863 0.2262876 -0.55905845 -0.06824044 0.35281483 -0.029946894
## 3   3 -0.31598605 0.2438226 -0.84400359 -0.29632804 0.09765073 -0.256962226
## 4   4  0.21681476 0.2400774 -0.21405414  0.19842004 0.73270009  0.145171878
## 5   5  0.05541629 0.2394591 -0.42072181  0.04861295 0.54975521  0.022742130
## 6   6  0.09820756 0.3123291 -0.48952622  0.07352890 0.78702693  0.019751308
## 7   7  0.15747742 0.2269242 -0.26414891  0.14217716 0.64103503  0.084553847
## 8   8  0.28362392 0.2440625 -0.13632874  0.26358616 0.81257666  0.219213093
## 9   9 -0.06126372 0.2324867 -0.54650688 -0.05229047 0.39354518 -0.023175363
## 10 10 -0.18719170 0.2245423 -0.66716850 -0.17107967 0.22165644 -0.120961645
## 11 11  0.11456874 0.2912521 -0.42793127  0.09021962 0.75478818  0.029011956
## 12 12  0.20648261 0.2986754 -0.31196508  0.17182697 0.87865841  0.061445945
## 13 13  0.02079936 0.2577052 -0.49423949  0.01573915 0.55533009  0.004298994
## 14 14 -0.05913412 0.2614629 -0.60758429 -0.04922351 0.45626288 -0.019586001
## 15 15  0.38726559 0.2891710 -0.08563289  0.36018695 1.02171430  0.305109145
## 16 16 -0.12615612 0.3266584 -0.86102883 -0.09422489 0.47144360 -0.026139129
## 17 17 -0.29097305 0.3224983 -1.02822765 -0.24842003 0.23464281 -0.102633381
## 18 18 -0.05976517 0.2420373 -0.55643584 -0.05383401 0.42546209 -0.026978992
## 19 19 -0.10987374 0.2545429 -0.64759864 -0.09560701 0.38135814 -0.042119284
## 20 20  0.12460890 0.2436790 -0.32492687  0.10532841 0.65278838  0.043287380
## 21 21 -0.08474342 0.3024644 -0.73868669 -0.06589786 0.49904458 -0.019675619
##             kld
## 1  1.090805e-04
## 2  4.802824e-05
## 3  2.187782e-04
## 4  1.047552e-04
## 5  6.100780e-05
## 6  3.979882e-05
## 7  5.926546e-05
## 8  1.602085e-04
## 9  6.041292e-05
## 10 9.532587e-05
## 11 4.556427e-05
## 12 7.053148e-05
## 13 5.515244e-05
## 14 5.193019e-05
## 15 2.894837e-04
## 16 3.136087e-05
## 17 6.867367e-05
## 18 6.044082e-05
## 19 6.592566e-05
## 20 4.562849e-05
## 21 5.363477e-05

Quantis dos efeitos aleatórios:

plot(1:nrow(dado), res1$summary.random$planta$`0.97`, col="red", ylim=c(-1,1),
     xlab="Identificador da Planta", ylab = "Quantis")
points(1:nrow(dado), res1$summary.random$planta$`0.5quant`)
points(1:nrow(dado), res1$summary.random$planta$`0.02`, col="blue")

Marginal de um efeito fixo:

m.beta1 = inla.tmarginal(fun = function(x) x, 
                         marginal = res1$marginals.fixed$x1)
#Coeficiente de x1
#Aqui pode se ver uma transformação do parâmetro, no nosso caso a identidade.

plot(m.beta1, type="l", xlab = expression(beta[1]), ylab = "Densidade de Probabilidade")

Marginal do Hiperparametro \(\sigma\)

m.sigma = inla.tmarginal(fun = function(x) exp(-1/2*x), 
                         marginal = res1$internal.marginals.hyperpar$`Log precision for planta`)
#m.sigma é a marginal do desvio padrão no efeito aleatório iid.
plot(m.sigma, type="l", xlab = expression(sigma[iid]), ylab = "Densidade de Probabilidade")

Marginal de um parâmetro

m.plate.7 = inla.tmarginal(fun = function(x) x, 
                           marginal = res1$marginals.random$planta$index.7)

# - m.plate.7 is one of the marginals for the parameters beloning to the plate iid effect
# - it is number 7, which corresponds to plate=7, which is our 7th row of data
plot(m.plate.7, type="l", xlab = "Marginal do Efeito da Sétima Planta", ylab = "Densidade de Probabilidade")

Distribuição das estimativas do efeito aleatório

plot(density(res1$summary.random$planta$mean), main = "Curva Suavizada das Medianas dos efeitos aleatórios")

Previsão

dado2 = rbind(dado, c(NA, 1, 0, 0, 22))
tail(dado2)
##     y  N x1 x2 planta
## 17  3 12  1  1     17
## 18 22 41  1  1     18
## 19 15 30  1  1     19
## 20 32 51  1  1     20
## 21  3  7  1  1     21
## 22 NA  1  0  0     22
res.pred = inla(formula = formula, data=dado2, 
            family = familia, Ntrials = N, 
            control.family = fam.contr, 
            control.predictor = list(compute=T, link = 1))# - to get the posterior of the predictor and fitted values


res.pred$summary.fitted.values[22, ]
##                          mean         sd 0.025quant  0.5quant 0.975quant
## fitted.Predictor.22 0.4070652 0.08791902  0.2402311 0.4025713   0.599076
##                          mode
## fitted.Predictor.22 0.3963464

Exemplo 2 - Pacientes

A base em questão se trata de dados sobre pacientes monitorados durante oito semanas sob um tratamento para crises epiléticas. As quantidades de crises foi medida em 4 tempos diferentes. As variáveis da base são:

O modelo proposto é o seguinte:

\[ y_{jk} \sim poisson(m_{jk}) \\ \text{log}(m_{jk}) = a_0 + a_{base}\text{log}(\text{Base}_j/4) + a_{trt}\text{Trt} + a_{BT}\text{Trt} \text{log}(\text{Base}_j/4) + a_{Age}\text{Age}_j + a_{V4}V4 + b_{1j} + b_{jk} \\ b_{1j} \sim N(0, t_{b_1}) \\ b_{jk} \sim N(0, t_{b}) \]

Foram dadas prioris vagas com distribuição normal para as \(a_i\) e gama para as quanridade \(1/\sqrt{t_{b1}}\) e \(1/\sqrt{t_{b}}\).

data(Epil)
head(Epil)
##   y Trt Base Age V4 rand Ind
## 1 5   0   11  31  0    1   1
## 2 3   0   11  31  0    2   1
## 3 3   0   11  31  0    3   1
## 4 3   0   11  31  1    4   1
## 5 3   0   11  30  0    5   2
## 6 5   0   11  30  0    6   2
centralizar = function(x) (x - mean(x))

##Centralizando covariaveis
Epil$CTrt    <- centralizar(Epil$Trt)
Epil$ClBase4 <- centralizar(log(Epil$Base/4))
Epil$CV4     <- centralizar(Epil$V4)
Epil$ClAge   <- centralizar(log(Epil$Age))
Epil$CBT     <- centralizar(Epil$Trt*Epil$ClBase4)

formula = y ~ ClBase4 + CTrt + CBT+ ClAge + CV4 +
          f(Ind, model="iid") + f(rand,model="iid")

result = inla(formula,family="poisson", data = Epil)
result$summary.fixed
##                   mean         sd  0.025quant   0.5quant  0.975quant       mode
## (Intercept)  1.5771961 0.07580877  1.42553006  1.5780121  1.72422579  1.5796293
## ClBase4      0.8793213 0.13395839  0.61517940  0.8793189  1.14314699  0.8793220
## CTrt        -0.3324507 0.15089200 -0.63172558 -0.3318023 -0.03711348 -0.3305241
## CBT          0.3502933 0.20719833 -0.05753924  0.3500886  0.75877250  0.3497193
## ClAge        0.4831725 0.35391801 -0.21790398  0.4843698  1.17686219  0.4867573
## CV4         -0.1032028 0.08516969 -0.27051014 -0.1032553  0.06420664 -0.1033527
##                      kld
## (Intercept) 2.938955e-06
## ClBase4     1.729027e-06
## CTrt        1.460274e-06
## CBT         1.580862e-06
## ClAge       1.329954e-06
## CV4         1.118419e-06
result$summary.random$Ind[1:5,]
##   ID       mean        sd 0.025quant   0.5quant 0.975quant       mode
## 1  1 0.04450906 0.2862575 -0.5265456 0.04727954  0.5995689 0.05265406
## 2  2 0.05534930 0.2858835 -0.5149913 0.05813135  0.6096276 0.06352566
## 3  3 0.28032835 0.3223808 -0.3574378 0.28139470  0.9111966 0.28327256
## 4  4 0.12849712 0.3056776 -0.4778603 0.13012234  0.7250564 0.13321291
## 5  5 0.01474608 0.2588728 -0.4953460 0.01490969  0.5233215 0.01526253
##            kld
## 1 2.856627e-08
## 2 3.063263e-08
## 3 5.940168e-07
## 4 2.274010e-07
## 5 9.606335e-07
#Posteriori para efeitos aleatórios dos 5 primeiros individuos

result$summary.random$rand[1:5,]
##   ID         mean        sd 0.025quant     0.5quant 0.975quant         mode
## 1  1  0.131889713 0.3012600 -0.4626794  0.132134233  0.7237947  0.132143663
## 2  2 -0.037800566 0.3053294 -0.6480746 -0.034880827  0.5553634 -0.029514042
## 3  3 -0.037800566 0.3053294 -0.6480746 -0.034880827  0.5553634 -0.029514042
## 4  4 -0.009395635 0.3077505 -0.6234848 -0.006844794  0.5895505 -0.002222897
## 5  5 -0.036316782 0.3054311 -0.6467071 -0.033429251  0.5571363 -0.028127359
##            kld
## 1 8.107996e-08
## 2 9.503580e-08
## 3 9.503580e-08
## 4 1.172704e-07
## 5 8.786394e-08
#plot(result, main = NULL)

a = inla.tmarginal(fun = function(x) exp(-1/2*x), 
                         marginal = result$marginals.fixed$ClBase4)
plot(a, main = "Aproximação para posteriori do coeficiente ClBase4")