Ajuste de modelos de regressão censurados em diferentes pacotes estatísticas

Iniciação Científica 2018-2019

Programa de Pesquisa para Alunos da Graduação

Aluno: Rafael Cabral Fernandez
Orientadores: Gustavo H. M. A. Rocha e Renata S. Bueno

Introdução

O documento a seguir apresenta a utilização do BUGS, utilizando a versão WinBugs.

Para utilização do programa, é necessário baixar o programa em: https://www.mrc-bsu.cam.ac.uk/software/bugs/the-bugs-project-winbugs/

A instalão do programa deve ser feito em uma pasta de fácil localização.

Para utilização interna, basta utilizar o pacote “R2WinBugs”.

O pacote foi rodado em 4 exemplos diferentes:

  1. Exemplo didático para treino
  2. Simulação de um modelo de regressão linear simples
  3. Modelo de regressão linear simples utilizando dados do ipea, modelando a taxa de homicídos através da taxa de desemprego
  4. Modelo de regressão linear simples utilizando dados do enem, modelando a nota de matemática através da nota de redação

Treino

Será utilizado o conjunto de dados Iris

Modelo de regressão linear simples

data(iris)
summary(iris)
##   Sepal.Length    Sepal.Width     Petal.Length    Petal.Width   
##  Min.   :4.300   Min.   :2.000   Min.   :1.000   Min.   :0.100  
##  1st Qu.:5.100   1st Qu.:2.800   1st Qu.:1.600   1st Qu.:0.300  
##  Median :5.800   Median :3.000   Median :4.350   Median :1.300  
##  Mean   :5.843   Mean   :3.057   Mean   :3.758   Mean   :1.199  
##  3rd Qu.:6.400   3rd Qu.:3.300   3rd Qu.:5.100   3rd Qu.:1.800  
##  Max.   :7.900   Max.   :4.400   Max.   :6.900   Max.   :2.500  
##        Species  
##  setosa    :50  
##  versicolor:50  
##  virginica :50  
##                 
##                 
## 
plot(iris$Sepal.Length ~ iris$Petal.Length, xlab = "Petal length", ylab="Sepal length")

sepal <- iris$Sepal.Length
petal <- iris$Petal.Length
lm1 <- lm(sepal ~ petal)
summary(lm1)
## 
## Call:
## lm(formula = sepal ~ petal)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.24675 -0.29657 -0.01515  0.27676  1.00269 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  4.30660    0.07839   54.94   <2e-16 ***
## petal        0.40892    0.01889   21.65   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.4071 on 148 degrees of freedom
## Multiple R-squared:   0.76,  Adjusted R-squared:  0.7583 
## F-statistic: 468.6 on 1 and 148 DF,  p-value: < 2.2e-16


Código para efetuar o MCMC, deixando em comentário pela falta de compatibilidade entre o RMarkdown e o WinBugs

sink("mod1.txt")        
cat("
MODEL LR1 {
    for(i in 1:N) {
    sepal[i] ~ dnorm(mu[i], tau)
    mu[i] <- alpha + beta*petal[i]
    }
    #...#...
    alpha ~ dnorm(0,0.001)
    beta ~ dnorm(0,0.001)
    
    tau <- pow(sigma, -2)
    sigma ~ dunif(0,10) 
    #...#...
} ", fill = TRUE)
sink()


N = length(sepal)
data = list("N","sepal","petal")
params = c("alpha", "beta")  

inits <- function () {list(alpha = rnorm(1),
                           beta = rnorm(1),
                           sigma = 1)}

nc <- 3      #number of MCMC chains to run
ni <- 5000  #number of samples for each chain     
nb <- 1000   #number of samples to discard as burn-in
nt <- 1      #thinning rate, increase this to reduce autocorrelation



result <- bugs(data = data, inits = inits, parameters.to.save = params,
                model.file = "mod1.txt", n.chains = 1, n.iter = 10000,
                n.burnin = 2000, bugs.directory = "C:\\Program Files\\WinBUGS14",
                debug = FALSE, save.history = FALSE, DIC = FALSE)


Contudo, é possível salvar os resultados e obter algumas informações resumo de forma simples. (Como não há compatilidade, não há como imprimir os resultados pelo r markdown diretamente)

#result
#result$summary

load("C:/Users/Rafael/Documents/result.RData")

result
## Inference for Bugs model at "mod1.txt", fit using WinBUGS,
##  1 chains, each with 10000 iterations (first 2000 discarded), n.thin = 8
##  n.sims = 1000 iterations saved
##       mean  sd 2.5% 25% 50% 75% 97.5%
## alpha  4.3 0.1  4.2 4.3 4.3 4.4   4.5
## beta   0.4 0.0  0.4 0.4 0.4 0.4   0.4
kable(result$summary)
mean sd 2.5% 25% 50% 75% 97.5%
alpha 4.3069970 0.0784842 4.1568499 4.256 4.30850 4.35925 4.4600000
beta 0.4092224 0.0189578 0.3726975 0.397 0.40875 0.42260 0.4456225 
kable(head(result$sims.matrix))
alpha beta
4.393 0.3793
4.248 0.4269
4.177 0.4297
4.320 0.4131
4.370 0.3867
4.401 0.4032

Onde alpha representa o intercepto, e beta o coeficiente de inclinação.

Simulação

Gerando os dados

phi_verdadeiro=5
beta0_verdadeiro=-2
beta1_verdadeiro=2
n=100
x=rnorm(n,0,1)
e=rnorm(n,0,sqrt(phi_verdadeiro))
y=beta0_verdadeiro+beta1_verdadeiro*x+e
plot(x,y)


Código para o modelo WinBugs

sink("mod1.txt")        
cat("
    MODEL LR1 {
    for(i in 1:N) {
    y[i] ~ dnorm(mu[i], tau)
    mu[i] <- beta0 + beta1*x[i]
    }
    beta0 ~ dnorm(0,0.001)
    beta1 ~ dnorm(0,0.001)
    sigma ~ dgamma(3,10)  
    tau <- pow(sigma, -2)
    } ", fill = TRUE)
sink()

N = length(y)
data = list("N","y","x")
params = c("beta0", "beta1","tau")  

inits <- function () {list(beta0 = rnorm(1),
                           beta1 = rnorm(1),
                           tau = 1)}


nc <- 3      #number of MCMC chains to run
ni <- 10000  #number of samples for each chain     
nb <- 1000   #number of samples to discard as burn-in
nt <- 1      #thinning rate, increase this to reduce autocorrelation

result <- bugs(data = data, inits = inits, parameters.to.save = params,
               model.file = "mod1.txt", n.chains = 1, n.iter = 10000,
               n.burnin = 2000, bugs.directory = "C:\\Program Files\\WinBUGS14",
               debug = TRUE, save.history = TRUE, DIC = FALSE)

Resultados

load("C:/Users/Rafael/Documents/result2.RData")

result
## Inference for Bugs model at "mod1.txt", fit using WinBUGS,
##  1 chains, each with 10000 iterations (first 2000 discarded), n.thin = 8
##  n.sims = 1000 iterations saved
##       mean  sd 2.5%  25%  50%  75% 97.5%
## beta0 -1.9 0.2 -2.3 -2.0 -1.9 -1.7  -1.4
## beta1  1.8 0.2  1.4  1.7  1.8  1.9   2.2
## tau    0.2 0.0  0.2  0.2  0.2  0.2   0.3
kable(result$summary)
mean sd 2.5% 25% 50% 75% 97.5%
beta0 -1.874223 0.2173828 -2.3070250 -2.01300 -1.87100 -1.7280 -1.440850
beta1 1.807659 0.2081142 1.3909500 1.67375 1.80750 1.9410 2.216075
tau 0.213182 0.0260787 0.1661925 0.19500 0.21215 0.2294 0.267605 
kable(head(result$sims.matrix))
beta0 beta1 tau
-1.658 2.015 0.2345
-1.615 1.832 0.1868
-1.893 1.875 0.2557
-1.953 1.666 0.2001
-2.247 1.839 0.2128
-1.920 1.620 0.1822


Homícidio

Ajustando os dados

y = homicidio$X__2
x = homicidio$X__1
lm1 <- lm(y ~ x)
summary(lm1)

Código para o modelo WinBugs

sink("mod1.txt")        
cat("
    MODEL LR1 {
    for(i in 1:N) {
    y[i] ~ dnorm(mu[i], tau)
    mu[i] <- beta0 + beta1*x[i]
    }
    beta0 ~ dnorm(0,0.001)
    beta1 ~ dnorm(0,0.001)
    sigma ~ dgamma(3,10)  
    tau <- pow(sigma, -2)
    } ", fill = TRUE)
sink()



N = length(y)
data = list("N","y","x")
params = c("beta0", "beta1","tau")  

inits <- function () {list(beta0 = rnorm(1),
                           beta1 = rnorm(1),
                           tau = 1)}


nc <- 3      #number of MCMC chains to run
ni <- 10000  #number of samples for each chain     
nb <- 1000   #number of samples to discard as burn-in
nt <- 1      #thinning rate, increase this to reduce autocorrelation

result <- bugs(data = data, inits = inits, parameters.to.save = params,
               model.file = "mod1.txt", n.chains = 1, n.iter = 10000,
               n.burnin = 2000, bugs.directory = "C:\\Program Files\\WinBUGS14",
               debug = TRUE, save.history = TRUE, DIC = FALSE)

Resultados

load("C:/Users/Rafael/Documents/result3.RData")

result
## Inference for Bugs model at "mod1.txt", fit using WinBUGS,
##  1 chains, each with 10000 iterations (first 2000 discarded), n.thin = 8
##  n.sims = 1000 iterations saved
##       mean  sd 2.5%  25%  50%  75% 97.5%
## beta0 15.7 1.8 12.0 14.4 15.8 16.9  19.3
## beta1  1.1 0.2  0.8  1.0  1.1  1.3   1.5
## tau    1.1 0.3  0.6  0.9  1.1  1.3   1.8
kable(result$summary)
mean sd 2.5% 25% 50% 75% 97.5%
beta0 15.668716 1.8361596 12.0179935 14.39000 15.780 16.91250 19.280250
beta1 1.136014 0.1992134 0.7580249 0.99250 1.125 1.27625 1.548125
tau 1.086771 0.3130508 0.5699574 0.85335 1.068 1.26900 1.792025 
kable(head(result$sims.matrix))
beta0 beta1 tau
16.01 1.108 1.0760
17.08 0.997 0.8784
19.81 0.686 0.9130
15.85 1.114 1.0680
11.61 1.573 0.6929
16.72 1.044 0.8514


Enem

Carregando os dados

y = enem$NU_NOTA_REDACAO
x = enem$NU_NOTA_MT
lm1 <- lm(y ~ x)
summary(lm1)


Código para o modelo WinBugs

sink("mod1.txt")        
cat("
    MODEL LR1 {
    for(i in 1:N) {
    y[i] ~ dnorm(mu[i], tau)
    mu[i] <- beta0 + beta1*x[i]
    }
    beta0 ~ dnorm(0,0.001)
    beta1 ~ dnorm(0,0.001)
    sigma ~ dgamma(0.1,0.001)  
    tau <- pow(sigma, -2)
    } ", fill = TRUE)
sink()

Infelizmente, o winbugs não foi capaz de rodar este modelo.