CE-227 - Primeiro semestre de 2014

No quadro abaixo será anotado o conteúdo dado em cada aula do curso.
São indicados os Capítulos e Sessões correspondentes nas referências bibliográficas, bem como os exercícios sugeridos.

Veja ainda depois da tabela as Atividades Complementares.

Observação sobre exercícios recomendados os exercícios indicados são compatíveis com o nível e conteúdo do curso.
Se não puder fazer todos, escolha alguns entre os indicados.

• Procurar materiais introdutórios sobre inferência Bayesiana na web, livros, etc
• Escrever um programa que encontre os parâmetros da distribuição Beta a partir de estimativas pontuais, intervalos e respectiva probabilidade (como visto na aula)
• Investigar como foram encontrados os parâmetros na distribuição que combina dados e informações subjetivas mostrada em aula.

Seguem alguns resultados obtidos com meu código para encontrar os parâmetros da priori beta a partir de um valor que considera-se ser a moda, uma faixa de valores e uma probabilidade associada a esta faixa

1. Considere o exemplo visto em aula (lembrando que consideramos n=80 e y=24) e a sua escolha de priori e mostre como obter resumos pontuais da posteriori do modelo Binomial-Beta de duas formas distintas:
   • usando resultados analíticos da distribuição,
   • através de algum algorítmo computacional.
2. Mostrar como obter intervalos de quantis e HPD para posteriori do modelo Binomial-Beta. Escrever rotinas computacionais e obter resultados para o exemplo de aula.
3. Repetir anteriores para o modelo Poisson-Gama

1. Obter amostras das posterioris nos modelo Binomial-Beta e Poisson-Gama. Extrair resultados e comparar com os analíticos obtidos anteriormente.
2. Obter a aproximação de Normal para as posterioris destes dois modelos
3. Escrever um algoritmo MCMC para obter amostras destes dois modelos (supondo - artificialmente - que a posteriori não é de forma conhecida).

1. Refazer a questão da 1a avaliação.
2. Encontrar a aproximação de Normal para a posteriori da questão da 1a avaliação. Fazer gráfico sobrepondo as distribuições, comparar resumos descritivos das distribuições. Avaliar a aproximação.
3. Reconsidere a questão da 1a avaliação utilizando o parâmetro de precisão . Encontre a priori correspondente à definida (por transformação de variáveis).
4. Mostre que para questão da 1a avaliação a distribuição gama inversa é conjugada. Idem para distribuição gama para o parâmetro de precisão.
5. Considere agora a distribuição indexada para parâmetro (reparametrização) . Encontre a priori equivalente na definida na avaliação, encontre a distribuição a posteriori e a aproximação da Laplace correspondente. Sobreponha em um gráfico e discuta.
6. Monte um algoritmo MCMC para amostrar de algumas das posterioris da avaliação ou itens anteriores. Sobreponha em um gráfico a densidade obtida e a estimada pelas amostras MCMC. Transforme os valores simulados de acordo com as transformações dos itens anteriores e novamente compare graficamente as posterioris obtidas analiticamente e pelas simulações.

1. Obter a aproximação normal para os exemplos vistos até agora no curso (modelos Poisson, Geométrico, Binomial Negativo, Binomial, Expoencial, Gamma, …)
2. Montar um algorítmo MCMC para os exemplos vistos até agora no curso (modelos Poisson, Geométrico, Binomial Negativo, Binomial, Exponencial, Gamma, …)

1. Considerar o problema da 3a avaliação semanal. Resolver analicamente e por métodos numéricos.

1. Considere o modelo binomial e uma amostra com n=15 e y=6. Obtenha os gráficos da "trinca" priori, verossimilhança e posteriori nos seguintes casos:
   • 
   • 
   • : priori de Jeffrey's

1. Implementar algorítmos de inferência para o problema da 3a avaliação semanal. Fazer gráficos das quantidades e funções de interesse.

1. Obter as expressões analíticas das posterioris conjunta, marginais e condicionais

do modelo de regressão linear com as prioris de referência.

1. Estender o código a seguir (visto em aula) para comparar os resultados analíticos com os obtidos pelo amostrador de Gibbs. Incluir na comparação as distribuições e resumos pontuais e intervalares.
2. Generalizar os resultados analíticos e códigos para priori Normal-InversaGamma

## simulando dados
x <- sort(runif(20, 0, 20))
y <- 2 + 0.25*x + rnorm(20, m=0, sd=2)
## visualizando dados e reta de regressão "usual"
plot(x,y)
reg <- lm(y ~ x)
abline(reg)
## Código ("ingênuo") para amostrador de Gibbs
reg.lin.GIBBS <- function(y, x, Nsim, iniBeta){ 
    require(MASS)
    n <- length(y)
    X <- cbind(1, x)
    XX <- crossprod(X)
    bhat <- solve(XX,crossprod(X,y))
    ##
    sim <- matrix(0, nrow = Nsim, ncol=3)
    sim[1,1:2] <- iniBeta
    sim[1, 3]  <- 1/rgamma(1, shape=(n-2)/2, scale=2/crossprod(y-X%*%sim[1,1:2]))
    for(i in 2:Nsim){
        sim[i, 3]  <- 1/rgamma(1, shape=(n-2)/2, scale=2/crossprod(y-X%*%sim[(i-1),1:2]))
        sim[i,1:2] <- mvrnorm(n=1, mu=bhat, Sigma = solve(XX)*sim[i, 3])
    }
    return(sim)
}
## Obtendo amostras
rlG0 <- reg.lin.GIBBS(y=y, x=x, Nsim=15000, iniBeta=c(0,0.6))
## Burn-in (3000) e thining (1/10)
rlG1 <- rlG0[-(1:3000),]
rlG2 <- rlG1[10*(1:1200),]
 
## comparando estimativas "usuais" e médias das posterioris
c(coef(reg), summary(reg)$sigma^2)
apply(rlG2, 2, mean)
 
## traços e posterioris (por simulação) com indicação das estimativas "usuais"
par(mfrow=c(2,3))
plot(rlG2[,1], type="l")
plot(rlG2[,2], type="l")
plot(rlG2[,3], type="l")
plot(density(rlG2[,1])); abline(v=coef(reg)[1])
plot(density(rlG2[,2])); abline(v=coef(reg)[2])
plot(density(rlG2[,3])); abline(v=summary(reg)$sigma^2)

1. Estender os códigos anteriores para incluir a predição de novos valores

#### código agora incluindo predições
reg.lin.GIBBS <- function(y, x, Nsim, iniBeta, x.pred = NULL){	
    require(MASS)
    n <- length(y)
    X <- cbind(1, x)
    XX <- crossprod(X)
    bhat <- solve(XX,crossprod(X,y))
    ##
    NC <- ncol(X) + 1
    if(!is.null(x.pred)){
        Xpred <- cbind(1, x.pred)
        NC <- NC + nrow(Xpred) 
    }
    sim <- matrix(0, nrow = Nsim, ncol=NC)
    sim[1,1:2] <- iniBeta
    sim[1, 3]  <- 1/rgamma(1, shape=(n-2)/2, scale=2/crossprod(y-X%*%sim[1,1:2]))
    ##
    for(i in 2:Nsim){
        sim[i, 3]  <- 1/rgamma(1, shape=(n-2)/2, scale=2/crossprod(y-X%*%sim[(i-1),1:2]))
        sim[i,1:2] <- mvrnorm(n=1, mu=bhat, Sigma = solve(XX)*sim[i, 3])
        if(!is.null(x.pred))
            sim[i,-(1:3)] <- rnorm(nrow(Xpred), m=Xpred%*%sim[i,1:2], sd=sqrt(sim[i,3]))
    }
    return(sim)
}
## agora com predição
rlG0 <- reg.lin.GIBBS(y=y, x=x, Nsim=15000, iniBeta=c(0,0.6), x.pred=1:20)
dim(rlG0)
rlG1 <- rlG0[-(1:3000),]
dim(rlG1)
rlG2 <- rlG1[10*(1:1200),]
dim(rlG2)
rbind(dim(rlG0), dim(rlG1), dim(rlG2))
 
reg
summary(reg)$sigma^2
 
par(mfrow=c(4,5))
apply(rlG2[,-(1:3)], 2, function(x) plot(density(x),type="l"))
par(mfrow=c(1,1))
 
y.pred <- apply(rlG2[,-(1:3)], 2, mean)
plot(y ~ x)
abline(reg)
lines(1:20, y.pred, col=2)
#
y.pred <- apply(rlG1[,-(1:3)], 2, mean)
plot(y ~ x)
abline(reg)
lines(1:20, y.pred, col=2)

• Incluir bandas de predição usuais e bayesianas no gráfico
• Generalizar código para família de priori Normal-GammaInversa
• Verificar o efeito/sensitividade à espeficicação das prioris

Exemplos JAGS/rjags

1. Média e variância (precisão) da distribuição normal

   n <- 20
   x <- rnorm(n, 70, 5)
    
   #write.table(x,
   #            file = 'normal.data',
   #            row.names = FALSE,
   #            col.names = FALSE)
    
   cat( "model {
   	for (i in 1:n){
   		x[i] ~ dnorm(mu, tau)
   	}
   	mu ~ dnorm(0, 0.0001)
   	tau <- pow(sigma, -2)
   	sigma ~ dunif(0, 100)
    }", file="normal.modelo"
   )
    
    
   require(rjags)
    
   ## OBS: valores iniciais são dispensáveis neste exemplo
   inis <- list(list(mu=10, sigma=2),
                list(mu=50, sigma=5),
                list(mu=70, sigma=10))
    
   jags <- jags.model('normal.modelo',
                      data = list('x' = x, 'n' = n),
                      n.chains = 3,
                      inits = inis,
                      n.adapt = 100)
    
   #update(jags, 1000)
    
   #sam <- jags.samples(jags, c('mu', 'tau'), 1000)
   sam <- coda.samples(jags, c('mu', 'tau'), n.iter=10000, thin=10)
    
   par(mfrow=c(2,2))
   plot(sam)
   str(sam)
   summary(sam)
   HPDinterval(sam)

2. regressão linear simples

   cat( "model {
         	for (i in 1:n){
   		y[i] ~ dnorm(mu[i], tau)
   		mu[i] <- b0 + b1 * x[i]
   	}			
   	b0 ~ dnorm(0, .0001)
   	b1 ~ dnorm(0, .0001)
   	tau <- pow(sigma, -2)
   	sigma ~ dunif(0, 100)
   }", file="reglin.modelo")
    
   ## alternativa:
   ## tau ~dgamma(0.001, 0.001)
   ## sigma2 <- 1/tau
    
   n <- 20
   x <- sort(runif(n, 0, 20))
   epsilon <- rnorm(n, 0, 2.5)
   y <- 2 + 0.5*x + epsilon
    
   #write.table(data.frame(X = x, Y = y, Epsilon = epsilon),
   #            file = 'reglin.dados',
   #            row.names = FALSE,
   #            col.names = TRUE)
    
   require(rjags)
   # require(R2jags)
   # require(dclone)
    
   inis <- list(list(b0=0, b1=1, sigma=1),
   #             list(b0=1, b1=0.5, sigma=2),
                list(b0=2, b1=0.1, sigma=5))
    
   jags <- jags.model('reglin.modelo',
                      data = list('x' = x,
                                  'y' = y,
                                  'n' = n),
                      n.chains = 2,
                      # inits=inits,
                      n.adapt = 100)
    
   #update(jags, 1000)
   class(jags)
    
   sam <- jags.samples(jags,
                c('b0', 'b1', 'sigma'),
                1000)
   class(sam)
    
   sam <- coda.samples(jags,
                c('b0', 'b1', 'sigma'),
                1000)
   class(sam)
   str(sam)
   plot(sam)
    
   sam <- coda.samples(jags,
                c('b0', 'b1', 'sigma', "y"),
                1000)
   str(sam)
   int <- HPDinterval(sam)
   str(int)
   ## complementar com gráficos, resumos, inferências de interesse, etc

1. Fazer uma análise Bayesiana de uma regressão de Poisson e comparar com resultados de um MLG (modelo linear generalizado) usual. Utilize um conjunto de dados qualquer da literatura (sugestão: conjunto gala do pacote "faraway")
2. Use os dados do artigo Confiabilidade e Precisão na Estimação de Médias de Singer et al na Revista Brasileira de Estatística, v73, n. 236, jan./jun 2012 e:
   1. ajuste um modelo "ingênuo" de média e resíduo
   2. ajuste um modelo de média, efeito de local e resíduo (medidas repetidas com efeitos fixos)
   3. ajuste um modelo de efeitos aleatórias
   4. ajuste um modelo Bayesiano.

Compare e discuta os resultados.

##
## Dados simulados do modelo:
## Y_{ij} \sim N(\mu_{i}, \sigma^2_y)
##     mu_{i} = theta + b_{i}
##     b_{i} \sim N(0, \sigma^2_b)
## que, por ser normal (com ligação identidade)
## pode ser escrito por:
## Y_{ij} = \beta_0 + b_{i} + \epsilon_{ij} 
##
Ngr <-  25
Nobs <- 10
set.seed(12)
sim <- data.frame(id  = Ngr*Nobs,
                  gr  = rep(1:Ngr, each=Nobs),
                  bs  = rep(rnorm(Ngr, m=0, sd=10), each=Nobs),
                  eps = rnorm(Ngr*Nobs, m=0, sd=4)
                  )
sim <- transform(sim, y = 100 + bs + eps)
sim
 
## estimativas "naive"
resumo <- function(x) c(media=mean(x), var=var(x), sd=sd(x), CV=100*sd(x)/mean(x))
(sim.res <- aggregate(y~gr, FUN=resumo, data=sim))
var(sim.res$y[,1])
mean(sim.res$y[,2])
mean(sim$y)
 
## Modelo de efeitos aleatórios
require(lme4)
fit.lme <- lmer(y ~ 1|gr, data=sim)
summary(fit.lme)
ranef(fit.lme)
coef(fit.lme)$gr - fixef(fit.lme)
print(VarCorr(fit.lme), comp="Variance")
 
## JAGS
require(rjags)
 
sim.lst <- as.list(sim[c("gr","y")])
sim.lst$N <- nrow(sim)
sim.lst$Ngr <- length(unique(sim$gr))
mean(sim.lst$y)
 
cat("model{
    for(j in 1:N){
        y[j] ~ dnorm(mu[gr[j]], tau.e)
     }
    for(i in 1:Ngr){
        mu[i] ~ dnorm(theta, tau.b)
    }
    theta ~ dnorm(0, 1.0E-6)
    tau.b ~ dgamma(0.001, 0.001)
    sigma2.b <- 1/tau.b
    tau.e ~ dgamma(0.001, 0.001)
    sigma2.e <- 1/tau.e
    cci <- sigma2.e/(sigma2.e+sigma2.b)
}", file="sim.jags")
 
sim.jags <- jags.model(file="sim.jags", data=sim.lst, n.chains=3, n.adapt=1000)
## inits = ...
 
fit.jags <- coda.samples(sim.jags, c("theta", "sigma2.b", "sigma2.e", "cci"), 10000, thin=10)
 
summary(fit.jags)
plot(fit.jags)
 
##
require(INLA)
 
fit.inla <- inla(y ~ f(gr) , family="gaussian", data=sim)
summary(fit.inla)
sqrt(1/fit.inla$summary.hyperpar[,1])

1. (Individual ou em duplas) Encontre um conjunto de dados para cada um dos casos na avaliação semanal de 16/04. Proceda análises não Bayesianas e Bayesianas e discuta os resultados. Trazer scripts para discussão na aula de 23/04 e apresentação em 25/04.
2. Considere o modelo de verossimilhança e a priori . Mostre como obter a densidade:
   .
   Como este resultado pode ser interpretado?

require(INLA)
##
## Visualizad dados
##
data(Tokyo)
head(Tokyo)
plot(y ~ time, data=Tokyo)
## colocando na forma de proporção de dias com chuva
plot(y/2 ~ time, data=Tokyo)
##
## 1. Modelo "Nulo": só intercepto  
## estimando a probabilidade de chuva como uma constante:
fit.glm <- glm(cbind(y, n-y) ~ 1, family=binomial, data=Tokyo)
abline(h=exp(coef(fit.glm))/(1+exp(coef(fit.glm))), col=2, lty=3, lwd=3)
## como, como este modelo, todos os valores preditos são iguais bastaria fazer:
abline(h=fitted(fit.glm)[1], col=2, lty=3, lwd=3)
##
## 2. Modelo com probabilidades variando no tempo (variável/processo latente)
## modelando o logito(probabilidade) como um efeito aleatório correlacionado no tempo
## segundo um "random walk" cíclico de ordm 2
modelo = y ~ 0 + f(time, model="rw2", cyclic=T, param=c(1, 0.0001))
fit <- inla(modelo, data=Tokyo, family="binomial", Ntrials=n, control.predictor=list(compute=TRUE))
##
names(fit)
head(fit$summary.fitted.values)
## sobrepondo ao gráfico dos dados (moda, mediana e médi são praticamente indistinguíveis
with(fit, matlines(summary.fitted.values[,c(1,3:6)], lty=c(1,2,2,2,3), col=1))
##
## 3. Agora o mesmo modelo nulo anterior porém jusado pelo INLA
##
modelo0 = y ~ 1
fit0 <- inla(modelo0, data=Tokyo, family="binomial", Ntrials=n, control.predictor=list(compute=TRUE))
summary(fit0)
fit0$summary.fitted.values
with(fit0, matlines(summary.fitted.values[,c(1,3:6)], lty=c(1,2,2,2,3), col=2))
##
## 4. Modelando usando GAM (generalised additive model)
##
require(mgcv)
fit.gam <- gam(cbind(y, n-y) ~ s(time), family=binomial, data=Tokyo)
names(fit.gam)
fitted(fit.gam, se=T)
pred.gam <- predict(fit.gam, type="response", se=T, newdata=Tokyo["time"])
names(pred.gam)
with(pred.gam, matlines(cbind(fit, fit+2*se.fit, fit-2*se.fit), lty=c(1,2,2), col=4))