Exercise 6.13

In 1986, the space shuttle Challenger exploded during takeoff, killing the seven astronauts aboard. The explosion was the result of an O-ring failure, a splitting of a ring of rubber that seals the parts of the ship together. The accident was believed to have been caused by the unusually cold weather (31°F or 0°C) at the time of launch, as there is reason to believe that the O-ring failure probabilities increase as temperature decreases.

Data on previous space shuttle launches and O-ring failures is given in the dataset challenger provided with the mcsm package. The first column corresponds to the failure indicators \(y_i\) and the second column to the corresponding temperature \(x_{i} \ \scriptsize 1 ≤ i ≤ 24\)

Solución

Primero procedamos a cargar los paquetes y otros datos necesarios.

# install.packages('mcsm')
library('MASS')
library('coda')
library('mcsm')
library('knitr')
library('VGAM')
library('plyr')
library('ggplot2')
library('rjags')
library('stats')
data(challenger)

De acuerdo con la documentación, disponible aquí, en sabemos que el conjunto de datos challenger registra la ocurencia de fallas en los O-anillos durante lanzamientos espaciales contra la temperatura al momento del lanzamiento.

De la columna oring, se tiene que:

• 1: Representa una falla
• 0: Representa un éxito

Exploremos, rápidamente los datos.

Primeras 5 observaciones del dataset challenger

oring	temp
1	53
1	57
1	58
1	63
0	66
#### Inci	so 1

Fit this dataset with a logistic regression, where:

\[ P(Y_{i} = 1 | \ x_{i}) = p(x_{i}) = \frac{exp(\alpha + \beta x_{i})}{1 + exp(\alpha + \beta x_{i})}\] using R glm function, as illustrated on page 21. Deduce the MLEs for \(\alpha\) and \(\beta\), along with standard errors.

Solución

Usando la función R glm

log_reg = glm(formula = oring~temp,family="binomial",data=challenger)
summary(log_reg)

## 
## Call:
## glm(formula = oring ~ temp, family = "binomial", data = challenger)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -1.0611  -0.7613  -0.3783   0.4524   2.2175  
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)  
## (Intercept)  15.0429     7.3786   2.039   0.0415 *
## temp         -0.2322     0.1082  -2.145   0.0320 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 28.267  on 22  degrees of freedom
## Residual deviance: 20.315  on 21  degrees of freedom
## AIC: 24.315
## 
## Number of Fisher Scoring iterations: 5

obtemenemos, los estimadores:

\[\begin{align*} \hat{\alpha} &= 15.0429 \\ \hat{\beta} &= -0.2322 \end{align*}\]

Y los errores estándar:

\[\begin{align*} se_{\alpha} &= 7.3786 \\ se_{\beta} &= 0.1082 \end{align*}\]

De acuerdo a los p-value's obtenidos, podemos afirmar que ambos estimadores son significativos en el modelo.

Inciso 2

Set up a Metropolis–Hastings algorithm with the likelihood as target using an exponential candidate for \(\alpha\) and a Laplace (double-exponential) candidate for \(\beta\). (Hint: Choose the parameters of the candidates based on the MLEs derived in 1))

Solución

Dado el modelo

\[ P(Y_{i} = 1 | \ x_{i}) = p(x_{i}) = \frac{exp(\alpha + \beta x_{i})}{1 + exp(\alpha + \beta x_{i})}\]

Se tiene que,

\[ Y_{i} \sim Bernoulli\left(p(x_{i})\right) \ \scriptsize 1 ≤ i ≤ 23\]

Así para cada \(1 ≤ i ≤ 23\),

\[ P[Y_{i} = y_{i} | x_{i}] = \left(\frac{exp(\alpha + \beta x_{i})}{1 + exp(\alpha + \beta x_{i})} \right)^{y_{i}} \left(1 - \frac{exp(\alpha + \beta x_{i})}{1 + exp(\alpha + \beta x_{i})} \right)^{1-y_{i}}\] O, equivalentemente,

\[ P[Y_{i} = y_{i} | x_{i}] = \left(\frac{exp(\alpha + \beta x_{i})}{1 + exp(\alpha + \beta x_{i})} \right)^{y_{i}} \left(\frac{1}{1 + exp(\alpha + \beta x_{i})} \right)^{1-y_{i}} \] Por lo que la función de verosimilitud es:

\[ \mathcal{L}(\alpha,\beta | \underline{y}) = \prod_{i=1}^{23} \left(\frac{exp(\alpha + \beta x_{i})}{1 + exp(\alpha + \beta x_{i})} \right)^{y_{i}} \left(\frac{1}{1 + exp(\alpha + \beta x_{i})} \right)^{1-y_{i}} \]

siguiendo la sugerencia, las distribuciones candidatas a utilizar serían:

\[\begin{align*} \alpha &\sim Exponencial(\lambda = 1/15.0429) \\ \beta &\sim Laplace(\mu = -0.2322, b = 0.0765) \end{align*}\]

Esto debido a que

\[\begin{align*} Var(\beta) &= 0.1082^{2} \\ 2 b^2 &= 0.1082^{2} \\ \sqrt{2} b & = 0.1082 \\ b &= 0.0765 \end{align*}\]

A continuación, se presenta la implementación del algoritmo:

# Cargando los datos
data(challenger)

# Fijamos la semilla
set.seed(210715)

# Parámetros globales
y = challenger[,1]
x = challenger[,2]


# Definimos la función de verosimilitud

verosimilitud = function(alpha, beta){

    # Función que devueve el valor de la función de 
    # verosimilitud para el modelo de regresión 
    # logística bajo la muestra para la muestra dada

    # Argumentos
    #   - alpha: Coeficiente alpha de la regresión (intercepto)
    #   - beta: Coeficiente beta de la regresión


    # Devuelve:
    # vero_prod: Valor de verosimilitud evaluada en alpha, beta

    vero_prod = 1

    for(i in 1:23){
        vero_prod = vero_prod * (exp(alpha + beta*x[i])/(1 + exp(alpha + beta*x[i])))^(y[i]) * 
                    (1/(1 + exp(alpha + beta*x[i])))^(1 - y[i])
    }
    return(vero_prod)
}

# Definimos la función del producto de las candidatas

producto = function(alpha, beta){

    # Función que devuelve el producto de las funciones de 
    # densidad candidatas para alpha y beta
    
    # Parámetros:
    #   - alpha: Valor de alpha tal que alpha ~ Laplace(,)
    #   - beta: Valor de beta tal que beta ~ Laplace(,)

    producto = dexp(alpha, 1/15.0429) * dlaplace(beta, −0.2322,0.0765)
    
    return(producto) 
}

# Implementación del algoritmo

MH_algorithm = function(Nsim = Nsim){

    # Función que implementa el método de Metrópolis Hasthings
    # para simular valores alpha, beta provenientes de la distribución 
    # de verosimilitud. 

    # Argumentos:
    #   - Nsim = Número de simulaciones

    # Retorna:
    #   - raplha: Valores para alpha de la simulación
    #   - rbeta: Valores de beta de la simulación

    # Inicializamos los parámetros a estimar
    ralpha = rep(rexp(1,1/15.0429),Nsim)
    rbeta = rep(rlaplace(1,−0.2322,0.0765),Nsim)


    # Simulaciones
    for(i in 2:Nsim){

        Y = c(rexp(1, 1/15.0429), rlaplace(1,−0.2322,0.0765))

        alpha_ratio = min(((verosimilitud(Y[1],Y[2]) * producto(ralpha[i-1],rbeta[i-1]))/ 
                   (verosimilitud(ralpha[i-1],rbeta[i-1]) * producto(Y[1],Y[2]))),1)

        ralpha[i] = ralpha[i-1] + (Y[1] - ralpha[i-1]) * (runif(1) < alpha_ratio)
        rbeta[i] = rbeta[i-1] + (Y[2] - rbeta[i-1]) * (runif(1) < alpha_ratio)
        
    }
    return(list(ralpha, rbeta))
}

# Usando el algortimo
Nsim = 5000
M_ejecucion = MH_algorithm(Nsim)

palpha = M_ejecucion[[1]]
pbeta = M_ejecucion[[2]]

par(mfrow=c(2, 2), cex = 0.6)
plot(palpha,type='l',
     main=expression(paste("Simulación para ", alpha)),
     ylab="",xlab="Simulaciones", col="chartreuse3")
plot(pbeta,type='l',
     main=expression(paste("Simulación para ", beta)),
     ylab="",xlab="Simulaciones", col="chocolate1")
hist(palpha,
     main=expression(paste("Histograma de ", alpha)),
     ylab="",xlab="Simulaciones", col="chartreuse3")
hist(pbeta,
     main=expression(paste("Histograma de ", beta)),
     ylab="",xlab="Simulaciones", col="chocolate1")

Inciso 3

Generate 5000 iterations of the Markov chain and construct a picture similar to Figure 6.6 to evaluate the variability of \(p(x)\) minus the observation dots.

Nsim = 5000
M_ejecucion = MH_algorithm(Nsim)

palpha = M_ejecucion[[1]]
pbeta = M_ejecucion[[2]]


# Graficando observaciones originales
plot(x,y,pch=16,cex=1.5,col="maroon4",xlab="Temperatura",ylab="Fallas",
     main="Variabilidad de p(x) menos los puntos observados")

# Curvas p(x) para los últimos 500 valores (palpha[i],pbeta[i])
for (i in 4500:Nsim){
  curve(exp(palpha[i]+pbeta[i]*x)/(1+exp(palpha[i]+pbeta[i]*x)),
        col="plum3",lwd=1.5, add=TRUE,)
}

# Curva usando los estimadores de la función glm 
curve(exp(15.0429−0.2322*x)/(1+exp(15.0429+−0.2322*x)),
      col="mediumorchid4",lwd=2.5, add=TRUE,)

Challenger data con curva \(p(x)\) (violeta oscuro) ajustada con los mínimos cuadrados de la función glm. Las curvas en violeta claro representan las muestras \(\left\lbrace \alpha^{(i)},\beta^{(i)}) \right \rbrace\) obtenidas a partir de la simulacion Monte Carlo y muestran la variabilidad de las curvas ajustadas en las últimas 500 iteraciones del total de simulaciones.

Inciso 4

Derive from this sample an estimate of the probability of failure at 60°, 50°, and 40° F along with a standard error.

Solución

Como estimadores de \(\alpha,\beta\) podemos utilizar:

\[\begin{align*} \hat{\alpha} &= \sum_{i=1}^{N_{sim}} \alpha^{(i)} \\ \hat{\beta} &= \sum_{i=1}^{N_{sim}} \beta^{(i)} \end{align*}\]

donde \(\alpha^{(i)}\) corresponde a la i-ésima simulación de \(\alpha\), mediante el algortimo de Metropólis-Hastings, análogamente para \(\beta\). Entonces, en general, para una temperatura \(t\), se tiene que la probabilidad de falla a dicha temperatura , es decir,

\[ P(Y = 1 | X = c) = \frac{exp(\hat{\alpha}+\hat{\beta} \cdot c)}{1 + exp(\hat{\alpha}+\hat{\beta} \cdot c)}\] En la siguiente tabla se muestra dicha probabilidad para los valores requeridos, usando los estimadores propuestos y los estimados mediante la función glm

Temperatura \((t)\)	\(P(Y=1\) | \(X=t)\) glm	\(P(Y=1\) | \(X=t)\) Metropolis-Hastings
40	0.9968428	0.9970414
50	0.9687171	0.9661424
60	0.7522969	0.7072793
70	0.2295065	0.1698440
80	0.0283850	0.0170288

De lo anterior, podemos conculir que para temperaturas menores a 60°F (15.5556°C) es más probable que ocurra una falla.

Implementación con RJAGS

Para la modelación con JAGS usaremos como distribuciones a priori, una gamma no informativa para \(\alpha\), respetando que la distribución de la candidata (exponencial) en los incisos anteriores tomaba valores en los reales positivos. Para \(\beta\) usaremos una distribución normal no informativa, puesto que el soporte de la distribución candidata (laplace) tiene soporte en todos los reales. En resumen utilizaremos:

\[\begin{align*} \alpha &\sim Gamma() \beta &\sim Normal() \end{align*}\]

## Hipérparámetros
y = challenger[,1]
x = challenger[,2]


# Información inicial
Nsim = 5000
alpha = 15.0429
beta = -0.2322
eta = alpha + beta * x
p = plogis(eta)
y_jags = rep(0,23)

for(i in 1:23){
  y_jags[i] = rbinom(1,size=1,prob=p[i])
}

# Carga de datos para el modelo
data <- list(y = y_jags, 
            x = x,
            n = 23)

# Especificando el nombre de los parámetros del modelo
param <- c('alpha','beta')

# Inicialización de los parámetros
inits = function(){
  list("alpha" = rexp(1,1/15.0429),
       "beta" = rlaplace(1,-0.2322,0.0765))
}


ruta = 'challenger.bug'

# Modelo
fit <- jags.model(ruta,data,inits(), n.chains = 3)

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

fit

## JAGS model:
## 
## /*
## Logit
## */
## 
## model{
## for(i in 1:n){   
##  y[i] ~ dbern(p[i])
##  logit(p[i]) <- eta[i]
##  eta[i] <- alpha + beta*x[i]
##  }
## 
## ### Prior
##  alpha ~ dnorm(0.0, 1.0E-4)
##  beta ~ dnorm(0.0, 1.0E-4)
## }
## Fully observed variables:
##  n x y

update(fit,1000)

sample <- coda.samples(fit, param, n.iter=5000, thin=1)
summary(sample)

## 
## Iterations = 2001:7000
## Thinning interval = 1 
## Number of chains = 3 
## Sample size per chain = 5000 
## 
## 1. Empirical mean and standard deviation for each variable,
##    plus standard error of the mean:
## 
##          Mean      SD  Naive SE Time-series SE
## alpha 12.5529 5.86940 0.0479235         0.9832
## beta  -0.1882 0.08466 0.0006912         0.0142
## 
## 2. Quantiles for each variable:
## 
##          2.5%     25%     50%    75%    97.5%
## alpha  2.1531  8.2596 12.3269 16.189 25.47368
## beta  -0.3733 -0.2408 -0.1849 -0.126 -0.03797

Gráficos de la simulación

plot(sample)

par(mfrow=c(2,2))
traceplot(sample)


dic.samples(fit, n.iter=2000,thin=1, type="pD")

## Mean deviance:  28.03 
## penalty 1.804 
## Penalized deviance: 29.84

# 
# dic.samples(fit, n.iter=2000,thin=1, type="popt")