Practical 7

class: center, middle, title-slide

.title[
# Practical 7
]
.author[
### Olivier Gimenez
]
.date[
### last updated: 2025-11-11
]

---

## DIC is a Bayesian version of AIC

---
## DIC in Jags

``` r
storks$BUGSoutput$DIC
```

```
## [1] 207.6708
```

``` r
storks$BUGSoutput$pV
```

```
## [1] 2.959272
```

---
# Problem

---
## Model selection with DIC

* Fit models with rainfall effect, temperature effect and without any covariate to the stork data.

* Rank them with DIC.

---
# Solution

---
## Data

``` r
nbchicks <- c(151,105,73,107,113,87,77,108,118,122,112,120,122,
              89,69,71,53,41,53,31,35,14,18)

nbpairs <- c(173,164,103,113,122,112,98,121,132,136,133,137,
             145,117,90,80,67,54,58,39,42,23,23)

temp <- c(15.1,13.3,15.3,13.3,14.6,15.6,13.1,13.1,15.0,11.7,
          15.3,14.4,14.4,12.7,11.7,11.9,15.9,13.4,14.0,13.9,
          12.9,15.1,13.0)

rain <- c(67,52,88,61,32,36,72,43,92,32,86,28,57,55,66,26,
          28,96,48,90,86,78,87)
```

---
## Model with both covariates

.tiny-font[

``` r
model <- 
paste("
model
{
	for( i in 1 : N) 
		{
		nbchicks[i] ~ dbin(p[i],nbpairs[i])
		logit(p[i]) <- a + b.temp * temp[i] + b.rain * rain[i]
		}
			
# priors for regression parameters
a ~ dnorm(0,0.001)
b.temp ~ dnorm(0,0.001)
b.rain ~ dnorm(0,0.001)
			
	}
")
writeLines(model,here::here("slides", "code", "logistic.txt"))
```
]

---
## MCMC details and data

``` r
init1 <- list(a = -0.5, b.temp = -0.5, b.rain = -0.5)
init2 <- list(a = 0.5, b.temp = 0.5, b.rain = 0.5)
inits <- list(init1,init2)
parameters <- c("a","b.temp","b.rain")
nb.burnin <- 1000
nb.iterations <- 2000

datax <- list(N = 23, 
              nbchicks = nbchicks, 
              nbpairs = nbpairs, 
              temp = (temp - mean(temp)) / sd(temp), 
              rain = (rain - mean(rain)) / sd(rain))
```

---
## Run Jags

.tiny-font[

``` r
storks <- jags(data  = datax,
               inits = inits,
               parameters.to.save = parameters,
               model.file = here::here("slides", "code", "logistic.txt"),
               n.chains = 2,
               n.iter = nb.iterations,
               n.burnin = nb.burnin)
```

```
## Compiling model graph
##    Resolving undeclared variables
##    Allocating nodes
## Graph information:
##    Observed stochastic nodes: 23
##    Unobserved stochastic nodes: 3
##    Total graph size: 181
## 
## Initializing model
```

``` r
dic <- c(storks$BUGSoutput$DIC, storks$BUGSoutput$pV)
```
]

---
## Model with temperature only

.tiny-font[

``` r
# model specification
model <- 
paste("
model
{
	for( i in 1 : N) 
		{
		nbchicks[i] ~ dbin(p[i],nbpairs[i])
		logit(p[i]) <- a + b * cov[i]
		}
			
# priors for regression parameters
a ~ dnorm(0,0.001)
b ~ dnorm(0,0.001)
	}
")
writeLines(model,here::here("slides", "code", "logtemp.txt"))
```
]

---
## MCMC details and data

.small-font[

``` r
# list of lists of initial values (one for each MCMC chain)
init1 <- list(a = -0.5, b = -0.5)
init2 <- list(a = 0.5, b = 0.5)
inits <- list(init1,init2)
# specify parameters that need to be estimated
parameters <- c("a","b")
# specify nb iterations for burn-in and final inference 
nb.burnin <- 1000
nb.iterations <- 2000

# read in data
datax <- list(N = 23, 
              nbchicks = nbchicks, 
              nbpairs = nbpairs, 
              cov = (temp - mean(temp)) / sd(temp))
```
]

---

## Load R2jags to run Jags through R

.tiny-font[

``` r
storks_temp <- jags(data  = datax,
               inits = inits,
               parameters.to.save = parameters,
               model.file = here::here("slides", "code", "logtemp.txt"),
               n.chains = 2,
               n.iter = nb.iterations,
               n.burnin = nb.burnin)
```

```
## Compiling model graph
##    Resolving undeclared variables
##    Allocating nodes
## Graph information:
##    Observed stochastic nodes: 23
##    Unobserved stochastic nodes: 2
##    Total graph size: 125
## 
## Initializing model
```

``` r
dic_temp <- c(storks_temp$BUGSoutput$DIC, storks_temp$BUGSoutput$pV)
```
]

---

## Model with rainfall only

``` r
# read in data
datax <- list(N = 23, 
              nbchicks = nbchicks, 
              nbpairs = nbpairs, 
              cov = (rain - mean(rain)) / sd(rain))
```

---
## Load R2jags to run Jags through R

.tiny-font[

``` r
storks_rain <- jags(data  = datax,
                    inits = inits,
                    parameters.to.save = parameters,
                    model.file = here::here("slides", "code", "logtemp.txt"),
                    n.chains = 2,
                    n.iter = nb.iterations,
                    n.burnin = nb.burnin)
```

``` r
dic_rain <- c(storks_rain$BUGSoutput$DIC, storks_rain$BUGSoutput$pV)
```
]

---
## Model with no effect of covariates

.tiny-font[

``` r
# model specification
model <- 
paste("
model
{
	for( i in 1 : N) 
		{
		nbchicks[i] ~ dbin(p[i],nbpairs[i])
		logit(p[i]) <- a
		}
			
# priors for regression parameters
a ~ dnorm(0,0.001)
	}
")
writeLines(model,here::here("slides", "code", "lognull.txt"))
```
]

---
## MCMC details and data

.small-font[

``` r
# list of lists of initial values (one for each MCMC chain)
init1 <- list(a = -0.5)
init2 <- list(a = 0.5)
inits <- list(init1,init2)
# specify parameters that need to be estimated
parameters <- c("a")
# specify nb iterations for burn-in and final inference 
nb.burnin <- 1000
nb.iterations <- 2000

# read in data
datax <- list(N = 23, 
              nbchicks = nbchicks, 
              nbpairs = nbpairs)
```
]

---
## Load R2jags to run Jags through R

.tiny-font[

``` r
storks_null <- jags(data  = datax,
                    inits = inits,
                    parameters.to.save = parameters,
                    model.file = here::here("slides", "code", "lognull.txt"),
                    n.chains = 2,
                    n.iter = nb.iterations,
                    n.burnin = nb.burnin)
```

```
## Compiling model graph
##    Resolving undeclared variables
##    Allocating nodes
## Graph information:
##    Observed stochastic nodes: 23
##    Unobserved stochastic nodes: 1
##    Total graph size: 51
## 
## Initializing model
```

``` r
dic_null <- c(storks_null$BUGSoutput$DIC, storks_null$BUGSoutput$pV)
```
]

---
## Compare DIC

.tiny-font[

``` r
data.frame(model = c('both_covariates', 'temp', 'rain', 'none'),
           dic = c(dic[1], dic_temp[1], dic_rain[1], dic_null[1]),
           pv = c(dic[2], dic_temp[2], dic_rain[2], dic_null[2])) %>%
  arrange(dic)
```

```
##             model      dic       pv
## 1            rain 206.7205 2.696616
## 2 both_covariates 207.9183 3.209029
## 3            none 209.9280 1.006815
## 4            temp 212.4213 2.295603
```
]

* Model with rainfall only seems to be better supported by the data.

* In case models have similar DIC values, model-averaging might be useful.