EJERCICIOS-APLICATIVOS.utf8

TRABAJO DE MODELOS LINEALES GENERALIZADOS

Castañeda Sixyel$^{1a}$, Coronel Karen$^{2a}$, Diaz Astrid$^{3a}$ & Guzmán Cristian$^{4a}$

$^{1}$ Ingeniero industrial, Universidad del Atlántico, sixyelc@uninorte.edu.co

$^{2}$ Matemática, Universidad del Atlántico, kdcoronell@uninorte.edu.co $^{3}$ Matemática, Universidad del Atlántico, cdastrid@uninorte.edu.co $^{4}$ Biólogo, Universidad del Magdalena, sotojc@uninorte.edu.co $^{a}$ Estudiante de Maestría en Estadística Aplicada, Universidad del Norte

PARA TENER EN CUENTA

A continuación se presentan un conjunto de definiciones y criterios, que serán de utilidad a lo largo del desarrollo del presente trabajo

DEFINICIÓN 1 (DOBSON, 1990): Considere una sola variable aleatoria $Y$ con función de probabilidad, si es discreta, o función de densidad de probabilidad, si es continua, depende de un solo parámetro $\theta$. La distribución pertenece a la familia exponencial si se puede escribir en la forma:

\[f(y;\theta)=s(y)t(\theta)e^{a(y)b(\theta)} \] Donde $a$, $b$, $s$ y $t$ son funciones conocidas. Observe la simetría entre $y$ y $\theta$. Enfatizado en la ecuación anterior por lo que esto se puede reescribir de la forma: \[f(y;\theta)=exp[a(y)b(\theta)+c(\theta)+d(y)]\] Donde $s(y)=exp(d(y))$ y $t(\theta)=exp(c(\theta))$.

Ahora, sea $l(y,\theta)$ la función log-verosimilitud de $y$ y sea $U$ la primera derivada de $l$ con respecto a $\theta$, i.e.

\[l(\theta,y)=logf(\theta,y)=a(y)b(\theta)+c(\theta)+d(y)\] \[U=\frac{dl}{d\theta}=a(y)b'(\theta)+c'(\theta); \hspace{0.5cm} y\] \[U'=\frac{d^2l}{d\theta^2}=a(y)b''(\theta)+c''(\theta)\] Sea $E(U)=b'(\theta)E[a(Y)]+c'(\theta)$, para el caso $E(U)=0$ tenemos que: \[E[a(Y)]=\frac{-c'(\theta)}{b'(\theta)} \] Sea $Var(U)=[b'(\theta)]^2Var[a(Y)]$ y $E(-U')=-b''(\theta)E[a(Y)-c''(\theta)]$ por lo que obtenemos: \[Var[a(Y)]=[b''(\theta)c'(\theta)b'(\theta)]/[b'(\theta)]^3 \]

DEFINICIÓN 2 (MYERS, MONTGOMERY, VINING & ROBINSON, 2010): Todos los miembros de la familia exponencial de distribuciones tienen funciones de densidad de probabilidad para una respuesta observada y que se puede expresar en la forma: \[f(y;\theta,\phi)=exp\left({\frac{y(\theta)-b(\theta)}{a(\phi)}+c(y,\phi)}\right) \] Donde $i=1,...,n$, $a (·)$, $b (·)$ y $c (·)$ son funciones específicas. El parámetro $\theta$ es un parámetro de ubicación natural, y $\phi$ a menudo se denomina parámetro de dispersión. $a(\phi) = w\phi$ , $w$ son los pesos, $\Phi$ es el parámetro de dispersión, $\theta$ el parámetro de localización.

Y sea \[E(y)=a(\phi)b'(\theta)\] y \[V(y)=a(\phi)b''(\theta)\]

Los siguientes criterios serán aplicados para el análisis de los modelos:

El criterio de información akaike (AIC): Es una herramienta para la selección del modelo. Dado un conjunto de datos, varios modelos pueden ser clasificados de acuerdo a su AIC, donde el modelo que tiene el mínimo AIC es el mejor candidato. A partir del valor AIC también se puede inferir si hay modelos más o menos “empatados”. Permite comparar modelos vía stepwise de regresión logística stepwise regresión Poisson.

Desviación residual: Es una prueba de bondad de ajuste del modelo de regresión logística, los valores pequeños de deviance (o un valor P grande implican que el modelo proporciona un ajuste satisfactorio a los datos, mientras que los valores grandes de deviance implican que el modelo actual no es adecuado.

SECCIÓN 1: DISTRIBUCIONES

DISTRIBUCIÓN RAYLEIGH

Sea $\textbf{y}$ una variable aleatoria con distribución Rayleigh y parámetro $\sigma$, de modo que $\textbf{y}=(y_{1},\cdots,y_{n})^T$ con $i=1,\cdots,n$; cuya función de densidad corresponde a: \[f(\mathbf{y},\sigma)=\frac{y\cdot exp(-\frac{y^2}{2\sigma^2})}{\sigma^2}, con \hspace{0.5cm}\sigma>0.\hspace{1cm}(*)\] Inicialmente veamos que la distribución Rayleigh pertenece a la familia exponencial. Para ello tengamos en cuenta la estructura de la familia exponencial que aparece en la definición 1.

Para realizar la transformación apliquemos logaritmo a la función de densidad de la distribución Rayleigh definida en $(*)$. \[\begin{align*} log f(\mathbf{y};\sigma)&=log \left(\frac{\mathbf{y}}{\sigma^2}\cdot e^{\frac{-\mathbf{y}^2}{2\sigma^2}}\right)\\ &=log\left(\mathbf{y}\cdot e^\frac{-\mathbf{y}^2}{2\sigma^2} \right)-log(\sigma^2)\\ &=log(\mathbf{y})-\frac{\mathbf{y}^2}{2\sigma^2}-log(\sigma^2) \end{align*}\] Luego, \[f(\mathbf{y};\sigma)=exp\left[\frac{-\mathbf{y}^2}{2\sigma^2}-log(\sigma^2)+log(\mathbf{y})\right]\hspace{0.5cm}(1)\] De modo que una parametrización para que esta distribución pertenezca a la familia exponencial queda definida como sigue: \[\theta=-\frac{1}{2\sigma^2}; \hspace{0.5cm}\sigma^2=\frac{-1}{2\theta}\] \[a(\mathbf{y})=\mathbf{y}^2;\hspace{0.5cm}d(\mathbf{y})=log(\mathbf{y})\] \[b(\theta)=\frac{1}{-2\sigma^2}=\theta\] \[c(\theta)=-log(\sigma^2)=log(-2\theta)\] Por la definición 1 $E(a(\mathbf{y}))$ y $Var(a(\mathbf{y}))$ están dadas por las expresiones que se muestran a continuación, además reemplazando las igualdades anteriores que permiten que esta distribución haga parte de la familia exponencial, se llega a los siguientes resultados: \[\begin{align*} E(a(\mathbf{y}))&=-\frac{c'(\theta)}{b'(\theta)}\\ &=-\frac{\frac{1}{\theta}}{1}=-\frac{1}{\theta}\\ E(a(\mathbf{y}))&=2\sigma^2 \end{align*}\]\[ \]\[\begin{align*} Var(a(\mathbf{y}))&=\frac{b''(\theta)c'(\theta)-c''(\theta)b'(\theta)}{b'(\theta)^3}\\ &=\frac{0(1/\theta)-(-1/\theta^2)(1)}{(1)^3}=\frac{1}{\theta^2}\\ &=\frac{1}{\theta}\cdot\frac{1}{\theta}=(-2\sigma^2)\cdot(-2\sigma^2)\\ Var(a(\mathbf{y}))&=4\sigma^4 \end{align*}\] Sea $y_{i}$ una sucesión de variables aleatorias con distribución Rayleigh tal que $y_i\geq 0$; $\forall$ $i=1,...,n$ se construye la función de máxima verosimilitud para $\textbf{y},$ como sigue: \[\begin{align*} L(\sigma,\mathbf{y}) &= \prod_{i=1}^n f(y_i,\sigma)\\ L(\sigma,\mathbf{y})&=\prod_{i=1}^n\frac{y_{i}}{\sigma^2}\cdot e^{\frac{-y_{i}^2}{2\sigma^2}}\\ &=\prod_{i=1}^{n}\frac{y_{i}}{\sigma^2}e^{\frac{-\sum_{i=1}^{n}y_{i}^2}{2\sigma^2}}\hspace{0.5cm}(2)\\ \end{align*}\] Aplicando logaritmo a la expresión (2) tenemos que la función log-verosimilitud para $\mathbf{y}$ es: \[\begin{align*} l(\sigma,\mathbf{y}) &= log\left(\frac{1}{\sigma^{2n}}\prod_{i=1}^{n}y_{i}e^{\frac{-\sum_{i=1}^{n}y_{i}^2}{2\sigma^2}}\right)\\ &=log\left(\frac{1}{\sigma^{2n}}\prod_{i=1}^{n}y_{i}\right)-\sum_{i=1}^{n}\frac{y_{i}^{2}}{2\sigma^2}\\ &=log\left(\prod_{i=1}^{n}y_{i}\right)-log (\sigma^{2n})-\sum_{i=1}^{n}\frac{y_{i}^{2}}{2\sigma^2}\\ &=\sum_{i=1}^{n}log(y_{i})-2nlog (\sigma)-\sum_{i=1}^{n}\frac{y_{i}^{2}}{2\sigma^2}\hspace{0.5cm}(3)\end{align*}\] Por definición tenemos que el vector score corresponde a $U=\frac{\partial l(\textbf{y},\theta)}{\partial \theta}$, así que procedemos a hallar la derivada parcial de la función log-verosimilitud con respecto al vector de parámetros, en este caso la derivada parcial de (3) con respecto a $\sigma$ es: \[U=\frac{\partial l(\sigma,\mathbf{y})}{\partial\sigma} =-\frac{2n}{\sigma}+\sum_{i=1}^{n}\frac{y_{i}^{2}}{\sigma^3}\hspace{0.5cm}\]

Para realizar la estimación de los parámetros por máxima verosimilitud igualamos a cero la expresión $\frac{\partial l(\textbf{y},\theta)}{\partial \sigma}$ donde realizando las operaciones aritméticas debidas, obtenemos: \[\begin{align*} 0&=\frac{\partial l(\textbf{y},\sigma)}{\partial \sigma}\\ 0&=-\frac{2n}{\sigma}+\sum_{i=1}^{n}\frac{y_{i}^{2}}{\sigma^3}\\ 0&=-2n\sigma^2+\sum_{i=1}^{n}y_{i}^{2}\\ 2n\sigma^2&=\sum_{i=1}^{n}y_{i}^{2}\\ \hat \sigma&=\sqrt{\frac{y_{i}^{2}}{2n}}\hspace{0.5cm}(4) \end{align*}\] Tenemos que (4) es el estimador por máxima verosimilitud del parámetro $\sigma$, para una distribución Rayleigh.

Para calcular la matriz de información de Fisher, procedemos a calcular las derivadas parciales de segundo orden de la función log-verosimilitud con respecto al vector de parámetros, en este caso la segunda derivada de (3) con respecto a $\sigma$ es: \[I=-\frac{\partial^{2}l(\textbf{y},\sigma)}{\partial \sigma^2}=-\frac{2n}{\sigma^2}+3\sum_{i=1}^{n}\frac{y_{i}^{2}}{\sigma^4}\hspace{0.5cm}(5)\] Luego $(5)$ corresponde a la matriz de información de la variable aleatoria $\textbf{y}$ con distribucón Rayleigh.

DISTRIBUCIÓN BINOMIAL NEGATIVA

Consideremos una variable aleatoria $y_{i}$ con distribución binomal negativa y parámetros $r$ y $\pi$; es decir $y_{i}\sim BN(r,\pi)$, donde $\mathbf{y}=(y_1,...,y_n)^T$ con $i=1,\cdots,n$.Se define la función de probabilidad asociada a $\textbf{y}$ como sigue: \[ (1)\hspace{0.5cm} f(\textbf{y},\pi)=\begin{pmatrix} \textbf{y}+r-1\\r-1\end{pmatrix}\pi^r(1-\pi)^{\textbf{y}};\hspace{0,5cm} r\in \mathbb{N}; \hspace{0.5cm}0<\pi<1; \hspace{0.2cm}\textbf{y}>\mathbf{0}\] Donde $\textbf{y}$ sería el número de fracasos necesarios antes de conseguir el $r$-ésimo éxito, con probabilidad de éxito $\pi$.

Inicialmente verifiquemos que esta distribución pertenece a la familial exponencial. Para ello utilizaremos la parametrización descrita en la definición 1:

Para realizar la transformación apliquemos logaritmo a la función de probabilidad de la distribución binomial negativa definida en $(1)$. \[\begin{align*} log(f(\mathbf{y},\pi))&=log\left \{ \begin{pmatrix} \mathbf{y}+r-1\\r-1\end{pmatrix}\pi^r(1-\pi)^{\mathbf{y}}\right \}\\ &= log\begin{pmatrix} \mathbf{y}+r-1\\ r-1\end{pmatrix}+log\pi^r(1-\pi)^{\mathbf{y}}\\ &=log\begin{pmatrix} \mathbf{y}+r-1\\ r-1\end{pmatrix} + r log\pi+\mathbf{y} log(1-\pi)\\ &= \mathbf{y} log(1-\pi)+(-rlog\pi)+log\begin{pmatrix} \mathbf{y}+r-1\\ r-1\end{pmatrix} \end{align*}\] Luego \[f(\mathbf{y},\pi)=exp\left \{\mathbf{y} log(1-\pi)+(-rlog\pi)+log\begin{pmatrix} \mathbf{y}+r-1\\ r-1\end{pmatrix}\right \}\] De modo que una parametrización, que hace que esta distribución pertenezca a la familia exponencial queda definida como sigue: \[\theta=log(1-\pi);\hspace{0.5cm}\pi=1-e^{\theta}\] \[a(\phi)=1;\hspace{0.5cm} c(\mathbf{y,\phi})=log\begin{pmatrix} y+r-1\\r-1 \end{pmatrix}\] \[b(\theta)=-r log\pi=-r log(1-e^{\theta})\] Para esta parametrización en particular $E(y_{i})$ y $Var(y_{i})$ están dadas por $E(y_{i})=a(\phi)b'(\theta)$ y la $Var(y_{i})=a(\phi)b''(\theta)$.

Se puede comprobar facilmente que: \[{b}'(\theta)=-r\frac{-e^{\theta}}{1-e^{\theta}}=\frac{re^{\theta}}{1-e^{\theta}} \hspace{5cm}(2)\] \[{b}''(\theta)=\frac{re^{\theta}(1-e^{\theta})+e^{\theta}(re^{\theta})}{(1-e^{\theta})^2}\hspace{4.5cm}(3)\] Luego reemplazando $(2)$ y $(3)$ en $E(y_{i})$ y $Var(y_{i})$ tenemos: \[E(y_{i})=1\frac{re^{\theta}}{1-e^{\theta}}=\frac{r(1-\pi)}{\pi}\] \[Var(y_{i})=1\frac{re^{\theta}}{(1-e^{\theta})^2}=\frac{r(1-\pi)}{\pi^2}\] Dado $y_{i}$ con distribución binomal negativa con f.d.p definida como $f(y_{i},\pi)=\begin{pmatrix} y_{i}+r-1\\r-1\end{pmatrix}\pi^r(1-\pi)^{y_{i}};\hspace{0,5cm} r\in \mathbb{N}; \hspace{0.5cm}0<\pi<1\hspace{0.5cm}y_{i}=1,2,3,\cdots,$ se construye su función por máxima verosimilitud como sigue: \[\begin{align*} L(\pi,\mathbf{y})&=\prod_{i=1}^nf(y_i,\pi)\\ &=\prod_{i=1}^n\begin{pmatrix}y_i+r-1\\ r-1\end{pmatrix}\pi^r(1-\pi)^{y_i}\\ &=\begin{bmatrix}\prod_{i=1}^n\begin{pmatrix}y_i+r-1\\ r-1\end{pmatrix}\end{bmatrix} \pi^{nr} (1-\pi)^{\sum_{i=1}^ny_i}\hspace{0.5cm}(4) \end{align*}\] Aplicando logaritmo a la expresión $(4)$ que corresponde a la función de verisimilitud de $\mathbf{y}$; obtenemos la función log-verosimilitud, note a continuación la construcción: \[\begin{align*} l(\pi,\mathbf{y} )&=log\left ( \begin{bmatrix}\prod_{i=1}^n\begin{pmatrix}y_i+r-1\\ r-1\end{pmatrix}\end{bmatrix} \pi^{rn}(1-\pi)^{\sum_{i=1}^n y_i}\right )\\ &=log\begin{bmatrix}\prod_{i=1}^n\begin{pmatrix}y_i+r-1\\ r-1\end{pmatrix}\end{bmatrix} +nrlog\pi+\sum_{i=1}^ny_ilog(1-\pi) \\ &=\sum_{i=1}^n log\begin{pmatrix}y_i+r-1\\ r-1\end{pmatrix}+nrlog\pi+\sum_{i=1}^ny_ilog(1-\pi) \end{align*}\] Para realizar la estimación del parámetro $\pi$ por máxima verosimilitud igualamos a cero la expresión $\frac{d}{d\pi} l(\pi,\mathbf{y})$ y realizamos las operaciones aritméticas concernientes, tenemos que si: \[\frac{d l(\pi,\mathbf{y})}{d\pi}=\frac{nr}{\pi}-\frac{\sum_{i=1}^ny_i}{(1-\pi)}.\hspace{0.5cm}(5)\] Haciendo $\frac{d l(\pi,\mathbf{y})}{d\pi}=0,$ se obtiene entonces que: \[0=\frac{nr}{\pi}-\frac{\sum_{i=1}^n y_{i}}{(1-\pi)}\] \[\frac{\sum_{i=1}^ny_i}{(1-\pi)}=\frac{nr}{\pi}\] \[\pi\sum_{i=1}^ny_i=nr-nr\pi\] \[\pi\left(\sum_{i=1}^ny_i+nr\right )=nr\] \[\hat{\pi}=\frac{nr} {{\sum_{i=1}^{n} y_{i}}+nr}\hspace{0.5cm}(6)\] Luego $(6)$ corresponde al estimador por máxima verosimilitud del parámetro $\pi.$ Además tenemos que el vector score corresponde a $U=\frac{\partial l(\textbf{y},\theta)}{\partial \theta}$, procediendo a hallar la derivada parcial de la función log-verosimilitud con respecto a $\pi$ llegamos al resultado en $(5).$

Para calcular la matriz de información de Fisher $I$, partimos de la definición que $Var(U)=I;$ donde $U$ es el vector score que corresponde a $(5)$, para realizar esta prueba usaremos las propiedades de la varianza para proceder aritméticamente y usaremos un resultado que ya está demostrado y es que $Var(y_i)=\frac{r(1-\pi)}{\pi^2}$, sigamos: \[\begin{align*} Var(U)&=Var\left ( \frac{nr}{\pi}-\frac{\sum_{i=1}^n}{(1-\pi)} \right )\\ &=Var\left ( \frac{nr}{\pi}\right )+Var\left (\frac{\sum_{i=1}^ny_i}{(1-\pi)} \right ) \\ &=0+\frac{1}{(1-\pi)^2}\sum_{i=1}^nVar(y_i)\\ &=\frac{1}{(1-\pi)^2}\frac{nr(1-\pi)}{\pi^2}\\ \end{align*}\] Así pues, llegamos a que $I=\frac{nr}{(1-\pi)\pi^2}.$

DISTRIBUCIÓN NORMAL INVERSA

Sea $\mathbf{y_i} \sim N(\mu,\lambda)$, y sea $\mathbf{y}=(y_1,...,y_n)^T$ una suscesión de variables para $i=1,...,n$, se define la función de densidad como: \[f(Y)= \left(\frac{\lambda}{2\pi y^3}\right)^\frac{1}{2}exp\left(\frac{-\lambda(y-\mu)^2}{2\mu^2 y}\right)\] Verifiquemos que la distribución normal inversa pertenece a la familia exponencial. Para ello se tendrá en cuenta la estrutura descrita en la definición 2.

\[\begin{align*} f(Y)&=exp\left(log\left( \left(\frac{\lambda}{2\pi y^3}\right)^\frac{1}{2}exp\left(\frac{-\lambda(y-\mu)^2}{2\mu^2 y}\right)\right)\right)\\ &=exp\left( \frac{-1}{2}log\lambda-\frac{1}{2}log(2\pi y^3)-\frac{\lambda y}{2\mu^2}+\frac{\lambda}{\mu}-\frac{\lambda}{2y} \right)\\ &=exp\left(-\frac{\lambda y}{2\mu^2}+\frac{\lambda}{\mu}-\frac{\lambda}{2y}- \frac{1}{2}log\lambda-\frac{1}{2}log(2\pi y^3) \right)\hspace{0.5cm}(1) \end{align*}\]

Tomando: $\frac{\theta y}{a(\phi)}=\frac{-\lambda}{2\mu^2}$

Entonces, tenemos que: \[\theta=\frac{1}{\mu^2};\hspace{0.5cm}\frac{1}{a(\phi)}=\frac{-\lambda}{2} \] Despejando $\theta$,$a(\Phi)$ obtenemos: \[\theta^{1/2} =\frac{1}{\mu};\hspace{0.5cm}{a(\phi)}=\frac{-2}{\lambda};\hspace{0.5cm}\phi=\frac{1}{\lambda}\] Entonces $a(\phi)=-2\phi$

Reemplazando estas equivalencias en la expresión (1) y reescribiendo la expresión llegamos a que:

\[\begin{align*} f(Y)&=exp\left[\frac{\theta y}{-2\phi}+\frac{\theta^{1/2}}{\phi}-\frac{1}{\phi 2y}-\frac{1}{2}log\phi-\frac{1}{2}log2\pi y^3 \right]\\ &=exp\left(\frac{\theta y-2\theta^{1/2}}{-2\phi}-\frac{1}{2\phi y}-\frac{1}{2}log(2\pi\phi y^3) \right) \end{align*} \] De modo que una parametrización para que esta distribución pertenezca a la familia exponencial queda definida como sigue: \[b(\theta)=2\theta^{1/2}; \hspace{0.5cm}\theta^{1/2}=\frac{1}{\mu} \] Por la definición 1 la $E(Y)$ y $Var(Y)$ están dadas por:

\[\begin{align*} E(Y)&=b'(\theta)\\ &=\theta^{-1/2}\\ &=\frac{1}{\theta^{1/2}}\\ &=\mu \end{align*}\]

\[\begin{align*} Var(Y)&=a(\phi)b''(\theta)\\ &=\frac{-2}{\lambda}.\frac{-1}{2}.\theta^{-3/2}\\ &=\frac{1}{\lambda}.\theta^{-3/2}\\ &=\frac{\mu^3}{\lambda} \end{align*}\] La función de máxima verosimilitud está dada por: \[\begin{align*} L(\theta,Y)&=\prod_{i=1}^n f(y_i,\theta)\\ &=\prod_{i=1}^n\left[\left(\frac{\lambda}{2\pi y^3}\right)^\frac{1}{2}exp\left(\frac{.\lambda(y-\mu)^2}{2\mu^2 y}\right) \right]\hspace{0.5cm}(2) \end{align*}\] Aplicando log a la expresión (2) obtenemos el log-verosimilitud

\[\begin{align*} l(\theta,y)&=\frac{n}{2}log\lambda-\frac{n}{2}log2\pi-\frac{3}{2}\sum_{i=1}^{n}logy_i-\frac{\lambda}{2\mu^2}\sum_{i=1}^n\frac{(y_i-\mu)^2}{y_i} \\ &=\frac{-\lambda}{2\mu^2}\sum_{i=1}^n \frac{(y_i-\mu)^2}{y_i}+\frac{n}{2}log\lambda-\frac{n}{2}log2\pi-\frac{3}{2}\sum_{i=1}^nnlogy_i \end{align*}\]

Por definición, el vector score se define como $U=∂l(y,θ)/∂θ$, entonces procedemos a hallar las derivadas parciales de la función log-verosimilitud con respecto a la de cada uno de los parámetros $\mu$ y $\lambda$ ; y las denotaremos por $U1$ y $U2$, respectivamente:

\[\begin{align*} U1&=\frac{∂(\mu,\lambda)}{∂\mu}\\ &=\frac{n\lambda}{\mu^3}(\bar{y}-\mu)\\ \end{align*}\]

Igualando la expresión anterior a cero tenemos que:

\[\begin{align*} \frac{n\lambda}{\mu^3}(\bar{y}-\mu)&=0\\ (\bar{y}-\mu)&=0 \end{align*}\]

\[\therefore\hat{\mu}=\bar{y}\] Luego, \[\begin{align*} U2&=\frac{∂(\mu,\lambda)}{∂\lambda}\\ &=\frac{n}{2\lambda}-\frac{n}{2\mu}\sum_{i=1}^{n}\frac{(y_i-\mu)^2}{y_i} \\ \end{align*}\]

Igualando a cero tenemos que: \[\begin{align*} \frac{n}{2\lambda}-\frac{1}{2\mu^2}\sum_{i=1}^{n}\frac{(y_i-\mu)^2}{y_i}&=0 \\ \frac{n}{2\lambda}-\frac{1}{2\bar{y}^2}\sum_{i=1}^{n}\frac{(y_i-\bar{y})^2}{y_i}&=0 \\ \frac{n}{2\lambda}-\frac{1}{2\bar{y}^2}\sum_{i=1}^{n}\left(y_i-2\bar{y}+\frac{y^2}{y_i} \right)&=0 \\ \frac{n}{2\lambda}-\frac{1}{2\bar{y}^2}\left(n\bar{y}-2n\bar{y}-\bar{y}^2\sum_{i=1}^n\frac{1}{y_i} \right)&=0\\ \end{align*}\]

Haciendo los despejes correspondientes llegamos a que:

\[\therefore\hat{\lambda}=\frac{1}{\sum_{i=1}^{n}y_i+\frac{n}{\bar{y}}}\] La matriz de información son las segundas derivadas parciales con respecto a $\mu$ y $\lambda$, entonces:

\[I=\begin{equation} \begin{pmatrix} -\frac{n\hat\lambda}{\bar{y}^6} & 0\\ 0 & \frac{-n\hat\lambda^{-2}}{2} \\ \end{pmatrix} \end{equation}\]

DISTRIBUCIÓN GAMMA

Sea $\textbf{y}\sim\Gamma(\alpha;\beta)$ con $\textbf{y}=(y_1,...,y_n)^T$ una sucesión de v.a. para $i=1,\cdots,n$ se define la función de densidad como: \[f(y;\alpha;\beta)=\frac{\beta^{\alpha}}{\Gamma(\alpha)}y^{\alpha-1}exp^{-\beta y}; con\hspace{0.5cm}\alpha,\beta,y>0\] Donde $\Gamma(\alpha)=\int_0^{\infty}y^{\alpha-1}exp^{-\alpha}dx$

Primero veamos que la distribución Gamma pertenece a la familia exponencial. Para ello tengamos en cuenta la estructura de la familia exponencial que aparece en la definición 2.

Para realizar la transformación apliquemos log a la función de densidad de la distribución gamma. \[\begin{align*} log f(y;\alpha;\beta)&=log(\frac{\beta^{\alpha}}{\Gamma(\alpha)}y^{\alpha-1}exp^{-\beta y})\\ &=\alpha log\beta-log[\Gamma(\alpha)]+(\alpha-1)log(y)-\beta y\\ &=-\beta y+ \alpha log(\beta)+(\alpha-1)log(y)-log[\Gamma(\alpha)] \end{align*}\] Luego \[f(y;\alpha;\beta)=exp\left \{ -\beta y + \alpha log(\beta)+(\alpha-1)log(y)-log[\Gamma(\alpha)]\right \}\hspace{0.5cm}(1)\] Reescribiendo la expresión obtenida en (1) tenemos que: \[\begin{align*} f(y;\alpha;\beta)&=exp\left \{ \frac{-\frac{\beta y}{\alpha}+\frac{\alpha log \beta}{\alpha}}{{\frac{1}{\alpha}}}+(\alpha-1)log(y)-log[\Gamma(\alpha)]\right \}\\ &=exp\left \{ \frac{{\frac{\beta y}{\alpha}-log(\beta)}}{{-\frac{1}{\alpha}}}+(\alpha-1)log(y)-log[\Gamma(\alpha)] \right \}\hspace{0.2cm}(**) \end{align*}\] Tomando que : \[\theta=\frac{\beta}{\alpha};\hspace{0.5cm} \phi=\frac{1}{\alpha}\] \[a(\phi)=-\frac{1}{\alpha}(*);\hspace{1cm}\alpha=\frac{1}{\phi}\] Despejando $\beta$ en (*) y reemplazando $\alpha$ llegamos a la siguiente expresión: \[\begin{align*} \beta&=\theta \alpha\\ \beta&=\theta\frac{1}{\phi}\\ \beta&=log(\theta)-log(\phi) \end{align*}\] Reemplazando estas equivalencias en la ecuación (**) y reescribiendo la expresión llegamos a que: \[\begin{align*} f(y;\alpha;\beta)&=exp\left \{ \frac{\theta y - log(\theta)}{-\phi}+\frac{log(\phi)}{\phi}+(\frac{1}{\phi}-1)log(y)-log[\Gamma(\frac{1}{\phi})]\right \}\\ &=exp\left \{ \frac{\theta y-log(\theta)}{-\phi} +c(y,\phi)\right \} \end{align*}\] De modo que una parametrización para que esta distribución pertenezca a la familia exponencial queda definida como sigue: \[b(\theta)=log(\theta); \hspace{0.5cm} a(\phi)=-\phi\] \[\theta=\frac{\beta}{\alpha};\hspace{0.5cm}\phi=\frac{1}{\alpha}\] \[c(y,\phi)=\frac{log(\phi)}{\phi}+(\frac{1}{\phi}-1)log(y)-log[\Gamma(\frac{1}{\phi})]\] Por la definición 2 $E(Y)$ y $Var(Y)$ están dadas por: \[\begin{align*} E(Y)&=b´(\theta)\\ &=\frac{1}{\theta}=\frac{1}{{\frac{\beta}{\alpha}}}\\ &=\frac{\alpha}{\beta} \end{align*}\]\[ \]\[\begin{align*} Var(Y)&=a(\phi)b´´(\theta)\\ &=-\phi(-\frac{1}{\theta^2}\\ &=\frac{\alpha}{\beta}=\frac{\phi}{\theta^2}=\frac{\frac{1}{\alpha}}{(\frac{\beta}{\alpha})^2}=\frac{\frac{1}{\alpha}}{\frac{\beta^{2}}{\alpha^{2}}}\\ &==\frac{\alpha^2}{\beta^{2}\alpha}=\frac{\alpha}{\beta^2} \end{align*}\] Función de máxima verosimilitud Sea $y_{i}\sim\Gamma(\alpha;\beta)$ con $y_i\geq 0$; $\forall$ $i=1,...,n$. su función de máxima verosimilitud se define como: \[\begin{align*} L(\alpha,\beta,Y) &= \prod_{i=1}^n f(y_i,\alpha,\beta)\\ L(\alpha,\beta,Y) &=\prod_{i=1}^n\frac{\beta^\alpha}{\Gamma(\alpha)}y_i^{\alpha-1exp^{-\beta y_i}}\\ L(\alpha,\beta,Y)&=\frac{\beta^{n\alpha}}{\Gamma(\alpha)^{n}}exp^{-\beta\sum_{i=1}^{n}y_i}\left \{\prod_{i=1}^ny_i \right \}^{\alpha-1}\hspace{0.3cm}(4)\end{align*}\] Aplicando log a la expresión (4) tenemos que el log-verosimilitud es: \[\begin{align*} l(\alpha,\beta,Y)&=log(\frac{\beta^{n\alpha}}{\Gamma(\alpha)^{n}}exp^{-\beta\sum_{i=1}^{n}y_i}\left \{\prod_{i=1}^ny_i \right \}^{\alpha-1}\\ l(\alpha,\beta,Y) &=n \alpha log(\beta)-nlog \Gamma(\alpha)-\beta\sum_{i=1}^ny_i+(\alpha-1)\sum_{i=1}^nlog(y_i)\end{align*}\] Por definición el vector score se define como $U=\frac{\partial l(\textbf{y},\theta)}{\partial \theta}$, entonces procedemos a hallar las derivadas parciales de la función log-verosimilitud con respecto a la de cada uno de los parámetros $\beta$ y $\alpha$ ; y las denotaremos por $U1$ y $U2$, respectivamente: \[U_{1}=\frac{\partial}{\partial\beta} l(\alpha,\beta,Y)=\frac{n\alpha}{\beta}-\sum_{i=1}^ny_i\hspace{0.5cm}(5)\] Para realizar la estimación de los parámetros por máxima verosimilitud igualamos a cero $\frac{\partial}{\partial\beta} l(\alpha,\beta,Y)$ donde realizando las operaciones aritméticas concernientes obtenemos: \[\begin{align*} 0&=\frac{n\widehat{\alpha}}{\beta}-\sum_{i=1}^ny_i\\ \frac{n\widehat{\alpha}}{\widehat{\beta}}&=\sum_{i=1}^ny_i\\ \frac{\widehat{\alpha}}{\widehat{\beta}}&=\frac{\sum_{i=1}^n y_{i}}{n}\\ \overline{Y}&=\frac{\widehat{\alpha}}{\widehat{\beta}} \end{align*}\] $\therefore \hat \beta=\frac{\widehat{\alpha}}{\overline{Y}}$

Por otro lado, \[U_{2}=\frac{\partial}{\partial\alpha}l(\alpha,\beta,Y)=nlog(\beta)-n\frac{1}{\Gamma(\alpha)}{\Gamma}'(\alpha)+\sum_{i=1}^nlog(y_i)\hspace{0.5cm}(6)\] Donde:

$\frac{{\Gamma}'}{\Gamma(\alpha)}=\Psi(\alpha)\rightarrow$ digamma función

entonces tenemos que:

$\therefore \Psi(\alpha)=log\beta+\frac{1}{n}\sum_{i=1}^{n-1}log(Y_n)$ así el vector score queda definido como $U=(U_{1},U_{2})$; donde $U_{1}$ y $U_{2}$ corresponden a (5) y (6), respectivamente.

Para calcular la matriz de información de Fisher, procedemos a calcular las derivadas parciales de segundo orden: \[I(\alpha,\beta)_{11}=\frac{\partial^{2}l(\alpha,\beta,Y)}{\partial \alpha^2}=\frac{n\partial^2}{\partial\alpha^2}log\Gamma(\alpha)\] \[I(\alpha,\beta)_{22}=\frac{\partial^{2}l(\alpha,\beta,Y)}{\partial \beta^2}=\frac{n\alpha}{\beta^2}\] Debido a que la matriz de Fisher es simétrica, se cumple que: \[I(\alpha,\beta){12}=I(\alpha,\beta){21}=\frac{\partial^{2}l(\alpha,\beta,Y)}{\partial\alpha \partial \beta}=-n\frac{1}{\beta}\] Luego, la matriz de información de Fisher queda definida de la siguiente manera: \[\begin{align*} I(\alpha,\beta)&=\begin{bmatrix} \frac{n\partial^2}{\partial\alpha^2}log\Gamma(\alpha)&-\frac{n}{\beta} \\-\frac{n}{\beta} &\frac{n\alpha}{\beta^2}\end{bmatrix}\\ I(\alpha,\beta)&=n\begin{bmatrix} \frac{\partial^2}{\partial\alpha^2}log\Gamma(\alpha)&-\frac{1}{\beta} \\-\frac{1}{\beta} &\frac{\alpha}{\beta^2}\end{bmatrix} \end{align*}\] DISTRIBUCIÓN CHI-CUADRADO

Sea $\mathbf{y}= [y_1,...,y_n]^t$ una sucesión de variables aleatorias independientes $\mathbf{y}\sim X_{n}^{2},\hspace{0.5cm} y>0,\hspace{0.5cm} n>0,\hspace{0.5cm} i=1,...,n$

Veamos si la distribución chi-cuadrado se puede expresar como familia exponencial.

Para ello tomamos la estructura de la familia exponencial dada en la definición 1. La función de densidad está dada por:

\[f(y)=\frac{1}{2^{n/2}\Gamma(n/2)}y^{n/2-1}e^{-y/2};\hspace{0.5 cm} donde \hspace{0.5cm} \Gamma(n/2)=\int_{0}^{\infty}y^{n/2-1}e^{-y}\ dy\] \[\begin{align*} logf(y)&=log\left(\frac{1}{2^{n/2}\Gamma(n/2)}y^{n/2-1}e^{-y/2}\right)\\ &=log\ 1-\frac{n}{2}log\ 2+log(\Gamma(n/2))+\left(\frac{n}{2}-1\right)log(y)-y/2\\ &=-\frac{y}{2}-\frac{n}{2}log(2)+\left(\frac{n}{2}-1\right)log(y)-log\ \Gamma(n/2)\\ f(y;n)&=exp\left(-\frac{y}{2}+\left(\frac{n}{2}-1\right)log(y)-\frac{n}{2}log(2)-log\ \Gamma(n/2)\right) \hspace{0.5cm}(1) \end{align*}\] tomando:

\[v=\frac{n}{2}-1\ \rightarrow v+1=\frac{n}{2} \ \rightarrow 2(v+1)=n\] Reemplazando en (1) las expresiones obtenidas en el paso anterior, tenemos que: \[f(y;n)=exp\left(-\frac{y}{2}+v\ log(y)-(v+1)log(2)-log\ \Gamma(v+1)\right)\] De modo que una parametrización para que esta distribución pertenezca a la familia exponencial queda definida como sigue: \[a(y)=log(y),\hspace{0.5cm} d(y)=-y/2\] \[b(\theta)=v, \hspace{0.5cm} c(\theta)=-[(v+1)log(2)-log\ \Gamma(v+1)]\] Por la definición 1 $E(Y)$ y $Var(Y)$ están dadas por: \[E[a(y)]=\frac{-c'(\theta)}{b'(\theta)} \] \[c'(\theta)=-(\Psi(v+1)+log2)=-\left(\Psi\left(\frac{n}{2}\right)+log\ 2\right)\] Donde, $\Psi\left(\frac{n}{2}\right)$ es la digamma de $\left(\frac{n}{2}\right)$ en $R$. \[b'(\theta)=v=1\] Entonces, \[E[a(y)]=\Psi\left(\frac{n}{2}\right)+log(2)\] sea $b''(\theta)=0$ y $c''(\theta)=\Psi(n-1)=\Psi'(\frac{n}{2})$, donde $\Psi'(\frac{n}{2})$ es la trigamma $\frac{n}{2}$ en $R$ \[Var[a(y)]=-\Psi'\left(\frac{n}{2}\right)\] Sea $y_i\sim X_{n}^{2}$ con: $yi≥0$; $∀ i=1,...,n.$, se construye su función de máxima verosimilitud \[\begin{align*} L(n,Y) &= \prod_{i=1}^n f(y_i,n)\\ L(n,y) &=\prod_{i=1}^n\left(\frac{1}{2^{n/2}\Gamma(n/2)}y^{n/2-1}e^{-y/2}\right)\\ L(n,Y)&=\frac{1}{2^{n/2}\Gamma(n/2)}\left \{\prod_{i=1}^ny_i \right \}^{{n/2}-1}exp^{-1/2\sum_{i=1}^{n}y_i}\hspace{0.3cm}(2)\end{align*}\] Aplicando log a la expresión (2) tenemos que el log-verosimilitud es: \[\begin{align*} l(\alpha,\beta,Y) &=-[n/2log(2)+nlog \Gamma(n/2)]+(n/2-1)\sum_{i=1}^nlog(y_i)-1/2\sum_{i=1}^n y_i\end{align*}\]

De la definición, el vector score se define como $U=∂l(y,n)/∂n$, entonces procedemos a hallar las derivadas parciales de la función log-verosimilitud con respecto al parámetro n ; y lo denotaremos por: \[\begin{align*} U=\frac{d(n,y)}{d(n)}&=-1/2log2+\frac{d(n,y)}{d(n)}log\Gamma(n/2)+1/2\sum_{i=1}^{n}logy_{i}\\ &=-1/2log2+\frac{1}{\Gamma(n/2)}\Gamma'(n/2)+1/2\sum_{i=1}^{n}logy_{i}=0 \end{align*}\] \[\therefore \Psi(n/2)=\frac{1}{2}log2-\frac{1}{2}l\sum_{i=1}^{n}log y_{i}\] Así el vector score definido como U es igual a la expresión anterior.

Para calcular la matriz de información de Fisher, procedemos a calcular las derivadas parciales de segundo orden, entonces tenemos que: \[I=\frac{d^2(n,y)}{dn}= \Psi'(n/2)\] 3. Solucionar el ejercicio de la normal considerando $\mu$ y $\sigma^2$ desconocidos

Sea $\mathbf{y}\sim N (\mu,\sigma^2)$ $\mu$ y $\sigma^2$ desconocidos La función de densidad de probabilidad para la distribución normal es:

\[f(y)=(2\pi\sigma^{2})^{-1/2} exp(-\frac{1}{2}\frac{{(y-\mu)^2}}{\sigma^2});\hspace{0.5cm} y,\sigma,\mu>0\] Primero probemos que hace parte de la familia exponencial. Para ello tomamos la estructura de la familia exponencial dada en la definición 2. La función de densidad está dada por:

\[\begin{align*} logf(Y)&=log((2\pi\sigma^{2})^{-1/2} exp(-\frac{1}{2}\frac{{(y_i-\mu)^2}}{\sigma^2}))\\ &=-\frac{1}{2}log(2\pi\sigma^2)-\frac{1}{2}\frac{(y-\mu)^2}{\sigma^2}\\ \end{align*}\]

\[f(Y)=exp[\frac{y\mu-(1/2)\mu^2}{\sigma^2}-\frac{1}{2}(\frac{y^2}{2}+log(2\pi\sigma^2))]\] De modo que una parametrización para que esta distribución pertenezca a la familia exponencial queda definida como sigue: \[a(\phi)=\sigma^2;\hspace{0.5cm} \theta=\mu\] \[b(\theta)=\frac{\mu^2}{2}=\frac{\theta^2}{2}\] \[c(y;\phi)=-\frac{1}{2}(\frac{y^2}{\sigma^2}+log(2\pi\sigma^2))\]

Sea $N (,^2); con: y_i≥0; ∀i=1,…,n. $, se construye la función de máxima verosimilitud: \[\begin{align*} L(\mu,\sigma^2,y_i)&=\prod_{i=1}^nf(y_i;\mu;\sigma^2)\\ &=\prod_{i=1}^n[(2\pi\sigma^2)^{-1/2}exp(\frac{-\frac{1}{2}(y_i-\mu)^2}{\sigma^2})]\\ &=(2\pi\sigma^2)^{-n/2}exp(-\frac{1}{2\sigma^2}\sum_{i=1}^n(y_i-\mu)^2) \hspace{0.5cm}(1) \end{align*}\] Aplicando log a la expresión (1) tenemos que el log-verosimilitud es: \[\begin{align*} l(\mu,\sigma^2;y_i)&=log((2\pi\sigma^2)^{-n/2}exp(-\frac{1}{2\sigma^2}\sum_{i=1}^n(y_i-\mu)^2))\\ &=-\frac{n}{2}log(2\pi\sigma^2)-\frac{1}{2\sigma^2}\sum_{i=1}^n(y_i-\mu)^2\\ &=-\frac{n}{2}log(2\pi)-\frac{n}{2}log(\sigma^2)-\frac{1}{2\sigma^2}\sum_{i=1}^n(y_i-\mu)^2 \end{align*}\] Ahora hallaremos los estimadores:

\[\begin{align*} \frac{\partial l(\mu,\sigma^2,y_i)}{\partial\mu}&=\frac{\partial }{\partial\mu}(-\frac{n}{2}log(2\pi)-\frac{n}{2}log(\sigma^2)-\frac{1}{2\sigma^2}\sum_{i=1}^n(y_i-\mu)^2)\\ &=\frac{1}{\sigma^2}\sum_{i=1}^n(y_i-\mu)\\ &=\frac{1}{\sigma^2}(\sum_{i=1}^ny_i-n\mu)\\ \end{align*}\] Igualando a cero tenemos que: \[\frac{1}{\sigma^2}\sum_{i=1}^ny_i-n\mu=0\]

\[\frac{1}{\sigma^2}\sum_{i=1}^ny_i=n\mu\]

\[\widehat{\mu }=\frac{1}{n}\sum_{i=1}^{n}y_i\] estimador para $\mu$ Ahora derivemos con respecto a $\sigma^2$ \[\begin{align*} \frac{\partial l(\mu,\sigma^2,y_i)}{d\sigma^2}&=-\frac{n}{2\sigma^2}-[\frac{1}{2}\sum_{i=1}^{n}(y_i-\mu)^2](-\frac{1}{(\sigma^2)^2})\\ &=\frac{1}{2\sigma^2}+[\frac{1}{2}\sum_{i=1}^{n}(y_{i}-\hat(\mu))^2](\frac{1}{(\sigma^2)^2}) \end{align*}\] Igualando a cero tenemos que: \[0=\frac{1}{2\sigma^2}[\frac{1}{\sigma^2}\sum_{i=1}^{n}(y_i-\mu)^2-n]\]

\[\widehat{\sigma}=\frac{1}{n}\sum_{i=1}^{n}(y_i-\widehat{\mu})^2\] estimador para $\sigma$

Función Score: Sea $[\mu_0,\sigma_0^2]$

\[U=\begin{bmatrix}\frac{\partial}{\partial\mu}\\ \frac{\partial}{\partial\sigma^2}\end{bmatrix}=\begin{bmatrix}\frac{(y_i-\mu)}{\sigma^2}\\ \frac{1}{2\sigma^2}+\frac{1}{2}\frac{(y_i-\mu)^2}{(\sigma^2)^2}\end{bmatrix}\]

Ahora sea $[\hat{\mu_{i}},\hat\sigma_{i}]^T$ la matriz de información

\[I=\begin{bmatrix}\frac{\partial^2}{\partial\mu^2} & \frac{\partial^2}{\partial\mu\partial\sigma^2} \\ \frac{\partial^2}{\partial\mu\partial\sigma^2}& \frac{\partial^2}{\partial(\sigma^2)^2}\end{bmatrix}\]

\[\frac{\partial}{\partial\mu^2}(\frac{1}{\sigma^2}(y_i-\mu))=-\frac{1}{\sigma^2}\]

\[\frac{\partial}{\partial\sigma^2}(-\frac{1}{2\sigma^2}+\frac{1}{2}\frac{(y_i-\mu)^2}{\sigma^2})=\frac{1}{2(\sigma^2)^2}-\frac{(y_i-\mu)^2}{(\sigma^2)^3}\]

\[\frac{\partial}{\partial\sigma^2}(\frac{\partial}{\partial\mu}(-\frac{n}{2}log(2\pi)-\frac{n}{2}log(\sigma^2)-\frac{1}{2\sigma^2}(y_i-\mu)^2))=-\frac{y_i-\mu}{(\sigma^2)^2}\]

\[I=\begin{bmatrix}-\frac{1}{\sigma^2} & -\frac{(y_i-\mu)}{(\sigma^2)^2} \\ -\frac{(y_i-\mu)}{(\sigma^2)^2}& \frac{1}{2(\sigma^2)^2}-\frac{(y_i-\mu)^2}{(\sigma^2)^3}\end{bmatrix}\]

DISTRIBUCIÓN EXPONENCIAL

Sea $\mathbf{y}$ una variable aleatoria con distribución exponencial, es decir $\mathbf{y}\sim exp(\lambda)$ donde $\mathbf{y}=(y_{1},y_{2},...,y_{n})^T$ es una sucesión de variables aleatorias. La función de densidad para $\mathbf{y}$ corresponde a: \[(*)\hspace{0.5cm}f(\mathbf{y};\lambda)=\lambda e^{-\lambda \mathbf{y}};\hspace{0.5cm} \lambda,\mathbf{y}> \mathbf{0}\]

Primero veamos que la distribución exponencial pertenece a la familia exponencial. Para ello tendremos en cuenta la estructura de la familia exponencial que aparece en la definición 2

Para realizar la transformación apliquemos logaritmo a la función de densidad de la distribución exponencial definida en $(*)$.

\[\begin{align*} log f(\mathbf{y};\lambda)&=log\left(\lambda e^{-\lambda \mathbf{y}}\right)\\ &=log(\lambda)+log (e^{-\lambda \mathbf{y}})\\ &=-\mathbf{y}\lambda +log(\lambda) \end{align*}\] Luego, \[f(\mathbf{y};\lambda)=exp\left \{ -\mathbf{y}\lambda-(-log\lambda) \right \}.\] De modo que una parametrización, para que esta distribución pertenezca a la familia exponencial queda definida como sigue: \[a(\phi)=1;\hspace{0.5cm} \theta=-\lambda\] \[b(\theta)=-log\lambda=-log(-\theta)\] \[c(\mathbf{y};\lambda)=0\] Por otro lado, se debe cumplir que $E(\mathbf{y})=\frac{1}{\lambda}$;$Var(\mathbf{y})=\frac{1}{\lambda^2}$

De la definición 2, $E(\mathbf{y})$ y $Var(\mathbf{y})$ están dadas por: \[E(\mathbf{y})=a(\phi){b}'(\theta)\] \[Var(\mathbf{y})=a(\phi){b}''(\theta)\] Se puede verificar que: \[{b}'(\theta)=-\frac{1}{\theta};\hspace{0.5cm}{b}''(\theta)=\frac{1}{\theta^2}\hspace{0,5cm}(1)\]

Luego, reemplazando $a(\phi)$ y (1) en las ecuaciones para la esperanza y la varianza tenemos: \[E(\mathbf{y})=(1)\left ( \frac{-1}{\theta}\right )=\frac{-1}{\theta}=\frac{-1}{-\lambda}=\frac{1}{\lambda}\] \[Var(\mathbf{y})=(1)\left ( \frac{1}{\theta^2}\right )=\left ( \frac{1}{\theta}\right )\left ( \frac{1}{\theta}\right )=\frac{1}{-\lambda}\frac{1}{-\lambda}=\frac{1}{-\lambda^2}\] Sea $y_{i}\sim\exp(\lambda)$ con $y_i\geq 0$; $\forall$ $i=1,...,n$, y f.d.p. definida como $f(y_{i};\lambda)=\lambda e^{-\lambda y_{i}}$ se construye su función por máxima verosimilitud como sigue: \[\begin{align*} L(\theta,\mathbf{y})&=\prod_{i=1}^n f(y_i,\lambda)\\ &=\prod_{i=1}^n \lambda e^{-\lambda y_{i}}\\ &=\lambda^{n} e^{-\sum_{i=1}^n}\lambda y_i\\ l(\theta,\mathbf{y})&=log(\lambda^ne^{-\sum_{i=1}^n}\lambda y_i)\\ &=nlog\lambda-\sum_{i=1}^n\lambda y_i\\ &=nlog\lambda-\lambda\sum_{i=1}^ny_i\hspace{0.5cm}(2) \end{align*}\] El vector score se define como $U=\frac{\partial l(\textbf{y},\theta)}{\partial \theta}$, entonces procedemos a hallar la derivada parcial de la función log-verosimilitud de la distribución Poisson con respecto al parámetro $\lambda$: \[\frac{d l(\lambda,\mathbf{y})}{d\lambda}=\frac{n}{\lambda}-\sum_{i=1}^ny_i \hspace{0.5cm}(3)\] Así, $(3)$ es el vector score para una variable aleatoria con distribucón Poisson.

Para realizar la estimación del parámetro $\lambda$ por máxima verosimilitud igualamos a cero la expresión obtenida en (3) y realizamos las operaciones aritméticas concernientes, veamos: \[\begin{align*} 0&=\frac{n}{\lambda}-\sum_{i=1}^ny_i\\ \sum_{i=1}^ny_i&=\frac{n}{\lambda}\\ \lambda\sum_{i=1}^ny_i&=n\\ \frac{\sum_{i=1}^ny_i}{n}&=\frac{1}{\lambda}\\ \bar{y}&=\frac{1}{\lambda}\\ \hat{\lambda}&=\frac{1}{\bar{y}}\hspace{0.5cm}(4) \end{align*}\] Luego, (4) es una estimación por máxima verosimilitud del parametro $\lambda$.

Para calcular la matriz de información de Fisher, procedemos a calcular la derivada de segundo orden de la función log-verosimilitud $(2)$ con respecto al parámetro $\lambda$: \[I=-\frac{d^2l(\lambda,\mathbf{y})}{d\lambda^2}=-\frac{d}{d\lambda}\left \{\frac{n}{\lambda}-\sum_{i=1}^n y_i \right \}\] \[I=\frac{n}{\lambda^2}\]

DISTRIBUCIÓN WEIBULL

Sea $y$ una variable aleatoria con distribución Weibull con función de densidad de probabilidad como sigue: \[f(y;\alpha;\beta)=\frac{\alpha}{\beta^{\alpha}}y^{\alpha-1}e^{ \left(\frac{y}{\beta}\right)^{\alpha}}; con\hspace{0.5cm}\alpha,\beta,y>0.\] Donde $\beta$ es el parámetro de escala y $\alpha$ es el parámetro de forma.

Podemos reescribir la función de densidad de probabilidad Weibull de dos parámetros como sigue: \[\begin{align*} f(y;\alpha;\beta)&=exp\left \{log\left(\frac{\alpha}{\beta^{\alpha}}y^{\alpha-1}e^{ \left(\frac{y}{\beta}\right)^{\alpha}}\right) \right \}\\ &=exp\left \{log(\alpha)-\alpha log(\beta)+(\alpha-1)log(y)-\left(\frac{y}{\beta}\right)^{\alpha} \right \}\hspace{0.5cm}(1) \end{align*}\] Para que la distribución Weibull pertenezca a la familia exponencial uniparamétrica, consideremos el parámetro $\alpha$ como constante, es decir $\alpha=a$, con $a>0$, así $(1)$ se reescribe como: \[f(y;\beta)=log(a)-a log(\beta)+(a-1)log(y)-\left(\frac{y}{\beta}\right)^{a}\] De modo que una parametrización, para que esta distribución pertenezca a la familia exponencial queda definida como sigue: \[\theta=\beta; \hspace{0.5 cm}a(y)= y^a\] \[d(y)=(a-1)log(y);\hspace{0.5 cm} c(\theta)=log(a)-alog(\theta)\] Así \[f(y;\beta)= \left\{ \begin{array}{lc} e^{a(y)b(\theta)+c(\theta)+d(y)} & si & a < y < b \\ \\ 0 & & En \hspace{0.2cm}caso\hspace{0.2cm} contrario \end{array} \right.\]

DISTRIBUCIÓN CAUCHY

Sea $y$ una variable aleatoria continua de distribución Cauchy con parámetros $\alpha,\beta$ y función de densidad de probabilidad: \[f(y,\alpha,\beta)=\frac{1}{\pi}\left[\frac{\beta}{(y-\alpha)^2+\beta^2}\right],-\infty<y<\infty\] Donde $\alpha$ es el parámetro de localización y $\beta$ es el parámetro escala, la distribución de Cauchy uniparamétrica puede ser reescrita cuando $\alpha=0$ como sigue: \[f(y,\alpha,\beta)=\frac{1}{\pi}\left[\frac{\beta}{(\mathbf{y}^2+\beta^2}\right],-\infty<y<\infty\] La razón de por qué no podemos representar la distribución de Cauchy en la forma exponencial es que no podemos determinar $a(y)b(\theta)$ o cualquier otro término de la familia exponencial, es decir, no hay alguna expresión que nos permita reparametrizarla.

Representemos la distribución de Cauchy en la forma exponencial \[f(y,\alpha,\beta)=exp\left(log\left(\frac{1}{\pi}\left(\frac{\beta}{\mathbf{y}^2+\beta^2}\right)\right)\right)\] \[exp[log(\beta)-log\pi-log(y^2+\beta^2)]\] donde \[a(y)b(\theta)=-log(y^2+\beta^2)\] \[c(\theta)=log(\beta),\hspace{0.5cm}d(y)=-log(\pi)\]

DISTRIBUCIÓN UNIFORME

Sea $y$ una v.a. continua con distribución uniforme y paramétro $\theta$, con función de densidad de probabilidad como sigue: \[f(y;\theta)=\frac{1}{\theta}; 0<y<\theta\] Intentemos parametrizar esta función de densidad de probabilidad de modo que identifiquemos si pertenece o no a la familia exponencial:

\[f(y;\theta)=exp[-log(\theta)]\] De modo que una parametrización para esta distribución es: con \[a(y)=0; \hspace{0.5cm}b(\theta)=0\] \[d(y)=0;\hspace{0,5cm}c(\theta)=-log(\theta)\] Notamos que no podemos representar esta distribución como parte de la familia exponencial porque $y$ depende de $\theta$ y por la definición $y$ debe ser independiente del parámetro $\theta$.

4.5 EJEMPLO DE REGRESIÓN LINEAL SIMPLE PARA RESPUESTAS POISSON (DOBSON, 1990)

Los datos de la tabla 4.1 son recuentos Y, observados en varios valores de una covariable X. Están graficados en la Fig. 4.1.

Cuadro 4.5
Y	2	3	6	7	8	9	10	12	15
X	-1	-1	0	0	0	0	1	1	1

Análisis del modelo de regresión Poisson, el objetivo de la actividad fue replicar el algoritmo propuesto en el texto empleando el software Rstudio, el método desarrollado corresponde al método del Scoring, tomando como valores iniciales $\mathbf{b}=(b_{1},b_{2})^T=(7,5)^T.$

#Parametros iniciales
dt=as.matrix(dat)
colnames(dt)=NULL
Y=dt[,1]#(variablerespuesta)
uno=rep(1,each=1,times=9)
X=as.matrix(data.frame(uno,dt[,2])) #(Matriz de datos)
colnames(X)=NULL

#Funciones
f <- function(b){
  M4=X%*%b
  C2=as.vector(1/M4)
  W=as.matrix(diag(C2))
  Xt=t(X)
  XtW=Xt%*%W 
  XtWz=XtW%*%Y
  XtWX=XtW%*%X
  iXtWX=solve(XtWX)
  b=(iXtWX%*%XtWz)
  print(b)
}

#Solución inicial
b=as.matrix(c(7,5))
delta <- 1e-8 # tolerancia
iterM <- 10 # numero maximo de iteraciones
tolera <- 1 # inicializar tolera
m <- 0 # inicializar itera
lista<- b # inicializar historial de iteraciones

while((tolera>delta)&(m<iterM)){
  xold <- b# b viejo
    b <-f(b) #b nuevo
  tolera <- abs( b - xold )
  lista<- data.frame(lista,b);
  m<- m + 1
}

##          [,1]
## [1,] 7.451389
## [2,] 4.937500
##          [,1]
## [1,] 7.451632
## [2,] 4.935314
##          [,1]
## [1,] 7.451633
## [2,] 4.935300
##          [,1]
## [1,] 7.451633
## [2,] 4.935300

length(lista)

## [1] 5

SECCIÓN 2. EJERCICIOS DEL CAPÍTULO 3 (DOBSON, 1990)

EJERCICIO 3.1 DOBSON (1990)

Si la variable aleatoria $Y$ tiene la distribución Gamma con un parámetro de escala $\theta$ el cual es el parámetro de interés y un parámetro de forma conocido $\phi$, entonces su función de densidad de probabilidad es:

\[f(y;\alpha;\beta)=\frac{\beta^{\alpha}}{\Gamma(\alpha)}y^{\alpha-1}exp^{-\beta y}\] Demuestre que esta distribución pertenece a la familia exponencial y encuentre el parámetro natural. También usando los resultados de este capítulo, encuentre $E(Y)$ y $Var(Y)$.

Verifiquemos que la distribución Gamma pertenece a la familia exponencial. Para ello tengamos en cuenta la estructura de la familia exponencial que aparece en la definición 2.

EJERCICIO 3.2 (DOBSON, 1990)

Muestre que las siguientes distribuciones de probabilidad pertenecen a la familia exponencial:

Distribución Pareto $f(y;\theta)=\theta y^{-\theta-1}.$

Sea $\mathbf{y}$ una variable aleatoria con distribución Pareto, es decir, $\mathbf{y}\sim Pr(\alpha)$ y función de probabilidad como sigue: \[f(y,\alpha)=\alpha y^{-\alpha-1};\hspace{0.5cm} y,\alpha>0 \] Mostremos que pertenece a la familia exponencial, para ello utilizaremos la estructura descrita en la definición 1. Tenemos que: \[\begin{align*} f(y,\alpha)&=exp[log(\alpha y)^{-\alpha-1}]\\ &=exp[log(\alpha)-(\alpha+1)logy] \end{align*}\] De modo que una parametrización, para que esta distribución pertenezca a la familia exponencial queda definida como sigue: \[a(y)=logy; \hspace{0.5cm}\theta=(-\alpha-1)\] \[b(\theta)=-(\alpha+1)=(-\alpha-1)=\theta\] \[c(\theta)=log\alpha=log(\theta+1)\]

Distribución exponencial $f(\mathbf{y};\lambda)=\lambda e^{-\lambda \mathbf{y}}.$

Sea $\mathbf{y}$ una variable aleatoria con distribución exponencial, es decir, $\mathbf{y}\sim exp(\lambda)$ donde $\mathbf{y}=(y_{1},y_{2},...,y_{n})^T$ es una sucesión de variables aleatorias. La función de densidad para $\mathbf{y}$ corresponde a: \[(*)\hspace{0.5cm}f(\mathbf{y};\lambda)=\lambda e^{-\lambda \mathbf{y}};\hspace{0.5cm} \lambda,\mathbf{y}> \mathbf{0}\]

Primero veamos que la distribución exponencial pertenece a la familia exponencial. Para ello tendremos en cuenta la estructura de la familia exponencial que aparece en la definición 2.

Para realizar la transformación apliquemos logaritmo a la función de densidad de la distribución exponencial definida en $(*)$. \[\begin{align*} log f(\mathbf{y};\lambda)&=log\left(\lambda e^{-\lambda \mathbf{y}}\right)\\ &=log(\lambda)+log (e^{-\lambda \mathbf{y}})\\ &=-\mathbf{y}\lambda +log(\lambda) \end{align*}\] Luego, \[f(\mathbf{y};\lambda)=exp\left \{ -\mathbf{y}\lambda-(-log\lambda) \right \}.\] De modo que una parametrización para que esta distribución pertenezca a la familia exponencial queda definida como sigue: \[a(\phi)=1;\hspace{0.5cm} \theta=-\lambda\] \[b(\theta)=-log\lambda=-log(-\theta)\] \[c(\mathbf{y};\lambda)=0\]

Distribución binomial negativa $f(\textbf{y},\pi)=\begin{pmatrix} \textbf{y}+r-1\\r-1\end{pmatrix}\pi^r(1-\pi)^{\textbf{y}}.$

Inicialmente verifiquemos que esta distribución pertenece a la familia exponencial. Para ello utilizaremos la parametrización descrita en la definición 1:

Para realizar la transformación apliquemos logaritmo a la función de probabilidad de la distribución binomial negativa definida en $(1)$. \[\begin{align*} log(f(\mathbf{y},\pi))&=log\left \{ \begin{pmatrix} \mathbf{y}+r-1\\r-1\end{pmatrix}\pi^r(1-\pi)^{\mathbf{y}}\right \}\\ &= log\begin{pmatrix} \mathbf{y}+r-1\\ r-1\end{pmatrix}+log\pi^r(1-\pi)^{\mathbf{y}}\\ &=log\begin{pmatrix} \mathbf{y}+r-1\\ r-1\end{pmatrix} + r log\pi+\mathbf{y} log(1-\pi)\\ &= \mathbf{y} log(1-\pi)+(-rlog\pi)+log\begin{pmatrix} \mathbf{y}+r-1\\ r-1\end{pmatrix} \end{align*}\] Luego, \[f(\mathbf{y},\pi)=exp\left \{\mathbf{y} log(1-\pi)+(-rlog\pi)+log\begin{pmatrix} \mathbf{y}+r-1\\ r-1\end{pmatrix}\right \}\] De modo que una parametrización que hace que esta distribución pertenezca a la familia exponencial queda definida como sigue: \[\theta=log(1-\pi);\hspace{0.5cm}\pi=1-e^{\theta}\] \[a(\phi)=1;\hspace{0.5cm} c(\mathbf{y,\phi})=log\begin{pmatrix} y+r-1\\r-1 \end{pmatrix}\] \[b(\theta)=-r log\pi=-r log(1-e^{\theta})\]

EJERCICIO 3.3 (DOBSON 1990)

Para la distribucuión binomial muestre que $E(U)=0$ y $Var(U)=E(U^2)=E(-U')$ donde: $U=dl/d\theta$ y $l$ es la función log-verosimilitud

Sea $\mathbf{y}$ una variable aleatoria con distribución binomial, es decir: $\mathbf{y}\sim B(n,\pi);$ donde $\mathbf{y}=(y_{1},y_{2},...,y_{n})^t$ es una sucesión de variables aleatorias, como ya se ha mostrado en clase la función de probabilidad de una variable con distribución binomial pertenece a la familia exponencial con parametrización como sigue: \[\theta=log \left(\frac{\pi}{1-\pi}\right);\hspace{0.5cm} \pi=\frac{e^\theta}{1+e^\theta}; \hspace{0.5cm}\] \[b(\theta)=n\log(1-p)=-n\log(1+e^\theta)\] \[a(\phi)=1;\hspace{0.5cm}c(y,\phi)=log\binom{n}{y}\] Dada una muestra aleatoria $y_1,...,y_n$ cada una con distribución binomial, construimos la función de verosimilitud para $\mathbf{y}$: \[\begin{align*} L(\theta,\mathbf{y}))&=\prod_{i=1}^{n}f(y_i,\theta)\\ &=\prod_{i=1}^{n}\binom{n}{y_i}\pi^{y_i}(1-\pi)^{n-y_i} \end{align*}\] Por independencia en la sucesión de variables aleatorias la función de verosimilitud se reescribe como: \[\begin{align*} L(\theta,\mathbf{y})=\prod_{i=1}^{n}\binom{n}{y_i}\pi^{\sum_{i=1}^{n} y_i}(1-\pi)^{\sum_{i=1}^{n}(n-y_i)} \hspace{0.5cm}(1)\\ \end{align*}\] Aplicando logaritmo a la función de verosimilitud $(1)$ tenemos que el log-verosimilitud para la variable $\mathbf{y}$ es: \[\begin{align*} l(\theta,\phi,\mathbf{y})&=log\left(\prod_{i=1}^{n}\binom{n}{y_i}\pi^{\sum_{i=1}^{n} y_i}(1-\pi)^{\sum_{i=1}^{n}(n-y_i)}\right)\\ &=log\left(\prod_{i=1}^{n}\binom{n}{y_i}\right)+\sum_{i=1}^{n}y_ilog(\pi)+\sum_{i=1}^{n}(n-y_i)log(1-\pi) \end{align*}\] Por otro lado, el vector score se define como $U=\frac{\partial l(\textbf{y},\theta)}{\partial \theta}$, es decir, la derivada parcial de primer orden de función log-verosimilitud con respecto al vector de parámetros $\theta$, para nuestro ejercicio esto sería equivalente a la derivada parcial de primer orden de función log-verosimilitud de la distribución binomial $(1)$ con respecto al parámetro $\pi$: \[\begin{align*} U&=\frac{\partial l(\textbf{y},\pi)}{\partial \pi}\\ &=\frac{1}{\pi}\sum_{i=1}^{n}y_i-\frac{1}{(1-\pi)}\sum_{i=1}^{n}(n-y_i)\\ &=\frac{1}{\pi}\sum_{i=1}^{n}y_i-\frac{\sum_{i=1}^{n}(n-y_i)}{(1-\pi)} \end{align*}\] Verifiquemos que en efecto $E(U)=0$ \[E(U)=E\left(\frac{1}{\pi}\sum_{i=1}^ny_i- \frac{\sum_{i=1}^n(n-y_i)}{n-\pi}\right)\hspace{0.5cm}(2)\] Como $y_i\sim B(n,\pi)$ para $i=1,...,n$; se cumple que $E(y_i)=n\pi;\hspace{0.5cm} Var(y_i)=n\pi(1-\pi).$ Aplicando las propiedades algebraicas de la esperanza en $(2)$ y al ser $y_i$ para $i=1,\cdots,n$ variables independientes, tenemos: \[\begin{align*} E(U)&=E\left(\frac{1}{\pi}\sum_{i=1}^ny_i\right)-E\left(\frac{\sum_{i=1}^n n}{1-\pi}\right)+ E\left(\frac{\sum_{i=1}^ny_i}{n-\pi}\right)\\ &=\frac{1}{\pi}\sum_{i=1}^{n}E(y_i)-\frac{n^2}{1-\pi}+\frac{1}{1-\pi}\sum_{i=1}^nE(y_i)\\ &=\frac{1}{\pi}\sum_{i=1}^{n}n\pi-\frac{n^2}{1-\pi}+\frac{1}{1-\pi}\sum_{i=1}^n n\pi\\ &=\frac{n^2\pi}{\pi}-\frac{n^2}{1-\pi}+\frac{n^2\pi}{1-\pi}\\ &=\frac{n^2(1-\pi)}{1-\pi}-\frac{n^2}{1-\pi}+\frac{n^2\pi}{1-\pi}\\ &=\frac{n^2-n^2\pi-n^2+n^2\pi}{(1-\pi)}\\ \therefore E(U)&=0 \end{align*}\] Dado una v.a. $y$ única, tal que $y\sim B(n,\pi)$ con $E(y)=n\pi$ y $V(y)=n\pi(1-\pi)$ se cumple que la función de log-verosimilitud para $y$ es: \[l(\theta,y)=log\binom{n}{y}+ylog(\pi)+(n-y)log(1-\pi)\] y su vector score para $y$ corresponde a: \[U=\frac{\partial l(\pi,y)}{\partial\pi}=\frac{y}{\pi}-\frac{n-y}{1-\pi}\] Notamos a continuación que se sigue cumpliendo el hecho que, $E(U)=0$ \[\begin{align*} E(U)&=E\left(\frac{y}{\pi}-\frac{n-y}{1-\pi}\right)\\ &=\frac{1}{\pi}E(y)-\frac{1}{(1-\pi)}E(n-y)\\ &=\frac{1}{\pi}n\pi-\frac{n}{1-\pi}+\frac{n\pi}{(1-\pi)}\\ E(U)&=\frac{n(1-\pi)-n+n\pi}{(1-\pi)}\\ \therefore E(U)=0\hspace{0.5cm}(*) \end{align*}\] Por otro lado tenemos que: \[\begin{align*} U^2&=\left(\frac{y}{\pi}-\frac{n-y}{1-\pi}\right)^2\\ &=\frac{(y-y\pi-n\pi+y\pi)^2}{\pi^2(1-\pi)^2}\\ &=\frac{y^2-2ny\pi+n^2\pi^2}{\pi^2(1-\pi)^2}\hspace{0.5cm}(3) \end{align*}\] Así haciendo $E(U^2)$ donde $U$ es como se definió en $(3)$, y teniendo en cuenta que $E(y^2)=n^2\pi^2+n\pi(n-\pi)$ obtenemos lo siguiente: \[\begin{align*} E(U^2)&=E\left(\frac{y^2-2ny\pi+n^2\pi^2}{\pi^2(1-\pi)^2}\right)\\ &=\frac{1}{\pi^2(1-\pi^2)}E(y^2)-\frac{2n\pi E(y)}{\pi^2(1-\pi)^2}+\frac{n^2\pi^2}{\pi^2(1-\pi)^2}\\ &=\frac{n^2\pi^2+n\pi(1-\pi)}{\pi^2(1-\pi)^2}-\frac{2n^2\pi^2}{\pi^2(1-\pi)^2}+\frac{n^2\pi^2}{\pi^2(1-\pi)^2}\\ &=\frac{n\pi(1-\pi)}{\pi^2(1-\pi)^2}\\ E(U^2)&=\frac{n}{\pi(1-\pi)}\hspace{0.5cm}(**) \end{align*}\] Por propiedades de la varianza tenemos $Var(U)=E(U^2)-(E(U))^2$, reemplazando $(*)$ y $(**)$ en esta expresión llegamos a que: \[Var(U)=\frac{n}{\pi(1-\pi)}-0^2=\frac{n}{\pi(1-\pi)}\hspace{0.5cm}(4)\] Calculando las derivadas parciales de segundo orden de l(, y) con respecto a $\pi$, llegamos a una expresión equivalente a $U'$ : \[U'=\frac{\partial^2l(y,\pi)}{\partial\pi^2}=\frac{y}{\pi^2}-\frac{(n-y)}{(1-\pi)^2}\hspace{0.5cm}(5)\] Determinando la $E(-U')$ como sigue, vemos que: \[\begin{align*} E(-U')=-E(U')&=-E(-\frac{y}{\pi^2}-\frac{(n-y)}{(1-\pi)^2})\\ &=\frac{E(y)}{\pi^2}+E(\frac{n}{(1-\pi)^2}-\frac{E(y)}{(1-\pi)^2}\\ &=\frac{n\pi}{\pi^2}+\frac{n}{(1-\pi)^2}-\frac{n\pi}{(1-\pi)^2}\\ &=\frac{n(1-\pi)^2+n\pi-n\pi^2}{\pi(1-\pi)^2}\\ &=n-2n\pi+n\pi^2+n\pi-n\pi^2;\hspace{0.5cm}luego\\ &=\frac{n-n\pi}{\pi(1-\pi)^2}=\frac{n(1-\pi)}{\pi(1-\pi)^2}\\ E(-U')&=\frac{n}{\pi(1-\pi)}\hspace{0.5cm}(6) \end{align*}\] Entonces por $(**)$, $(4)$ y $(6)$ queda mostrado que $Var(U)=E(U^2)=E(-U')$.

EJERCICIO 3.4 (DOBSON, 1990) Utilice las ecuaciones (3.4) y (3.5) para verificar estos resultados:

E(y)=V(y) distribución de Poisson

Recordemos las ecuaciones 3.4 y 3.5 del (Dobson, 1990)

\[E[a(y)]=-\frac{c'(\theta)}{b'(\theta)} \hspace{0.5cm} (3.4)\]

\[V[a(y)]=\frac{b''(\theta)c'(\theta)-b'(\theta)c''(\theta)}{[b'(\theta)]^3} \hspace{0.5cm} (3.5)\] La Función de probabilidad de Poisson está dada por:

\[f(y;\lambda)=\frac{\lambda^y\ exp(-\lambda)}{y!};\hspace{0.5cm} Y\sim P(\lambda);\hspace{0.5cm}y,\lambda>0\] Primero verifiquemos que la distribución de Poisson pertenece a la familia exponencial utilizando la siguiente estructura:

\[f(y;\theta)=exp[a(y)b(\theta)+c(\theta)+d(y)]\] \[\begin{align*} \ f(y;\lambda)&=Exp\left[log\left(\frac{\lambda^y\ exp(-\lambda)}{y!}\right)\right]\\ &=Exp[y\ log(\lambda)-log(y!)-\lambda]\\ &=Exp[y\ log(\lambda)-\lambda-log(y!)]\hspace{0.5cm}(*) \\ \end{align*}\]

De la expresión (*) tenemos que:

\[a(y)=y\] \[b(\theta)=log(\lambda)\] \[c(\theta)=-\lambda\] \[d(y)=-log(y!)\] Calculemos la esperanza utilizando la fórmula (3.4)

\[E[a(y)]=-\frac{c'(\theta)}{b'(\theta)}, c'(\theta)=-1,b'(\theta)=1/\lambda \] \[E(y)=-\frac{-(-1)}{1/\lambda}=\frac{1}{1/\lambda}=\lambda\]

Ahora verifiquemos la varianza utilizando la ecuación (3.5)

\[V[a(y)]=\frac{b''(\theta)c'(\theta)-b'(\theta)c''(\theta)}{[b'(\theta)]^3} \ (3.5)\] Del paso anterior conocemos

\[c'(\theta)=-1,b'(\theta)=1/\lambda, a(y)=y \] reemplazamos en (3.5)

\[b''(\theta)=-1/\lambda^2, c''(\theta)=0\] tenemos que:

\[V(y)=\frac{(-1/\lambda^2)(-1)-(0)(1/\lambda)}{[1/\lambda]^3}=\frac{\frac{1}{\lambda^2}}{\frac{1}{\lambda^3}}=\frac{\lambda^3}{\lambda^2}=\lambda\] Por lo tanto $E[y]=Var[y]$

Si \[y\sim N (\mu,\sigma^2),E(y)=\mu,Var(y)=\sigma^2\] La función de densidad de la distribución normal está dada por:

\[f(y;\mu)=\frac{1}{(2\pi\sigma^2)^{1/2}}exp\left(-\frac{1}{2\sigma^2}(y-\mu)^2\right); \hspace{0.5cm} \mu,\sigma,y>0 \] Expresemos la normal en forma de familia exponencial utilizando la parametrización:

\[f(y;\theta)=exp[a(y)b(\theta)+c(\theta)+d(y)]\] \[\begin{align*} f(y;\mu)&=Exp\left[log\left(\frac{1}{(2\pi\sigma^2)^{1/2}}exp\left(-\frac{1}{2\sigma^2}(y-\mu)^2\right)\right)\right]\\ &=Exp\left[log(1)-\frac{1}{2}log(2\pi\sigma^2)-\frac{y^2}{2\sigma^2}+\frac{y\mu}{\sigma^2}-\frac{\mu^2}{2\sigma^2}\right]\\ &=exp\left[-\frac{1}{2}log(2\pi\sigma^2)-\frac{y^2}{2\sigma^2}+\frac{y\mu}{\sigma^2}-\frac{\mu^2}{2\sigma^2}\right] \end{align*}\]

De la expresión anterior, tenemos que: \[a(y)=y;\hspace{0.5cm}b(\theta)=\frac{\mu}{\sigma^2};\hspace{0.5cm}c(\theta)=\frac{-\mu^2}{2\sigma^2};\hspace{0.5cm}d(y)=\frac{-y^2}{2\sigma^2}-\frac{1}{2}log(2\pi\sigma^2);\hspace{0.5cm}y\hspace{0.5cm}\theta=\mu \] Calculemos la esperanza y la varianza utilizando las ecuaciones (3.4) y (3.5), respectivamente:

\[\begin{align*} E[Y]&=-\frac{\frac{-\mu}{\sigma^2}}{\frac{1}{\sigma^2}}\\ &=\mu \end{align*}\] \[E[Y]=\mu\] Luego,

\[\begin{align*} Var(y)&=\frac{\frac{1}{\sigma^4}}{\frac{1}{\sigma^6}}\\ &=\frac{\sigma^6}{\sigma^4}\\ &=\sigma^2 \end{align*}\]

\[var[Y]=\sigma^2\] C. \[y_i\sim b(n,\pi),\hspace{0.5cm} E(Y)=n\pi,\hspace{0.5cm} Var(y)=n\pi(1-\pi)\]

La función de probabilidad de la distribución binomial está dada por:

\[f(y;\pi)={n\choose y}\ \pi^y(1-\pi)^{n-y}\]

Expresemos la distribución como familia exponencial utilizando la parametrización dada en la definición 1:

\[\begin{align*} f(y;\pi)&=exp\left[log\left({n \choose y }\ \pi^y(1-\pi)^{n-y}\right)\right]\\ &=exp[log({n \choose y })+ylog(\pi)+(n-y)log(1-\pi))]\\ &=exp[ylog(\pi)+(n-y)log(1-\pi))+log({n \choose y })]\\ &=exp[ylog(\pi)+n\ log(1-\pi))-y\ log(1-\pi))+log({n \choose y })]\\ &=exp[y\ log\left(\frac{\pi}{1-\pi}\right)+n\ log(1-\pi))+log({n \choose y })]\hspace{0.5cm}(**) \end{align*} \]

De la expresión (**), tenemos que:

\[c(\theta)=n\ log(1-\pi)\] \[d(y)=log\left({n \choose y }\right)\] \[a(y)=y\] \[b(\theta)=log\left(\frac{\pi}{1-\pi}\right)\]

\[\theta=\pi\] Calculemos la esperanza utilizando la ecuación (3.4) \[c'(\theta)=\frac{n}{1-\pi}(-1)=-\frac{n}{1-\pi}\] \[b'(\theta)=\frac{1}{\frac{\pi}{1-\pi}}\frac{1}{(1-\pi)^2}=\frac{1-\pi}{\pi}\frac{1}{(1-\pi)^2}=\frac{1}{\pi(1-\pi)}\] \[a(y)=y\] \[E[y]=-\left(\frac{-\frac{n}{1-\pi}}{\frac{1}{\pi(1-\pi)}}\right)=\frac{n\pi(1-\pi)}{1-\pi}=n\pi\rightarrow E[Y]=n\pi\] Calculemos la varianza utilizando la ecuación (3.5)

Del paso anterior conocemos:

\[b'(\theta)=\frac{1}{\pi(1-\pi)}, c'(\theta)=-\frac{n}{1-\pi}\] \[b''(\theta)=\frac{-(1-\pi)+\pi}{\pi^2(1-\pi)^2}=\frac{-1+2\pi}{\pi^2(1-\pi)^2}\] \[c''(\theta)=-\frac{n}{(1-\pi)^2}\] Reemplazamos en (3.5)

\[\begin{align*} V[y]&=\frac{\left(\frac{-1+2\pi}{\pi^2(1-\pi)^2}\right)\left(-\frac{n}{1-\pi}\right)-\left(-\frac{n}{(1-\pi)^2}\right)\left(\frac{1}{\pi(1-\pi)}\right)}{\left[\frac{1}{\pi(1-\pi)}\right]^3}\\ &=\frac{\frac{n(1-2\pi)}{\pi^2(1-\pi)^3}+\frac{n}{\pi(1-\pi)^3}}{\frac{1}{\pi^3(1-\pi)^3}}\\ &=\frac{\frac{n}{\pi(1-\pi)^3}\left(\frac{1-2\pi}{\pi}+1\right)}{\frac{1}{\pi^3(1-\pi)^3}}\\ &=\left(\frac{n\pi^3(1-\pi)^3}{\pi(1-\pi)^3}\right)\left(\frac{1-\pi}{\pi}\right)\\ &=\frac{n\pi^2(1-\pi)}{\pi}=n\pi(1-\pi) \end{align*}\]

EJERCICIO 3.5 (DOBSON, 199O)

Grafique los datos de la Tabla 3.2 en una escala adecuada para que pueda estimar aproximadamente el parámetro $\theta$ en el modelo propuesto con $\lambda_{i}=i^{\theta}$ Utilice este valor para estimar los valores esperados $E(y_{i})$ para cada $i$ compare las estimaciones con los valores observados $y_{i}$ ¿Este modelo parece ajustarse a los datos?

#Gráfico de los datos sin escala adecuada y con escala adecuada

i<-c(1,2,3,4,5,6,7,8,9,10,11,12,13,14)
yi<-c(0,1,2,3,1,4,9,18,23,31,20,25,37,45)
logi<-log(i)
par(mfrow=c(1,2))
plot(x=logi,y=yi,main="Gráfico con ajuste de escala",xlab="Logi",ylab="yi",pch=9,ylim = c(0,50),cex.main=1.0,col="blue",xlim = c(0,5))
plot(x=i,y=yi,main="Gráfico sin ajuste de escala",xlab="i",ylab="yi",pch=9,ylim =c(0,50),cex.main=1.0,col="red",xlim = c(0,15))

Claramente, el número de muertes aumenta con $i$. Para estos datos, un modelo posible es la distribución de Poisson con $\lambda_{i}=i^{\theta}$, donde $\theta$ es un parámetro a estimar. Esto puede describirse mediante un modelo lineal generalizado en el que la función de enlace es: \[g(\lambda_{i})=log(\lambda_{i})=\theta logi\] De modo que $\textbf{x}_{i}=[logi]$, $\beta=[\theta]$

Modelo de regresión lineal

i<-c(2,3,4,5,6,7,8,9,10,11,12,13,14)
yi<-c(1,2,3,1,4,9,18,23,31,20,25,37,45)
logi<-log(i)
logyi=log(yi)
r_Model=lm(logyi~logi)
summary(r_Model)

## 
## Call:
## lm(formula = logyi ~ logi)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.48957 -0.13769  0.08492  0.39243  0.46623 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  -1.9427     0.5395  -3.601  0.00416 ** 
## logi          2.1326     0.2669   7.991  6.6e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.5536 on 11 degrees of freedom
## Multiple R-squared:  0.853,  Adjusted R-squared:  0.8397 
## F-statistic: 63.85 on 1 and 11 DF,  p-value: 6.603e-06

Se concluye que a un nivel de significancia del 5% el parámetro para “logi” y el término constante son significativos, con un p-valor=6.6e-06 y 0.00416, respectivamente. En este caso, concluimos que el modelo es adecuado.

#Modelo tomando las variables yi e logi

Modelo=glm(yi~logi,family = poisson(link="log"))                              
summary(Modelo)

## 
## Call:
## glm(formula = yi ~ logi, family = poisson(link = "log"))
## 
## Deviance Residuals: 
##      Min        1Q    Median        3Q       Max  
## -2.11278  -1.24843   0.01996   0.38310   1.93539  
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  -1.9055     0.5186  -3.674 0.000239 ***
## logi          2.1587     0.2180   9.900  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 174.81  on 12  degrees of freedom
## Residual deviance:  16.80  on 11  degrees of freedom
## AIC: 73.727
## 
## Number of Fisher Scoring iterations: 4

Se concluye que a un nivel de significancia del 5% el parámetro para “logi” y el término constante son significativos con un p-valor=2e-16 y 0.000145, respectivamente. En este caso, concluimos que un modelo Poisson es adecuado dado que realiza un buen ajuste a los datos, hecho que se puede evidenciar gráficamente, con los valores predichos.

EJERCICIO 3.6 (DOBSON, 1990)

Sea $y_1,....y_n$ n variables aleatorias binarias e independientes. Además, $y_i\sim \beta(N,\pi_i)\ donde \ i=1,...,N$ \[y_i= 1\hspace{0.5cm} con\ probabilidad\ \pi_i\] \[y_i=0\hspace{0.5cm} con\ probabilidad\ 1-\pi_i\] La función de probabilidad de densidad es

\[f(y_i,\pi_i)=\pi_i^{y_i}(1-\pi_i)^{1-y_{i}},\hspace{0.5cm} y_i=0\ ó\ y_i=1\]

La estructura de la fórmula exponencial viene dada por

\[f(y;\theta;\phi)=exp\left[\frac{y\theta-b(\theta)}{a(\phi)}+C(y,\phi)\right]\] \[\begin{align*} \Rightarrow log\ f(y_i,\pi_i)&=log\ (\pi_i^{y_i}(1-\pi_i)^{1-y_i})\\ &=y_{i}\ log\ \pi_i+(1-y_i)log(1-\pi_i) \end{align*}\]

\[\begin{align*}f(y_i,\pi_i)&=exp\ [{{y_{i}\ log\ \pi_i+(1-y_i)log(1-\pi_i)}}]\\ &=exp [y_ilog\pi_i-y_ilog(1-\pi_i)+log(1-\pi_i)]\\ &=exp[y_ilog(\frac{\pi_i}{1-\pi_i})+log(1-\pi_i)] \end{align*}\]

De la ecuación anterior:

\[\theta=log \left(\frac{\pi_i}{1-\pi_i}\right),\hspace{0.5cm} C(y,\phi)=0\] \[b(\theta)=-log(1-\pi_i)\] \[a(\phi)=1\] \[\theta=log\left(\frac{\pi_i}{1-\pi_i}\right)\] \[e^\theta=\left(\frac{\pi_i}{1-\pi_i}\right)\rightarrow\ e^\theta(1-\pi_i)=\pi_i\] \[e^\theta-e^\theta\pi_i=\pi_i\] \[e^\theta=\pi_i+e^\theta\pi_i\] \[e^\theta=(1+e^\theta)\pi_i\] \[\pi_i=\frac{e^\theta}{1+e^\theta}\] \[\begin{align*} b(\theta)&=-log(1-\pi_i)\\ &=-log\left(1-\frac{e^\theta}{1+e^\theta}\right)\\ &=-log\left(\frac{1+e^\theta-e^\theta}{1+e^\theta}\right)=-log\left(\frac{1}{1+e^\theta}\right)\\ &=-log(1)+log(1+e^\theta)=log(1+e^\theta) \end{align*}\]

\[b(\theta)=log(1+e^\theta)\] b. Sea

\[g(\pi_i)=\sum_{i=1}^{n}\beta_iy_i\] \[g(\pi_i)=log\left(\frac{\pi_i}{1-\pi_i}\right)\]

La función log mapea el intervalo $(0,1)$. Podemos escribir la probabilidad de éxito como:

\[\frac{\pi_i}{1-\pi_i}=e^{Y_i\beta_i}\] \[\pi_i=e^{Y_i\beta_i}(1-\pi_i)\] \[\pi_i=e^{Y_i\beta_i}-\pi_ie^{Y_i\beta_i}\] \[\pi_i+\pi_ie^{Y_i\beta_i}=e^{Y_i\beta_i}\] \[\pi_i(1+e^{Y_i\beta_i})=e^{Y_i\beta_i}\] \[\pi_i=\frac{e^{Y_i\beta_i}}{1+e^{Y_i\beta_i}}\]

Como $\pi_i$ es una probabilidad, entonces

\[log\left(\frac{\pi_i}{1-\pi_i}\right)\] es el log(odds ratio)

Donde odds ratio de un evento $A=\frac{p_r(A)}{1-p_r(A)}$

\[E(y_i)=\pi_i,\hspace{0.5cm} b(\theta)=log(1+e^\theta)\]

\[\begin{align*} E(y_i)&=a(\phi)b'(\theta)\\ &=(1)b'(\theta)=b'(\theta)\\ &=\frac{1}{1+e^\theta}\frac{\partial }{\partial \theta}(1+e^\theta)\\ &=\frac{e^\theta}{1+e^\theta} \end{align*}\]

Sabemos que $\theta=log\left(\frac{\pi_i}{1-\pi_i}\right)$

\[\rightarrow\frac{exp\ \left[log\left(\frac{\pi_i}{1-\pi_i}\right)\right]}{1+exp\ \left[log\left(\frac{\pi_i}{1-\pi_i}\right)\right]}=\frac{\frac{\pi_i}{1-\pi_i}}{1+\frac{\pi_i}{1-\pi_i}}\] \[=\frac{\frac{\pi_i}{1-\pi_i}}{\frac{1+\pi_i-\pi_i}{1-\pi_i}}=\frac{\frac{\pi_i}{1-\pi_i}}{\frac{1}{1-\pi_i}}=\frac{\pi_i(1-\pi_i)}{1-\pi_i}=\pi_i\]

\[\ si\ g(\pi)=log(\frac{\pi}{1-\pi})=X^{T}\beta\]

\[log\left(\frac{\pi}{1-\pi}\right)=X^{T}\beta\] \[exp\ \left[log\left(\frac{\pi}{1-\pi}\right)\right]=exp(X^{T}\beta)\] \[\left(\frac{\pi}{1-\pi}\right)=exp(X^{T}\beta)\] \[\pi=exp(X^{T}\beta)(1-\pi)\] \[\pi=exp(X^{T}\beta)-\pi \exp(X^{T}\beta)\] \[\pi(1-exp(X^{T}\beta))=exp(X^{T}\beta)\] \[\pi=\frac{exp(X^{T}\beta)}{(1-exp(X^{T}\beta))}\] e. Sean $\beta_1,\beta_2$ constantes, consideremos el caso particular en el que $\beta_1=1.5,\hspace{0.5cm} \beta_2=2.5$ Teniendo en cuenta que: \[\pi=\frac{exp(\beta_1+\beta_2x)}{1-exp(\beta_1+\beta_2x)} \] Obtenemos los siguientes valores:

X	0.2000000	0.5000000	0.750000	1.0000000	2.0000000	3.0000000	4.0000000	5.0000000	6.0000000	7	8	9	10	11	12	13
Pi	0.8807971	0.9399134	0.966914	0.9820138	0.9984988	0.9998766	0.9999899	0.9999992	0.9999999	1	1	1	1	1	1	1

De modo que la gráfica queda definida como sigue:

En donde notamos que a medida que la dosis, (x), aumenta la probabilidad de muerte (π) aumenta desde casi cero hasta un valor asintótico de 1.

EJERCICIO 3.7 (DOBSON, 1990)

Es la distribución de valor extremo (Gumbel) \[f(y;\theta)=\frac{1}{\theta}exp\left \{ \frac{x-\theta}{\phi}-exp\frac{x-\theta}{\phi} \right \} \hspace{1cm}(1)\] Donde $\phi>0$ es considerado como un parámetro molestia es esta distribución un miembro de la familia exponencial.

Sea $\textbf{y}$ una variable aleatoria con distribución de valor extremo (Gumbel) y parámetro $\phi>0$, cuya función de densidad de probabilidad corresponde a (1). Verifiquemos que esta distribución pertenece a la familia exponencial. Para ello tengamos en cuenta la estructura de la familia exponencial que aparece en la definición 1.

Para realizar la transformación apliquemos log a la función de densidad de la distribución Rayleigh. \[\begin{align*} log f(\mathbf{y};\phi)&=log \left(\frac{1}{\theta} exp^{{\frac{y-\theta}{\phi}}-exp^{{\frac{y-\theta}{\phi}}}}\right)\\ &=log\left(\frac{1}{\theta}\right)+log( exp^{{\frac{y-\theta}{\phi}}-exp^{{\frac{y-\theta}{\phi}}}})\\ &=log(\frac{1}{\theta}+\frac{(y-\theta)}{\phi}-exp^{{\frac{(y-\theta)}{\phi}}})\\ &=\frac{1}{\theta}+\frac{y}{\phi}-\frac{\theta}{\phi}-exp^{{\frac{y}{\theta}}}exp^{{-\frac{\theta}{\phi}}}\hspace{0.5cm}(1) \end{align*}\] Luego, \[f(y;\theta)=exp\left(log\left(\frac{1}{\theta}\right)+\frac{y}{\phi}-\frac{\theta}{\phi}-exp^{{\frac{y}{\theta}}}exp^{{-\frac{\theta}{\phi}}}\right)\] De modo que una parametrización para que esta distribución pertenezca a la familia exponencial queda definida como sigue: \[=exp \left \{-exp^{{\frac{y}{\phi}}}exp^{{-\frac{\theta}{\phi}}}-log(\theta)-\frac{\theta}{\phi}+\frac{y}{\phi}\right \}\] Donde: \[a(y)=-exp^{\frac{y}{\phi}}\] \[b(\theta)=exp^{-\frac{\theta}{\phi}}\] \[c(\theta)=-log(\theta)-\frac{\theta}{\phi}\] \[d(y)=\frac{y}{\phi}\] EJERCICIO 3.8 DOBSON (1990) Sea $y_{1},\ldotp\ldotp\ldotp,y_{n}$ variables aleatorias independientes $Y_{i}=\beta_{0}+log(\beta{1}+\beta_{2}X_{i})+e_{i}$, donde $e_{i}~N(0,\sigma^{2})$, para todo $i=1,\ldotp\ldotp\ldotp,n$, ¿Es esto un modelo lineal generalizado? Justifica la respuesta.

Asumamos que lo descrito anteriormente es un modelo lineal generalizado, entonces: Sea \[\begin{align*} Y_{1}=\beta_{0}+log(\beta_{1}+\beta_{2}X_{1})\\ Y_{2}=\beta_{0}+log(\beta_{1}+\beta_{2}X_{2})\\ \end{align*}\] \[\vdots\\\] \[Y_{n}=\beta_{0}+log(\beta_{1}+\beta_{2}X_{n})\] Ya demostramos previamente que la distribución normal pertenece a la familia exponencial.

Tomemos la función exponencial como la función de enlace, entonces:

\[exp(\mu_{i})=exp[\beta_{0}+ log(\beta_{1}+\beta_{2}X_{n}]\]

\[=\beta_0^{*}(\beta_0+\beta_{2}x_{i})\] \[=\beta_{1}^{*}+\beta_2^{*}x_{i}\] Por definición, sabemos que $g(\mu_{i})=X^t(\beta)$ lo que no es posible obtener ya que el logaritmo natural no es un operador lineal por lo cual no podemos obtener $X^t$ como combinación lineal de $(\mathbf{\beta})$ por lo que el modelo definido anteriormente no es un modelo lineal generalizado.

SECCIÓN 4. EJERCICIOS CAPÍTULO 4 (DOBSON, 1990)

EJERCICIO 4.1 (DOBSON, 1990)

Los datos del Cuadro 4.1 muestran el número de muertes por SIDA en Australia durante periodos sucesivos de tres meses desde 1983 a 1986.

Cuadro 4.1 Número de muertes por SIDA en Australia por trimestre desde enero- marzo de 1983 y abril-junio de 1986; yi denota el número de muertes $x_{i} = logi$ donde i = 1, …, 14. indica el trimestre
yi	0	1.000	2.000	3.000	1.000	4.000	9.000	18.000	23.000	31.000	20.000	25.000	37.000	45.000
logi	0	0.693	1.099	1.386	1.609	1.792	1.946	2.079	2.197	2.303	2.398	2.485	2.565	2.639

Suponga que las variables aleatorias $y_{i}$, son variables de Poisson con $E (y_{i})=\mu_{i}; con \hspace{0.5cm} i=1,\cdots,14.$, Donde: $g(\mu_{i})=log\mu_{i}=\beta{1}+\beta_{2}x_{i}$

y $x_{i}=logi$. La función de enlace utilizada en este caso es la función logarítmica (que es el enlace ‘natural’ para la distribución de Poisson en el sentido de que corresponde al parámetro).

Utilice las ecuaciones (4.9) y (4.10) para obtener expresiones para los elementos de $W$ y $z$ para este modelo.

Sea $\textbf{y}\sim P(\mu)$ con $\textbf{y}=(y_{1},y_{2},...,y_{n})^T$ una sucesión de variables aleatorias para $i=1,\cdots,n$ se define la distribución de probabilidad como sigue: $f(\textbf{y},\mu)=\frac{e^{-\mu}\mu^{\textbf{y}}}{\textbf{y}!}; con\hspace{0.5cm}\textbf{y}>0;\mu=\lambda n$

donde:

$E(y_{i})=Var(y_{i})=\mu_{i}$

Para este ejercicio tenemos que $x_{i}=logi$, con función de enlance igual a log, en efecto $g(\mu_{i})=log\mu_{i}=\eta_{i}$, el modelo se define como:

\[log\mu_{i}=\eta_{i}=\beta_{1}+\beta_{2}x_{i}\] \[e^{log\mu_{i}}=e^{ \eta_{i}} = e^{\beta_{1}+\beta_{2}x_{i}}\]

\[\mu_{i}=e^{ \eta_{i}}=e^{\beta_{1}+\beta_{2}x_{i}}\hspace{0.2cm}(1) \]

De la ecuación (4.9) (Dobson, 1997), tenemos que $W$ es una matriz diagonal de orden $14\times14$ cuya estructura es: \[W_{ii}=\frac{1}{Var(\textbf{y}_{i})} \left( \frac{\partial\mu_{i}}{\partial\eta_{i}} \right)^2\hspace{0.5cm} (*)\] para nuestro ejercicio
\[\frac{\partial\mu_{i}}{\partial\eta_{i}}=e^{\eta_{i}}=e^{\beta_{1}+\beta_{2}x_{i}}\hspace{0.5cm} (2)\] Reemplazando (1) y (2) en (*) llegamos a que $W_{ii}$equivale a la siguiente expresión: \[W_{ii}=\frac{1}{\mu_{i}}\left(e^{\beta_{1}+\beta_{2}x_{i}}\right)^2\] \[W_{ii}=\frac{\left(e^{\beta_{1}+\beta_{2}x_{i}}\right)^2}{e^{\beta_{1}+\beta_{2}x_{i}}}\] \[W_{ii}=e^{\beta_{1}+\beta_{2}x_{i}}\] \[W_{ii}=[\mu_{i}]_{\beta= \hat{b}}\] Siendo $\hat{b}$ el estimador por máxima verosimilitud de $\beta$.

En general la ecuación de $z_{i}$ equivale a: \[z_i=\underbrace{\sum_{k}x_{ik}b_{k}^{(m-1)}}_{A}+\underbrace{(y_{i}-\mu_{i})\left(\frac{\partial\eta_{i}}{\partial\mu_{i}}\right)}_{B}; con \hspace{0.05cm} k=1,2.\hspace{1cm}(3) \]

Luego, A corresponde a: \[A=\sum_{k}x_{lk}b_{k}^{(m-1)}=b_{1}+b_{2}x_{i}\]

Si $\frac{\partial\mu_{i}}{\partial\eta_{i}}=e^{\eta_{i}}$, luego $\frac{\partial\eta_{i}}{\partial\mu_{i}}=\frac{1}{e^{\eta_{i}}}$, por lo tanto $\frac{\partial\eta_{i}}{\partial\mu_{i}}=\frac{1}{e^{b_{1}+b_{2}x_{i}}}$, reemplazando(1) y $\frac{\partial\eta_{i}}{\partial\mu_{i}}$ tenemos que: \[B=(y_{i}-\mu_{i})\left(\frac{\partial\eta_{i}}{\partial\mu_{i}}\right)\] \[B=(y_{i}-e^{b_{1}+b_{2}x_{i}})\frac{1}{e^{b_{1}+b_{2}x_{i}}}\], distribuyendo obtenemos que: \[B=(y_{i}e^{-(b_{1}+b_{2}x_{i})}-1)\] reemplazando A y B en (3), se tiene una expresión para determinar $z_{i}$ \[z_{i}=b_{1}+b_{2}x_{i}+(y_{i}e^{-(b_{1}+b_{2}x_{i})}-1)\]

Para los datos de la Tabla 4.3, estime los parámetros del modelo adaptando la macro Rstudio.

Para realizar la estimación de los parámetros del modelo, empleamos un algoritmo que permite a través del método de aproximación del Scoring hallar la estimación de los parámetros.

dat4.1 <- read.csv("/cloud/project/dat4.1.csv", sep=";")
dt=as.matrix(dat4.1);colnames(dt)=NULL
Y=dt[,1]#(Variable respuesta)
uno=rep(1,each=1,times=14)
X=as.matrix(data.frame(uno,dt[,2]))#(Matriz de datos)
colnames(X)=NULL
f <- function(b){
  M4=(X%*%b)#c1 
  M=as.vector(1/exp(M4))
  C2=as.vector(exp(M4))#varia
  W=as.matrix(diag(C2))#new M4 
  Xt=t(X)#M5
  XtW=Xt%*%W#M6 
  M9=Y*M
  z=M4+M9-uno
  XtWz=XtW%*%z#M7
  XtWX=XtW%*%X#M8 # Matriz de fisher
  iXtWX=solve(XtWX)#M9
  b=(iXtWX%*%XtWz)
  print(b)
}
b=as.matrix(c(-1.5,1.5))
delta <- 1e-8 # tolerancia
iterM <- 15 # numero maximo de iteraciones
tolera <- 1 # inicializar tolera
m <- 0 # inicializar itera
lista<- b # inicializar historial de iteraciones
while((tolera>delta)&(m<iterM)){
  xold <- b# b viejo
  b <-f(b) #b nuevo
  tolera <- abs( b - xold )
  lista<- data.frame(lista,b)
  m<- m + 1
}

##           [,1]
## [1,] -3.160700
## [2,]  3.139537
##           [,1]
## [1,] -2.803736
## [2,]  2.718729
##           [,1]
## [1,] -2.166138
## [2,]  2.306425
##           [,1]
## [1,] -1.953811
## [2,]  2.180810
##           [,1]
## [1,] -1.943680
## [2,]  2.174564
##           [,1]
## [1,] -1.943654
## [2,]  2.174548
##           [,1]
## [1,] -1.943654
## [2,]  2.174548

length(lista)

## [1] 8

Notamos que el algoritmo por el método de Fisher Scoring converge en 7 iteraciones y la estimación de los parámetros por máxima verosimilitud corresponden a $\hat{b_{1}}=-1.943654$ y $\hat{b_{2}}=2.174548$.

EJERCICIO 4.3 (DOBSON, 1990): Los datos de la Tabla 4.4 son los tiempos de muerte, en $Y$, semanas desde el diagnóstico y $log_{10}$ (recuento inicial de glóbulos blancos), $X$, para diecisiete pacientes que padecen leucemia. (Este es el Ejemplo U de Cox y Snell, 1981).

Tabla 4.4. Tiempo de sobrevivencia $Y$, en semanas y $log_{10}$ (recuento inicial de glóbulos blancos) $X$, para diecisiete pacientes que padecen leucemia
Y	65.00	156.00	100.00	134.00	16.00	108.00	121	4.00	39.00	143.00	56.00	26.00	22.00	1	1	5.00	65
X	3.36	2.88	3.63	3.41	3.78	4.02	4	4.23	3.73	3.85	3.97	4.51	4.54	5	5	4.72	5

Grafique $y_{i}$ contra $x_{i}$. ¿Los datos muestran alguna tendencia?

No se muestra ninguna tendencia o relación gráfica entre la variable de respuesta y la predictora

Una posible especificación para E(Y) es $E(Y_{i})=exp(\beta_{1}+\beta_{2}x_{i})$ que asegurará que E(Y) no sea negativo para todos los valores de los parámetros y todos los valores de x. ¿Qué función de enlace es apropiada en este caso?

Las funciones exponenciales, son las funciones que tienen la variable independiente x en el exponente, es decir, son de la forma:

\[f(x)=(ab)^{x}\]

donde la base correspone a un número real mayor que cero y diferente de uno, encontramos además el caso particular en el que la base corresponde a la constante euler, independientemente de cual sea la escogencia de la base siempre y cuando se mantengan las condiciones se cumplirá que su dominio es $\mathbb{R}$ y el recorrido es de $(0,\infty)$. En este caso, esto nos asegurará que E(Y) no sea negativa para todos los valores de los parámetros y todos los valores de x. El gráfico no sugiere una tendencia lineal, por lo que se recomienda utilizar un modelo lineal generalizado, y enlace “log” para el experimento se estudia tiempo de supervivecia antes del deceso de un paciente con leucemia por lo que es recomendable utilizar un distribución exponencial debido a que nuestra variable respuesta se obtine a partir de un intervalo entre el diagnóstico y el deceso.

La distribución exponencial se usa a menudo para describir los tiempos de supervivencia. Demuestre que este es un caso especial de la distribución Gamma (véanse los Ejercicios 3.1 y 3.2 (b)).

Veamos que la distribución exponencial es un caso especial de la distribución Gamma cuando el parámetro $\alpha=1$.

Esta prueba puede ser equivalente a mostrar que la función de momentos de la distribución Gamma es igual a la distribución exponencial cuando $\alpha=1$

Sea $y$ una variable aleatoria con distribución exponencial tal que la función de densidad de probabilidad viene dada por:

\[ f(Y)=\lambda\ e^{(-y\lambda)}, \hspace{1cm} y>0,\lambda>0\] \[\begin{align*} M(t)=E(e^{(Ty)})=&\int_{0}^{\infty}e^{(ty)}\lambda*e^{(-y\lambda)}dy\\ &= \int_{0}^{\infty}\lambda*e^{(t-\lambda)y}dy \end{align*}\] Resolviendo la integral tenemos que: \[\begin{align*} &=\dfrac{\lambda}{(t-\lambda)}(e^{(t-y)})|_{0}^{\infty}\\ &=\dfrac{\lambda}{(t-\lambda)}[\lim\limits_{y\rightarrow infty}e^{(t-\lambda)y}-e^{(t-\lambda)0}]\\ \end{align*}\]

Evaluando la integral en el intervalo definido \[\begin{align*} &=\dfrac{\lambda}{(t-\lambda)}[0-1]=\dfrac{\lambda}{(t-\lambda)}\hspace{1cm} t-\lambda<0,t<\lambda \end{align*}\]

así \[M(T)=\dfrac{\lambda}{(t-\lambda)};\hspace{1cm}t-\lambda<0,t<\lambda\]

Sea X una variable aleatoria con distribución Gamma sobre el intervalo $[0,\infty)$ con parámetros $\alpha$ y $\lambda$, se define la función de densidad de probabilidad como sigue:

\[f(X)=\dfrac{\lambda^{(\alpha)}}{(\Gamma(\alpha))}x^{(\alpha-1)}e^{(-\lambda\alpha)},\hspace{1cm} y>0,\lambda>0\]

Por definición \[\Gamma(\alpha)=\int_{0}^{\infty}t^{(\alpha-1)}e^{(-\lambda)t}dt\]

y para valores enteros de $\alpha$ tenemos que $\Gamma(\alpha)=(\alpha-1)!$

La función de momentos se define como:

\[\begin{align*} M(t)=E(e^{(tx)})&=\int_{0}^{\infty}e^{(tx)}\dfrac{\lambda^{(\alpha)}x^{(\alpha-1)}e^{(-\lambda)x}}{\Gamma(\alpha)}dx\\ &=\int_{0}^{\infty}\dfrac{\lambda^{(\alpha)}}{\Gamma(\alpha)}x^{(\alpha-1)}e^{(t-\lambda)}x dx \end{align*}\]

Hagamos la siguiente sustitución $y=(\lambda-t)x$ de modo que la integral se transforma en

\[M(t)=\int_{0}^{\infty}\dfrac{\lambda^{(\alpha)}}{\Gamma(\alpha)}\dfrac{y^{(\alpha-1)}}{(\lambda-t)^{(\alpha-1)}}e^{(-y)}\dfrac{1}{(\lambda-t)}dy\] \[=\dfrac{\lambda^{(\alpha)}}{(\Gamma(\alpha))(\lambda-t)^{(\alpha)}}\int_{0}^{\infty}y^{(\alpha-1)}dy\] \[=\dfrac{\lambda^{(\alpha)}}{(\Gamma(\alpha))(\lambda-t)^{(\alpha)}}\Gamma(\alpha)\] \[M(t)=\dfrac{\lambda^{(\alpha)}}{(\lambda-t)^{(\alpha)}}\]

En donde para el caso específico $\alpha=2$ se cumple que es equivalente a la función de momentos de la distribución exponencial.

Utilice R para ajustar el modelo sugerido en las partes (b) y (c) Trace el modelo ajustado en el gráfico obtenido en el inciso (a). ¿Considera que el modelo es una descripción adecuada de los datos?

Para realizar la estimación de los parámetros del modelo, adoptamos el macro Rstudio, empleando la función glm del paquete {stats},

## 
## Call:
## glm(formula = datos1$Y ~ datos1$X, family = Gamma(link = "log"), 
##     data = datos1)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -1.9922  -1.2102  -0.2242   0.2102   1.5646  
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)   8.4775     1.6548   5.123 3.01e-07 ***
## datos1$X     -1.1093     0.3997  -2.776  0.00551 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for Gamma family taken to be 1)
## 
##     Null deviance: 26.282  on 16  degrees of freedom
## Residual deviance: 19.457  on 15  degrees of freedom
## AIC: 173.97
## 
## Number of Fisher Scoring iterations: 8

Del análisis se concluye que el modelo es una representación adecuada de los datos, dado que a un nivel de significancia del 5% los parámetros $X_{i}$ e intercepto son significativos con p-valor 0.00551 y 0.00551, respectivamente.

EJERCICIO 4.4 (DOBSON, 1990)

Se puede obtener una derivación alternativa de la ecuación de Newton-Raphson (4.5) aproximando la función de verosimilitud logarítmica $l(\mathbf{\beta},\mathbf{y})$ mediante una expansión en serie de Taylor sobre $\mathbf{\beta}=\mathbf{\beta^*}$.La ecuación utilizada es: \[l(\mathbf{\beta},\mathbf{y})=l(\mathbf{\beta}^*,\mathbf{y})+(\mathbf{\beta}-\beta^*)^TU+(\mathbf{\beta}-\mathbf{\beta}^*)^TH(\mathbf{\beta}-\mathbf{\beta}^*)\hspace{0.5cm}\] donde $U$, el vector $p\times 1$ con elementos $U_{j}=\frac{\partial l}{\partial\beta_j}$, y $H$, la matriz $p\times p$ con elementos $\frac{\partial^2l}{\partial\beta_j\partial\beta_k}$ evaluados en $\beta=\beta^*$.

Escriba la versión de parámetro único de esta aproximación y utilícela para obtener una expresión para el estimador de máxima verosimilitud para $\beta$. Si $\beta^*$ se considera la (m - l) aproximación y $\beta$ como la m-ésima aproximación muestre que la ecuación corresponde a la versión de un solo parámetro de (4.5).
Nota: La aproximación a la serie de Taylor para una función $f(x)$ para una variable $x$ alrededor del valor $t$ está dada a partir de la siguiente expresión: \[f(x)=\sum_{n=0}^\infty \begin{bmatrix}\frac{d^nf}{dx^n}\end{bmatrix}_{x=t}(x-t)\] Una aproximación de la función de log-verosimilitud $l(\mathbf{\beta},\mathbf{y})$ para una variable aleatoria $\mathbf{y}$ con $\mathbf{y}=(y_{1},\cdots,y_{n})^T$ usando la expansión de Taylor alrededor de $\mathbf{\beta}=\mathbf{\beta^*}$ es la siguiente. Los primeros tres términos son: \[l(\mathbf{\beta},\mathbf{y})=l(\mathbf{\beta}^*,\mathbf{y})+(\mathbf{\beta}-\beta^*)^TU+(\mathbf{\beta}-\mathbf{\beta}^*)^TH(\mathbf{\beta}-\mathbf{\beta}^*)\hspace{0.5cm}(1)\] Donde $U$ es un vector de $p\times 1$ cuyos elementos son de la forma $U_{j}=\frac{\partial l}{\partial\beta_j}$ y $H$ es una matriz de $p\times p$ con elementos de la forma $\frac{\partial^2l}{\partial\beta_j\partial\beta_k}$ evaluados en $\beta=\beta^*$

Si $\beta$ es un parámetro único, entonces $U=\frac{d l(\beta)}{d\beta}\left.\begin{matrix}\end{matrix}\right|_{\beta=\beta^*}$ y $H=\frac{d^2 l(\beta)}{d\beta ^2}\left.\begin{matrix}\end{matrix}\right|_{\beta=\beta^*}$

Por lo tanto, la expresión $(1)$ se reescribe como: \[l(\beta,\mathbf{y})=l(\beta^*,\mathbf{y})+(\beta-\beta^*)\begin{bmatrix}\frac{dl(\beta)}{d\beta}\end{bmatrix}_{\beta=\beta^*}+\frac{1}{2}(\beta-\beta^*)^2\begin{bmatrix}\frac{d^2 l(\beta)}{d\beta}\end{bmatrix}_{\beta=\beta^*}\hspace{0.5cm}(2)\]

Como nos interesa obtener una expresión para el estimador por máxima verosimilitud de $\beta$, derivamos la expresion $(2)$ con respecto a $\beta$, llegando así a la siguiente expresión. \[\frac{dl (\beta,\mathbf{y})}{d\beta}=\frac{dl(\beta^*,\mathbf{y})}{d\beta}+ \left [ \frac{dl(\beta)}{d\beta} \right ]_{\beta=\beta^*}+(\beta-\beta^*)\left [ \frac{d ^2 l(\beta)}{d \beta^2} \right ]_{\beta=\beta^*}\] Donde $\frac{dl(\beta^*,y)}{d\beta}=0$ por ser derivada de un término constante, teniendo así lo siguiente: \[\begin{align*} 0&=\begin{bmatrix}\frac{d l(\beta)}{d\beta}\end{bmatrix}_{\beta=\beta^*}+(\beta-\beta^*)\begin{bmatrix}\frac{d^2l(\beta)}{d\beta^2}\end{bmatrix}_{\beta=\beta^*}\\ - \left [\frac{dl(\beta)}{d \beta}\right ]_{\beta=\beta^*}\left ( \left [\frac{d^2l(\beta)}{d\beta^2}\right ]_{\beta=\beta^*}\right )^{-1}&=(\beta-\beta^*)\hspace{0.5cm}(3) \end{align*}\] Reemplazando $b$, donde $b$ es el estimador por máxima verosimilitud de $\beta$. Así, $(3)$ se reescribe como: \[-\left [ \frac{d l(\beta)}{d \beta} \right ]_{\beta=b^*}\left [ \frac{d ^2 l(\beta)}{d \beta^2} \right ]_{\beta=b^{*}}^{-1}=b-b^*\] $\Rightarrow$ \[b=b^*-\left [\frac{d l(\beta)}{d l\beta} \right ]_{\beta=b^*}\left [ \frac{d^2l(\beta)}{d \beta^2}\right ]_{\beta=b^*}^{-1}\]

Demuestre el resultado correspondiente para el caso general.

Ahora si suponemos que $\beta$ no es un parámetro único, sino un vector de parámetros; y procedemos a hacer una aproximación de la función log-verosimilitud $l(\mathbf{\beta},\mathbf{y})$ usando la expansión de Taylor alrededor de $\mathbf{\beta}=\mathbf{\beta^*}$, los primeros tres términos de la serie quedan definidos como en $(1)$.

Donde, $U$ es un vector $p\times 1$ de la forma $U_j=\frac{\partial l}{\partial\beta_j}$ que contiene todas las derivadas parciales de $l(\mathbf{\beta})$ con respecto a los parámetros $\mathbf{\beta}$, este vector recibe el nombre de score. Además, $H$ es una matriz de $p\times p$ que contiene las derivadas parciales de segundo orden evaluadas en $\mathbf{\beta}=\mathbf{\beta*}$, también llamada (Matriz Hessiana).

Procedemos a obtener una expresión que permita estimar los parámetros $\mathbf{\beta}$ por máxima verosimilitud, realizamos las derivadas parciales de $(1)$ con respecto a $\mathbf{\beta}$, las cuales tienen la estructura de un vector de derivadas parciales de orden $p\times 1$,la derivada queda definida como sigue. \[\frac{\partial l(\mathbf{\beta},\mathbf{y})}{\partial \mathbf{\beta}}=\frac{\partial l(\mathbf{\beta^*},\mathbf{y})}{\partial\mathbf{\beta}}+U+ \frac{\partial}{\partial\mathbf{\beta}}\left[\frac{(\mathbf{\beta}-\mathbf{\beta}^*)^TH(\mathbf{\beta}-\mathbf{\beta}^*)}{2} \right ]\hspace{0.5cm}(4)\] Desde luego \[\frac{\partial l(\mathbf{\beta^*},\mathbf{y})}{\partial\mathbf{\beta}}=0\hspace{0.5cm}(5)\] dado que $l(\mathbf{\beta}^*,\mathbf{y})$ es un vector numérico al estar evaluado en $\mathbf{\beta}^*.$

: Si $A$ es una matriz simétrica de orden $n$ y $\overrightarrow{x}\in \mathbb{R}^n$ y $f(x)=\overrightarrow{x}^TA\overrightarrow{x}$, entonces $\frac{\partial f(\overrightarrow{x})}{\partial \overrightarrow{x}}=2A\overrightarrow{x}$

Usando este argumento, como $H$ es una matriz simétrica, tenemos que si $g(\mathbf{\beta}-\mathbf{\beta}^*)=(\mathbf{\beta}-\mathbf{\beta}^*)^TH(\mathbf{\beta}-\mathbf{\beta}^*)$ con $\mathbf{\beta},\mathbf{\beta}^*\in \mathbb{R}^{p\times1}$, entonces \[{g}'(\mathbf{\beta}-\mathbf{\beta}^*)=2H(\mathbf{\beta}-\mathbf{\beta}^*)\hspace{0.5cm}(6)\]

Por definición, para obtener le estimación de $\mathbf{\beta}$ que es $\widehat{\mathbf{\beta}}$ hacemos $\frac{\partial l(\mathbf{\beta},\mathbf{y})}{\mathbf{\beta}}=0$, luego (4) queda como se muestra a continuación, después de reemplazar (5) y (6):

\[0=U+H(\mathbf{\beta}-\mathbf{\beta}*)\] \[-UU=H(\mathbf{\beta}-\mathbf{\beta}^*)\] \[-H^{-1}U=H^{-1}H(\mathbf{\beta}-\mathbf{\beta}^*)\] \[-H^{-1}U=(\mathbf{\beta}-\mathbf{\beta}^*)\] \[\widehat{\mathbf{\beta}}=\widehat{\mathbf{\beta}}^*-H^{-1}U\]

SECCIÓN 4. EJERCICIOS CAPÍTULO 4 (MYERS, MONTGOMERY, VINING & ROBINSON, 2010)

EJERCICIO 4.5 (MYERS, MONTGOMERY, VINING & ROBINSON, 2010)

Anand (1997) analiza un experimento para mejorar el rendimiento de gel de sílice en una planta química. Los factores fueron (A) densidad de silicato de sodio, (B) pH, (C) tiempo de fraguado y (D) temperatura de secado. La respuesta fue el porcentaje de lotes que se consideró material de segundo grado de 1000 lotes. Los datos son los siguientes.

Cuadro 4.5
A	B	C	D	P
1.125	3.5	40	150	14.0
1.125	3.5	50	120	13.5
1.125	3.5	30	135	18.3
1.125	3.0	40	120	17.4
1.125	3.0	50	135	16.3
1.125	3.0	30	150	13.9
1.125	4.0	40	135	14.1
1.125	4.0	50	150	12.0
1.125	4.0	30	120	12.0
1.121	3.5	40	120	17.0
1.121	3.5	50	135	10.5
1.121	3.5	30	150	15.3
1.121	3.0	40	135	21.0
1.121	3.0	50	150	20.5
1.121	3.0	30	120	19.8
1.121	4.0	40	150	12.3
1.121	4.0	50	120	11.3
1.121	4.0	30	135	10.5
1.123	3.5	40	135	15.2
1.123	3.5	50	150	12.3
1.123	3.5	30	120	12.0
1.123	3.0	40	150	11.1
1.123	3.0	50	120	12.8
1.123	3.0	30	135	13.4
1.123	4.0	40	120	10.3
1.123	4.0	50	135	10.4
1.123	4.0	30	150	11.7

Utilice el enlace logit y la regresión logística para analizar estos datos.

Notamos que la variable respuesta $p_{i}$ con $i=1,\cdots,27.$ es el porcentaje de lotes en el que se consideró material de segundo grado de 1000 lotes, entonces construimos una variable $y_{i}$ equivalente a la expresión decimal de estos porcentajes y además se toma la función de enlace $g(\mu_{i})=log(\frac{y}{1-y})$.

Como estamos frente a un modelo con varias variables regresoras, determinamos la matriz de correlación, para tener una idea de cuáles variables probablemente no van a tener algún efecto sobre la variable de respuesta.

##      A    B     C     D     y
## A  1.0  0.0  0.00  0.00 -0.10
## B  0.0  1.0  0.00  0.00 -0.60
## C  0.0  0.0  1.00  0.00 -0.11
## D  0.0  0.0  0.00  1.00 -0.04
## y -0.1 -0.6 -0.11 -0.04  1.00

De la matriz de correlación, notamos que la variable pH es la única que tiene efecto sobre el porcentaje de lotes, debido a que tiene un alto grado de correlación, siendo exactos de -0.6.

Para realizar la estimación de los parámetros del modelo logístico logit, adoptamos el macro Rstudio, empleando la función glm del paquete {stats}, en donde se especifica family=“binomial” y la función de enlace link=“logit”.

M_logistico <- glm(y ~A+B+C+D,family="binomial"("logit"),data=bdatos)
summary(M_logistico)

## 
## Call:
## glm(formula = y ~ A + B + C + D, family = binomial("logit"), 
##     data = bdatos)
## 
## Deviance Residuals: 
##       Min         1Q     Median         3Q        Max  
## -0.145350  -0.055786   0.005221   0.034967   0.120635  
## 
## Coefficients:
##               Estimate Std. Error z value Pr(>|z|)
## (Intercept)  1.716e+01  3.818e+02   0.045    0.964
## A           -1.547e+01  3.399e+02  -0.046    0.964
## B           -3.850e-01  1.367e+00  -0.282    0.778
## C           -3.372e-03  6.798e-02  -0.050    0.960
## D           -9.227e-04  4.531e-02  -0.020    0.984
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 0.21144  on 26  degrees of freedom
## Residual deviance: 0.12661  on 22  degrees of freedom
## AIC: 18.179
## 
## Number of Fisher Scoring iterations: 5

De manera general, vemos que el modelo logístico logit no es el más apropiado, lo que permite concluir que a un nivel de significancia del 5% los parámetros (A, B, C y D) no son significativos con p-valores de 0.964, 0.778, 0.960 y 0.984, respectivamente.

Repita este análisis utilizando el enlace probit.

Para realizar la estimación de los parámetros del modelo logístico probit empleamos nuevamente la función glm del paquete {stats} en el macro Rstudio, en donde se especifica family=“binomial” y la función de enlace link=“probit”.

M_probit <- glm(y ~A+B+C+D,family=binomial(link="probit"),data=bdatos)
summary(M_probit)

## 
## Call:
## glm(formula = y ~ A + B + C + D, family = binomial(link = "probit"), 
##     data = bdatos)
## 
## Deviance Residuals: 
##       Min         1Q     Median         3Q        Max  
## -0.145570  -0.055366   0.004363   0.034427   0.118986  
## 
## Coefficients:
##               Estimate Std. Error z value Pr(>|z|)
## (Intercept)  8.485e+00  2.067e+02   0.041    0.967
## A           -7.749e+00  1.840e+02  -0.042    0.966
## B           -2.082e-01  7.379e-01  -0.282    0.778
## C           -1.850e-03  3.680e-02  -0.050    0.960
## D           -4.672e-04  2.453e-02  -0.019    0.985
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 0.21144  on 26  degrees of freedom
## Residual deviance: 0.12674  on 22  degrees of freedom
## AIC: 18.179
## 
## Number of Fisher Scoring iterations: 5

De manera general, vemos que el modelo logístico probit tampoco es el más apropiado, se concluye que a un nivel de significancia del 5% los parámetros (A, B, C y D) no son significativos con p-valores de 0.966, 0.778, 0.960 y 0.985, respectivamente.

Comente los dos análisis.

##        Intercept          A          B            C             D Residual dev.
## Logit  17.160265 -15.472185 -0.3850023 -0.003371718 -0.0009227031     0.1266080
## Probit  8.485077  -7.749073 -0.2082121 -0.001850028 -0.0004671764     0.1267426
##             AIC
## Logit  18.17922
## Probit 18.17925

De los dos análisis, inferimos que los modelos de regresión “logit” y "probit no son significativos, ambos proporcionan un ajuste similar, no hay diferencias entre los valores de los criterios (Deviance y AIC), se recomienda emplear otro tipo de modelo que se ajuste mejor a los datos.

La distribución límite para el binomio es el de Poisson. Repita la parte (a) usando regresión de Poisson. Discuta las diferencias.

Para realizar la estimación de los parámetros del modelo, adoptamos nuevamente el macro Rstudio, empleando la función glm del paquete {stats}, en donde se especifica la family=poisson y su función de enlace correspondiente. Además, la variable respuesta a considerar es P.

## 
## Call:
## glm(formula = P ~ A + B + C + D, family = poisson(link = "log"), 
##     data = bdatos)
## 
## Deviance Residuals: 
##      Min        1Q    Median        3Q       Max  
## -1.34560  -0.52425   0.05306   0.32427   1.11607  
## 
## Coefficients:
##               Estimate Std. Error z value Pr(>|z|)   
## (Intercept)  18.896383  35.341322   0.535  0.59287   
## A           -13.254240  31.465183  -0.421  0.67358   
## B            -0.330869   0.126700  -2.611  0.00902 **
## C            -0.002889   0.006293  -0.459  0.64623   
## D            -0.000789   0.004195  -0.188  0.85081   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 18.037  on 26  degrees of freedom
## Residual deviance: 10.746  on 22  degrees of freedom
## AIC: Inf
## 
## Number of Fisher Scoring iterations: 4

Se concluye que a un nivel de significancia del 5% los parámetros (A, C y D) no son significativos con p-valores de 0.67358, 0.64623 y 0.85081, respectivamente, al igual que el intercepto. Notamos además, que la variable B que corresponde al pH es la única que tiene efecto sobre la variable respuesta Porcentaje, entonces construimos nuevamente un modelo Poisson pero empleando solamente pH como variable regresora.

Modelo reducido

## 
## Call:
## glm(formula = P ~ B, family = poisson(link = "log"), data = bdatos)
## 
## Deviance Residuals: 
##      Min        1Q    Median        3Q       Max  
## -1.39274  -0.44092  -0.02502   0.35477   1.12318  
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)   3.7904     0.4395   8.624  < 2e-16 ***
## B            -0.3309     0.1267  -2.611  0.00901 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 18.037  on 26  degrees of freedom
## Residual deviance: 11.170  on 25  degrees of freedom
## AIC: Inf
## 
## Number of Fisher Scoring iterations: 4

## Waiting for profiling to be done...

##                 2.5 %      97.5 %
## (Intercept)  2.926283  4.65020010
## B           -0.580209 -0.08323441

Del análisis se concluye que a un nivel de significancia del 5% el parámetro correspondiente a la variable (pH) y el término independiente de la regresión son significativos con p-valores de 0.00901 y 2e-16, respectivamente, hecho que se confirma al observar los intervalos de confianza, ya que ninguno contiene el cero. A continuación, calculamos la significancia del modelo logístico completo (p-valor del modelo), a través del test Likelihood ratio.

Evaluación del Modelo

##   Diferencia.de.residuos Grados.de.libertad p.value
## 1                 6.8666                  1  0.0088

En este caso, el estadístico con distribución chi-cuadrado corresponde a 6.8666, con un grado de libertad; y un p-valor=0.0088, que comparado con el nivel de significancia del 0.05 afirma que el modelo Poisson con una variable regresora es significativo.

EJERCICIO 4.6 (MYERS, MONTGOMERY, VINING & ROBINSON, 2010)

Nelson (1982, págs. 407-409) analiza un experimento para determinar la relación del tiempo de uso con el número de fisuras que se desarrollan en las ruedas de furbina. Los datos son los siguientes.

Cuadro 4.6
Hours	Turbines	Fissures
400	39	0
1000	53	4
1400	33	2
1800	73	7
2200	30	5
2600	39	9
3000	42	9
3400	13	6
3800	34	22
4200	40	21
4600	36	21

Utilice el enlace logit y la regresión logística para analizar estos datos. Para realizar la estimación de los parámetros del modelo, adoptamos el macro Rstudio, empleando la función glm del paquete {stats}, en donde se especifica la familia y la función de enlace correspondiente.

Inicialmente notamos que la variable respuesta $y_{i}$ es una variable conteo. Si deseamos realizar un análisis de regresión logística, esta variable debe tener distribución binomial, por lo que es necesario una transformación, para que represente la proporción de fisuras por número de turbinas; y además se toma la función de enlace $g(\mu_{i})=log(\frac{p}{1-p})$.

tur.m1 <- glm( Fissures/Turbines ~ Hours, family=binomial,        weights=Turbines, data=turbines)
summary(tur.m1)

## 
## Call:
## glm(formula = Fissures/Turbines ~ Hours, family = binomial, data = turbines, 
##     weights = Turbines)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -1.5055  -0.7647  -0.3036   0.4901   2.0943  
## 
## Coefficients:
##               Estimate Std. Error z value Pr(>|z|)    
## (Intercept) -3.9235966  0.3779589 -10.381   <2e-16 ***
## Hours        0.0009992  0.0001142   8.754   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 112.670  on 10  degrees of freedom
## Residual deviance:  10.331  on  9  degrees of freedom
## AIC: 49.808
## 
## Number of Fisher Scoring iterations: 4

#Intervalos de confianza
confint(object = tur.m1, level = 0.95 )

## Waiting for profiling to be done...

##                    2.5 %       97.5 %
## (Intercept) -4.704479229 -3.219091395
## Hours        0.000783381  0.001231879

Del análisis se concluye que a un nivel de significancia del 5% los parámetros Horas e intercepto son significativos con un p-valor de 2e-16, hecho que que se confirma al observar los intervarlos de confianza, ya que ninguno contiene el cero. A continuación, calculamos la significancia del modelo logístico completo (p-valor del modelo), a través del test Likelihood ratio.

Evaluación del Modelo

##   Diferencia.de.residuos Grados.de.libertad p.value
## 1               102.3386                  1       0

Para ello, construirmos un estadístico con la diferencia de residuos entre el modelo 1 (modelo con predictores) y el modelo 2 (modelo sin predictores).

El estadístico anterior tiene una distribución chi-cuadrado con grados de libertad equivalentes a $m_{1}-m_{2}$, es decir, la diferencia de grados de libertad entre el modelo 1 y el modelo 2. Dado que el p-valor es 0, comparado con el nivel de significancia afirmamos que el modelo 1 es altamente significativo.

Repita este análisis utilizando el enlace probit.

En este caso, empleamos en la función de enlace probit $\eta_{i}=\phi^{-1}(\mu_{i})=probit(\mu_{i})$, donde $\Phi(.)$ es la función de densidad acumulada para la distribución normal. Esta función de enlace se especifica en R como link=“probit”.

## 
## Call:
## glm(formula = Fissures/Turbines ~ Hours, family = binomial(link = "probit"), 
##     data = turbines, weights = Turbines)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -1.2696  -0.6399  -0.2447   0.3859   2.0894  
## 
## Coefficients:
##               Estimate Std. Error z value Pr(>|z|)    
## (Intercept) -2.2758075  0.1974187 -11.528   <2e-16 ***
## Hours        0.0005783  0.0000626   9.239   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 112.6700  on 10  degrees of freedom
## Residual deviance:   9.8148  on  9  degrees of freedom
## AIC: 49.292
## 
## Number of Fisher Scoring iterations: 4

##   (Intercept)         Hours 
## -2.2758074623  0.0005783211

## Waiting for profiling to be done...

##                     2.5 %        97.5 %
## (Intercept) -2.6760191509 -1.9006694833
## Hours        0.0004583096  0.0007036015

Del análisis obtenido luego de emplear la función de enlace “probit”, se concluye que a un nivel de significancia del 5% los parámetros Horas e intercepto son significativos con un p-valor de 2e-16, hecho que se confirma al observar los intervarlos de confianza ya que ninguno contine el cero. A continuación calculamos la significancia del modelo logístico completo con probit (p-valor del modelo), a través del test Likelihood ratio.

Evaluación del Modelo

De manera similar a lo hecho con el modelo logit, procedemos con el modelo probit. En este caso, el estadístico con distribución chi-cuadrado corresponde a 102.8552, con un grado de libertad; y un p-valor de 0, que comparado con el nivel de significancia del 0.05 afirma que el modelo completo es significativo.

##   Diferencia.de.residuos Grados.de.libertad p.value
## 1               102.8552                  1       0

Comente los dos análisis.

##        Intercept        Hours Residual dev.      AIC
## Logit  -3.923597 0.0009992372     10.331466 49.29156
## Probit -2.275807 0.0005783211      9.814837 49.80819

Para comparar el modelo logístico completo, resultante de usar la función de enlace “logit” y “probit”, tenemos en cuenta los criterios Akaike y Deviance definidos con anterioridad.

En este caso, inferimos que los modelos de regresión “logit” y “probit” proporcionan un ajuste similar, debido a que el modelo logit tiene mayor deviance pero menor AIC que el modelo probit, y el modelo probit tiene menor deviance pero mayor AIC que el modelo logit. Cabe resaltar que la diferencia entre los valores de los criterios no es muy significativa, lo que genera dificultad para decidir cuál modelo es el más adecuado.

El análisis gráfico confirma lo dicho anteriormente, ambos modelos proporcionan un ajuste de los datos similar.

La distribución límite para el binomio es el de Poisson. Repita la parte (a) usando regresión de Poisson. Discuta las diferencias.

Como investigadores, obtenemos dos modelos, un modelo 1 en donde consideramos el número de horas como variable predictora, y un modelo dos en donde consideramos horas y turbinas como variables influyentes.

Modelo 1

## 
## Call:
## glm(formula = Fissures ~ Hours, family = poisson(link = "log"), 
##     data = turbines)
## 
## Deviance Residuals: 
##      Min        1Q    Median        3Q       Max  
## -1.92935  -0.82584  -0.08907   0.68502   1.63581  
## 
## Coefficients:
##              Estimate Std. Error z value Pr(>|z|)    
## (Intercept) 3.742e-01  3.308e-01   1.131    0.258    
## Hours       6.175e-04  9.069e-05   6.809 9.81e-12 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 69.252  on 10  degrees of freedom
## Residual deviance: 13.745  on  9  degrees of freedom
## AIC: 57.394
## 
## Number of Fisher Scoring iterations: 5

##  (Intercept)        Hours 
## 0.3742167410 0.0006175167

## [1] 13.7447

## Analysis of Deviance Table
## 
## Model: poisson, link: log
## 
## Response: Fissures
## 
## Terms added sequentially (first to last)
## 
## 
##       Df Deviance Resid. Df Resid. Dev  Pr(>Chi)    
## NULL                     10     69.252              
## Hours  1   55.507         9     13.745 9.313e-14 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Se concluye que a un nivel de significancia del 5%, el parámetro Hora es significativo con un p-valor de 9.81e-12 y un deviance de 13.7447.

Evaluación del modelo

##   Diferencia.de.residuos Grados.de.libertad p.value
## 1                 55.507                  1       0

En este caso, el estadístico con distribución chi-cuadrado corresponde a 55.507, con un grado de libertad; y un p-valor=0, que comparado con el nivel de significancia del 0.05 afirma que el modelo completo es significativo.

El análisis gráfico confirma que con una sola variable predictiva, el modelo proporciona un buen ajuste de los datos.

Modelo 2

tr.poiss=glm(Fissures~Hours+Turbines, family=poisson, data=turbines)
summary(tr.poiss)

## 
## Call:
## glm(formula = Fissures ~ Hours + Turbines, family = poisson, 
##     data = turbines)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -1.7711  -0.5171  -0.2580   0.3614   1.8168  
## 
## Coefficients:
##               Estimate Std. Error z value Pr(>|z|)    
## (Intercept) -5.604e-01  5.665e-01  -0.989    0.322    
## Hours        6.828e-04  9.602e-05   7.111 1.15e-12 ***
## Turbines     1.891e-02  8.660e-03   2.183    0.029 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 69.2517  on 10  degrees of freedom
## Residual deviance:  8.9842  on  8  degrees of freedom
## AIC: 54.634
## 
## Number of Fisher Scoring iterations: 4

coef(tr.poiss)

##   (Intercept)         Hours      Turbines 
## -0.5604432333  0.0006827818  0.0189076391

deviance(tr.poiss)

## [1] 8.984185

anova(tr.poiss,test="Chisq")

## Analysis of Deviance Table
## 
## Model: poisson, link: log
## 
## Response: Fissures
## 
## Terms added sequentially (first to last)
## 
## 
##          Df Deviance Resid. Df Resid. Dev  Pr(>Chi)    
## NULL                        10     69.252              
## Hours     1   55.507         9     13.745 9.313e-14 ***
## Turbines  1    4.761         8      8.984   0.02912 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Se concluye que a un nivel de significancia del 5%, los parámetros Hora y Turbinas son significativos con un p-valor de 9.313e-14 y 0.02912, respectivamente, y un deviance de 8.984185.

Para proceder con la comparación de los tres modelos (Modelo logístico “logit”, Modelo logístico “probit” y Modelo Poisson), es necesario elegir cuál modelo Poisson nos brinda un mejor ajuste; en este caso elegimos el modelo 2 que presenta un AIC y Deviance menor.

Comparación de los tres modelos (Análisis)

##          Intercept        Hours   Turbines Residual dev.      AIC
## Logit   -3.9235966 0.0009992372 0.00000000     10.331466 49.29156
## Probit  -2.2758075 0.0005783211 0.00000000      9.814837 49.80819
## Poisson -0.5604432 0.0006827818 0.01890764      8.984185 54.63361

En la práctica, las funciones de enlace logit y probit son muy similares. Como se concluyó, ambas brindan un buen ajuste, en tanto que el modelo Poisson presenta un deviance mucho más bajo pero su AIC es superior al de los demás modelos, por lo que se podría decir que es mejor usar un modelo logit o probit.

EJERCICIO 4.7 (MYERS, MONTGOMERY, VINING & ROBINSON, 2010): Un importante fabricante de aviones estudió el número de fallas de un sujetador de aleación después de que cada sujetador se sometió a una carga de presión. Los datos son los siguientes.

Tabla 4.7. Muestra el número de fallas, T. de sujetadores y presión
presion	sujeta	fallas
2500	50	10
2700	70	17
2900	100	30
3100	60	21
3300	40	18
3500	85	43
3700	90	54
3900	50	33
4100	80	60
4300	65	51

Utilice el enlace logit y la regresión logística para analizar estos datos.

Inicialmente notamos que la variable respuesta $y_{i}$ es una variable conteo. Si deseamos realizar un análisis de regresión logística, esta variable debe tener distribución binomial, por lo que es necesario una transformación, para que represente la proporción de fallas por número de sujetadores; y además se toma:$g(\mu_{i})=log(\frac{p}{1-p})$

tur.m1 <- glm( fallas/sujeta ~ presion, family=binomial,
               weights=sujeta, data=dato1)

summary(tur.m1)

## 
## Call:
## glm(formula = fallas/sujeta ~ presion, family = binomial, data = dato1, 
##     weights = sujeta)
## 
## Deviance Residuals: 
##      Min        1Q    Median        3Q       Max  
## -0.29475  -0.11129   0.04162   0.08847   0.35016  
## 
## Coefficients:
##               Estimate Std. Error z value Pr(>|z|)    
## (Intercept) -5.3397115  0.5456932  -9.785   <2e-16 ***
## presion      0.0015484  0.0001575   9.829   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 112.83207  on 9  degrees of freedom
## Residual deviance:   0.37192  on 8  degrees of freedom
## AIC: 49.088
## 
## Number of Fisher Scoring iterations: 3

Del análisis se concluye que a un nivel de significancia del 5% los parámetros Presión e intercepto son significativos con un p-valor de 2e-16, hecho que que se confirma al observar los intervarlos de confianza, ya que ninguno contiene el cero. A continuación, calculamos la significancia del modelo logístico completo (p-valor del modelo) a través del test Likelihood ratio.

# Diferencia de residuos
dif_residuos <- tur.m1$null.deviance - tur.m1$deviance
# Grados libertad
df <-tur.m1$df.null - tur.m1$df.residual
# p-value
p_value <- pchisq(q = dif_residuos,df = df, lower.tail = FALSE)
Analisis=data.frame("Diferencia de residuos"=round(dif_residuos, 4),"Grados de libertad"=df,"p-value"=round(p_value, 4));Analisis

##   Diferencia.de.residuos Grados.de.libertad p.value
## 1               112.4602                  1       0

Para ello, construirmos un estadístico con la diferencia de residuos entre el modelo 1 (modelo con predictores) y el modelo 2 (modelo sin predictores). El estadístico anterior tiene una distribución chi-cuadrado con grados de libertad equivalentes a m1−m2,es decir, la diferencia de grados de libertad entre el modelo 1 y el modelo 2. Dado que el p-valor es 0, comparado con el nivel de significancia afirmamos que el modelo 1 es altamente significativo. De manera similar a lo hecho con el modelo logit, procedemos con el modelo probit. En este caso, el estadístico con distribución chi-cuadrado corresponde a 102.8552 con un grado de libertad y un p-valor de 0, que comparado con el nivel de significancia del 0.05 afirma que el modelo completo es significativo.

Repita este análisis utilizando el enlace probit.

En este caso, empleamos en la función de enlace probit $\eta_{i} =\phi^{-1} (\mu_{i})=probit(\mu_{i})$, donde $\phi(.)$ es la función de densidad acumulada para la distribución normal. Esta función de enlace se específica en R como link=“probit”.

tur.m2<- glm( fallas/sujeta ~ presion, family=binomial(link="probit"),
               weights=sujeta, data=dato1)

summary(tur.m2)

## 
## Call:
## glm(formula = fallas/sujeta ~ presion, family = binomial(link = "probit"), 
##     data = dato1, weights = sujeta)
## 
## Deviance Residuals: 
##      Min        1Q    Median        3Q       Max  
## -0.33473  -0.11319   0.01759   0.09600   0.36474  
## 
## Coefficients:
##               Estimate Std. Error z value Pr(>|z|)    
## (Intercept) -3.271e+00  3.213e-01  -10.18   <2e-16 ***
## presion      9.488e-04  9.281e-05   10.22   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 112.83207  on 9  degrees of freedom
## Residual deviance:   0.43735  on 8  degrees of freedom
## AIC: 49.154
## 
## Number of Fisher Scoring iterations: 3

Del análisis obtenido luego de emplear la función de enlace “probit”, se concluye que a un nivel de significancia del 5% los parámetros presión e intercepto son significativos con un p-valor de 2e-16, hecho que se confirma al observar los intervarlos de confianza, ya que ninguno contine el cero. A continuación, calculamos la significancia del modelo logístico completo con probit (p-valor del modelo), a través del test Likelihood ratio.

# Diferencia de residuos
dif_residuos <- tur.m2$null.deviance - tur.m2$deviance
# Grados libertad
df <-tur.m2$df.null - tur.m2$df.residual
# p-value
p_value <- pchisq(q = dif_residuos,df = df, lower.tail = FALSE)
Analisis=data.frame("Diferencia de residuos"=round(dif_residuos, 4),"Grados de libertad"=df,"p-value"=round(p_value, 4));Analisis

##   Diferencia.de.residuos Grados.de.libertad p.value
## 1               112.3947                  1       0

Evaluación del modelo

De manera similar a lo hecho con el modelo logit, procedemos con el modelo probit. En este caso el estadístico con distribución chi-cuadrado corresponde a 112.3947, con un grado de libertad; y un p-valor de 0, que comparado con el nivel de significancia del 0.05 afirma que el modelo completo es significativo.

Comente los dos análisis.

&nbsp

##        Intercept      presion Residual dev.      AIC
## Logit  -5.339712 0.0015484336     0.3719169 49.08842
## Probit -3.271209 0.0009488532     0.4373459 49.15385

Parte c: graficos de predicción para el modelo probit y logit Para comparar el modelo logístico completo, resultante de usar la función de enlace “logit” y “probit”, tenemos en cuenta los criterios Akaike y Deviance definidos con anterioridad.

En este caso, inferimos que los modelos de regresión “logit” y "probit proporcionan un ajuste similar, debido a que el modelo logit tiene menor deviance pero mayor AIC que el modelo probit, y el modelo probit tiene mayor deviance pero menor AIC que el modelo logit, cabe resaltar que la diferencia entre los valores de los criterios no es muy significativa, lo que genera dificultad para decidir cuál modelo es el más adecuado.

El análisis gráfico confirma lo dicho anteriormente, ambos modelos proporcionan un ajuste de los datos similar.

La distribución límite para el binomio es el de Poisson. Repita el inciso a) utilizando la regresión de Poisson. Discutir las diferencias

Como investigadores, obtenemos dos modelos, un modelo 1 en donde consideramos la presión como variable predictora, y un modelo dos en donde consideramos presión y sujeta como variables influyentes.

Modelo 1

## 
## Call:
## glm(formula = fallas ~ presion, family = poisson, data = dat1)
## 
## Deviance Residuals: 
##    Min      1Q  Median      3Q     Max  
## -2.110  -1.452  -0.495   1.494   2.347  
## 
## Coefficients:
##              Estimate Std. Error z value Pr(>|z|)    
## (Intercept) 0.8584212  0.3690003   2.326     0.02 *  
## presion     0.0007549  0.0001003   7.530 5.09e-14 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 84.993  on 9  degrees of freedom
## Residual deviance: 25.052  on 8  degrees of freedom
## AIC: 81.243
## 
## Number of Fisher Scoring iterations: 4

Se concluye que a un nivel de significancia del 5%, el parámetro Presión es significativo con un p-valor de 5.09e-14, y un deviance de 25.052.

# Diferencia de residuos
dif_residuos1 <- tr.poiss$null.deviance - tr.poiss$deviance
# Grados libertad
df1 <-tr.poiss$df.null - tr.poiss$df.residual
# p-value
p_value <- pchisq(q = dif_residuos1,df = df1, lower.tail = FALSE)
Analisis=data.frame("Diferencia de residuos"=round(dif_residuos1, 4),"Grados de libertad"=df,"p-value"=round(p_value, 4));Analisis

##   Diferencia.de.residuos Grados.de.libertad p.value
## 1                59.9408                  1       0

En este caso, el estadístico con distribución chi-cuadrado corresponde a 59.9408 con un grado de libertad y un p-valor de 0, que comparado con el nivel de significancia del 0.05 afirma que el modelo completo es significativo.

Gráficos de predicción para el modelo probit, logit y Poisson

El análisis gráfico confirma que con una sola variable predictiva, el modelo proporciona un buen ajuste de los datos.

Modelo 2

## 
## Call:
## glm(formula = fallas ~ presion + sujeta, family = poisson, data = dato1)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -0.5548  -0.2232   0.1045   0.1883   0.3617  
## 
## Coefficients:
##               Estimate Std. Error z value Pr(>|z|)    
## (Intercept) -0.3256388  0.4678987  -0.696    0.486    
## presion      0.0007703  0.0001031   7.471 7.98e-14 ***
## sujeta       0.0157207  0.0032516   4.835 1.33e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 84.99313  on 9  degrees of freedom
## Residual deviance:  0.93859  on 7  degrees of freedom
## AIC: 59.129
## 
## Number of Fisher Scoring iterations: 3

Se concluye que a un nivel de significancia del 5% los parámetros presión y sujeta son significativos con un p-valor de 7.98e-14 y 1.33e-06, respectivamente, y un deviance de 0.93859.

Para proceder con la comparación de los tres modelos (Modelo logístico “logit”, Modelo logistico “probit” y Modelo Poisson), es necesario elegir cuál modelo Poisson nos brinda un mejor ajuste. En este caso, elegimos el modelo 2 que presenta un AIC y Deviance menor.

Comparación de los tres modelos

En la práctica, las funciones de enlace logit y probit son muy similares, como se concluyó ambas brindan un buen ajuste, mientras que el modelo Poisson presenta un deviance mayor y AIC es superior al de los demás modelos, por lo que se podría decir que es mejor usar un modelo logit o probit.

tr.logit <-  glm( fallas/sujeta ~ presion, family=binomial,
                  weights=sujeta, data=dato1)
tr.probit <- update( tr.logit, family=binomial(link="probit") )
tr.poiss <-glm(fallas~presion+sujeta, family=poisson, data=dato1)
tr.array <- rbind( c(coef(tr.logit),0), c(coef(tr.probit),0), coef(tr.poiss))
tr.array <- cbind( tr.array, c(deviance(tr.logit), deviance(tr.probit), deviance(tr.poiss)),c((tr.logit$aic),(tr.probit$aic),(tr.poiss$aic)) )
colnames(tr.array) <- c("Intercept", "Hours","Turbines","Residual dev.","AIC")
rownames(tr.array) <- c("Logit","Probit","Poisson")
tr.array

##          Intercept        Hours  Turbines Residual dev.      AIC
## Logit   -5.3397115 0.0015484336 0.0000000     0.3719169 49.08842
## Probit  -3.2712088 0.0009488532 0.0000000     0.4373459 49.15385
## Poisson -0.3256388 0.0007703048 0.0157207     0.9385927 59.12939

EJERCICIO 4.8 (MYERS, MONTGOMERY, VINING & ROBINSON, 2010)

Maruthi y Joseph (1999-2000) llevaron a cabo un experimento para mejorar el rendimiento de circuitos de capa interna altos y densos en una operación de placa de circuito impreso. Los factores fueron (A) preparación de la superficie, (B) precalentamiento, (C) velocidad de laminación, (D) presión de laminación, (E) temperatura de laminación, (F) paso de exposición, (G) velocidad del revelador y (H) oxidación-reducción potencial. Midieron dos respuestas. La primera respuesta, y1 fue el porcentaje de cortos sobre 400 oportunidades. La segunda respuesta, y2, fue el porcentaje de apertura de 800 oportunidades. Los datos se muestran en el cuadro 4.8.

Cuadro 4.8. Resultados experimento Maruthi y Joseph (1999-2000)
A	B	C	D	E	F	G	H	y1	y2	P1	P2
-1	-1	-1	-1	-1	-1	-1	-1	0.260	0.122	26.0	12.2
-1	-1	0	1	0	0	0	0	0.190	0.166	19.0	16.6
-1	-1	1	-1	1	1	1	1	0.126	0.164	12.6	16.4
0	-1	-1	1	0	0	1	1	0.164	0.230	16.4	23.0
0	-1	0	-1	1	1	-1	-1	0.118	0.186	11.8	18.6
0	-1	1	1	-1	-1	0	0	0.169	0.077	16.9	7.7
1	-1	-1	-1	-1	1	0	1	0.128	0.178	12.8	17.8
1	-1	0	1	0	-1	1	-1	0.190	0.194	19.0	19.4
1	-1	1	-1	1	0	-1	0	0.175	0.205	17.5	20.5
-1	1	-1	1	1	0	0	-1	0.119	0.212	11.9	21.2
-1	1	0	-1	-1	1	1	0	0.098	0.172	9.8	17.2
-1	1	1	1	0	-1	-1	1	0.133	0.125	13.3	12.5
0	1	-1	-1	1	-1	1	0	0.169	0.186	16.9	18.6
0	1	0	1	-1	0	-1	1	0.116	0.133	11.6	13.3
0	1	1	-1	0	1	0	-1	0.092	0.205	9.2	20.5
1	1	-1	1	0	1	-1	0	0.075	0.169	7.5	16.9
1	1	0	-1	1	-1	0	1	0.212	0.200	21.2	20.0
1	1	1	1	-1	0	1	-1	0.164	0.172	16.4	17.2

Para ambas respuestas, individualmente, determine el modelo de regresión utilizando el enlace logit y logistic.

Notamos que la variable respuesta P1 y P2 son el porcentaje de cortos sobre 400 oportunidades y el porcentaje de apertura de 800 oportunidades, entonces construimos las variable y1 e y2 que equivalen a la expresión decimal de estos porcentajes, respectivamente. Para esta primera parte tomemos la función de enlace $g(\mu_{i})=log(\frac{y}{1-y})$.

Modelo logistico “Logit” con variable respuesta y1 e y2

Model1 <- glm(y1 ~A+B+C+D+E+F+G+H,family="binomial"("logit"),data=cr)
summary(Model1)

## 
## Call:
## glm(formula = y1 ~ A + B + C + D + E + F + G + H, family = binomial("logit"), 
##     data = cr)
## 
## Deviance Residuals: 
##      Min        1Q    Median        3Q       Max  
## -0.08821  -0.03909  -0.01200   0.02168   0.11080  
## 
## Coefficients:
##              Estimate Std. Error z value Pr(>|z|)   
## (Intercept) -1.770950   0.680939  -2.601   0.0093 **
## A            0.009856   0.823608   0.012   0.9905   
## B           -0.143017   0.677241  -0.211   0.8327   
## C           -0.023041   0.825580  -0.028   0.9777   
## D           -0.093269   0.727315  -0.128   0.8980   
## E           -0.027537   0.857824  -0.032   0.9744   
## F           -0.353556   0.844413  -0.419   0.6754   
## G            0.032249   0.828268   0.039   0.9689   
## H           -0.021936   0.823398  -0.027   0.9787   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 0.288281  on 17  degrees of freedom
## Residual deviance: 0.054466  on  9  degrees of freedom
## AIC: 23.889
## 
## Number of Fisher Scoring iterations: 5

De manera general, vemos que el modelo logístico logit para la variable y1 no es el más apropiado. Se concluye que a un nivel de significancia del 5% los parámetros (A,B,C,D,E,F,G,H) no son significativos con p-valores (0.9905,0.8327,0.9777,0.8980,0.9744,0.6754,0.9689 y 0.9787), respectivamente.

Model2<-glm(y2 ~A+B+C+D+E+F+G+H,family="binomial"("logit"),data=cr)
summary(Model2)

## 
## Call:
## glm(formula = y2 ~ A + B + C + D + E + F + G + H, family = binomial("logit"), 
##     data = cr)
## 
## Deviance Residuals: 
##      Min        1Q    Median        3Q       Max  
## -0.14097  -0.02535   0.01265   0.02502   0.09419  
## 
## Coefficients:
##              Estimate Std. Error z value Pr(>|z|)  
## (Intercept) -1.588687   0.634490  -2.504   0.0123 *
## A            0.102635   0.772750   0.133   0.8943  
## B            0.019155   0.631261   0.030   0.9758  
## C           -0.088085   0.772365  -0.114   0.9092  
## D            0.004995   0.687143   0.007   0.9942  
## E            0.183934   0.814571   0.226   0.8214  
## F            0.116035   0.806251   0.144   0.8856  
## G            0.115885   0.774657   0.150   0.8811  
## H           -0.038032   0.772253  -0.049   0.9607  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 0.18346  on 17  degrees of freedom
## Residual deviance: 0.05962  on  9  degrees of freedom
## AIC: 24.819
## 
## Number of Fisher Scoring iterations: 5

De manera general, vemos que el modelo logístico logit para la variable y2 no es el más apropiado. Se concluye que a un nivel de significancia del 5% los parámetros (A,B,C,D,E,F,G,H) no son significativos con p-valores (0.8943,0.9758,0.9092,0.9942, 0.8214,0.8856, 0.8811 y 0.9607), respectivamente.

## Warning in eval(family$initialize): non-integer #successes in a binomial glm!

## Warning in eval(family$initialize): non-integer #successes in a binomial glm!

Resultados
	Residual dev.	AIC
ModeloLogit1(y1)	0.0544657	23.88891
ModeloLogit2(y2)	0.0596197	24.81898

Ambos modelos “logit” presentan una desviance similar, pero el modelo con la variable y1 presenta un AIC de 23.88891 y el modelo con la variable y2 presenta un AIC de 24.81898, de lo cual podemos inferir que ambos modelos tienen un AIC similar. (b) Repita este análisis utilizando el enlace probit.

Modelo logistico “Probit” con variable respuesta y1 e y2

ModelP1 <- glm(y1 ~A+B+C+D+E+F+G+H,family="binomial"("probit"),data=cr)
summary(ModelP1)

## 
## Call:
## glm(formula = y1 ~ A + B + C + D + E + F + G + H, family = binomial("probit"), 
##     data = cr)
## 
## Deviance Residuals: 
##      Min        1Q    Median        3Q       Max  
## -0.08980  -0.03352  -0.01272   0.02094   0.10556  
## 
## Coefficients:
##              Estimate Std. Error z value Pr(>|z|)   
## (Intercept) -1.052872   0.367434  -2.865  0.00416 **
## A            0.006085   0.447852   0.014  0.98916   
## B           -0.078303   0.367691  -0.213  0.83136   
## C           -0.012619   0.448338  -0.028  0.97755   
## D           -0.050586   0.396668  -0.128  0.89852   
## E           -0.016811   0.466205  -0.036  0.97124   
## F           -0.194891   0.461927  -0.422  0.67309   
## G            0.018578   0.449043   0.041  0.96700   
## H           -0.014475   0.447397  -0.032  0.97419   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 0.288281  on 17  degrees of freedom
## Residual deviance: 0.053196  on  9  degrees of freedom
## AIC: 23.889
## 
## Number of Fisher Scoring iterations: 5

De manera general, vemos que el modelo logístico probit para la variable y1 tampoco es el más apropiado. Se concluye que a un nivel de significancia del 5% los parámetros (A,B,C,D,E,F,G,H) no son significativos con p-valores (0.98916, 0.83136, 0.97755, 0.89852, 0.97124, 0.67309, 0.96700 y 0.97419), respectivamente.

ModelP2<-glm(y2 ~A+B+C+D+E+F+G+H,family="binomial"("probit"),data=cr)
summary(ModelP2)

## 
## Call:
## glm(formula = y2 ~ A + B + C + D + E + F + G + H, family = binomial("probit"), 
##     data = cr)
## 
## Deviance Residuals: 
##      Min        1Q    Median        3Q       Max  
## -0.13736  -0.02565   0.01355   0.02532   0.09470  
## 
## Coefficients:
##              Estimate Std. Error z value Pr(>|z|)   
## (Intercept) -0.954118   0.351938  -2.711  0.00671 **
## A            0.055836   0.429970   0.130  0.89668   
## B            0.011704   0.352475   0.033  0.97351   
## C           -0.050077   0.430213  -0.116  0.90734   
## D            0.002316   0.382795   0.006  0.99517   
## E            0.102873   0.451367   0.228  0.81971   
## F            0.063946   0.448593   0.143  0.88665   
## G            0.063088   0.430580   0.147  0.88351   
## H           -0.020828   0.430041  -0.048  0.96137   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 0.183458  on 17  degrees of freedom
## Residual deviance: 0.058621  on  9  degrees of freedom
## AIC: 24.819
## 
## Number of Fisher Scoring iterations: 5

De manera general, vemos que el modelo logístico probit para la variable y2 tampoco es el más apropiado. Se concluye que a un nivel de significancia del 5% los parámetros (A,B,C,D,E,F,G,H) no son significativos con p-valores (0.89668, 0.97351, 0.90734, 0.99517, 0.81971, 0.88665, 0.88351 y 0.88351), respectivamente.

Resultados
	Residual dev.	AIC
MddeloProbit1(y1)	0.0531964	23.88890
ModeloProbit2(y2)	0.0586215	24.81901

Ambos modelos “probit” presentan un desviance similar, pero el modelo con la variable y1 presenta un AIC de 23.88890 y el modelo con la variable y2 presenta un AIC de 24.81901, de lo cual podemos inferir que ambos modelos tienen un AIC parecido. (c) Comente los dos análisis. De los dos análisis, inferimos que los modelos de regresión “logit” y "probit no son significativos, ambos proporcionan un ajuste similar, no hay diferencias entre los valores de los criterios (Deviance y AIC) y se recomienda emplear otro tipo de modelo para obtener un buen ajuste de los datos.

La distribución límite para el binomio es el de Poisson. Repita el inciso a) utilizando la regresión de Poisson. Discuta las diferencias.

Modelo Poisson para las variables P1 y P2

Modelo Completo para P1

ModeloPoiss1 <- glm(P1 ~A+B+C+D+E+F+G+H,family=poisson(link = "log"), data=cr)
summary(ModeloPoiss1 )

## 
## Call:
## glm(formula = P1 ~ A + B + C + D + E + F + G + H, family = poisson(link = "log"), 
##     data = cr)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -0.7994  -0.3968  -0.1091   0.2049   1.0516  
## 
## Coefficients:
##              Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  2.670783   0.063091  42.332  < 2e-16 ***
## A            0.008565   0.075843   0.113 0.910081    
## B           -0.121701   0.062498  -1.947 0.051499 .  
## C           -0.017284   0.076070  -0.227 0.820256    
## D           -0.077547   0.066802  -1.161 0.245707    
## E           -0.019795   0.079119  -0.250 0.802434    
## F           -0.298211   0.077429  -3.851 0.000117 ***
## G            0.028843   0.076378   0.378 0.705706    
## H           -0.014399   0.075943  -0.190 0.849623    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 24.4222  on 17  degrees of freedom
## Residual deviance:  4.6569  on  9  degrees of freedom
## AIC: Inf
## 
## Number of Fisher Scoring iterations: 4

Se concluye que a un nivel de significancia del 5% el intercepto y el parámetro F son significativos con p-valores (2e-16, 0.000117), respectivamente. La variable F es la única que tiene efecto sobre la variable respuesta y1, entonces construimos nuevamente un modelo Poisson, pero empleando solamente F (paso de exposición) como variable regresora.

Modelo Reducido para P1

ModeloRed1<- glm(P1 ~F,family=poisson(link = "log"), data=cr)
summary(ModeloRed1 )

## 
## Call:
## glm(formula = P1 ~ F, family = poisson(link = "log"), data = cr)
## 
## Deviance Residuals: 
##      Min        1Q    Median        3Q       Max  
## -1.45054  -0.56954   0.07631   0.46013   1.44439  
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  2.68147    0.06246  42.928  < 2e-16 ***
## F           -0.27933    0.07602  -3.674 0.000238 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 24.422  on 17  degrees of freedom
## Residual deviance: 10.657  on 16  degrees of freedom
## AIC: Inf
## 
## Number of Fisher Scoring iterations: 4

confint(ModeloRed1)

## Waiting for profiling to be done...

##                  2.5 %     97.5 %
## (Intercept)  2.5563522  2.8013327
## F           -0.4294837 -0.1312168

Del análisis se concluye que a un nivel de significancia del 5% los parámetros F (paso de exposición) e intercepto son significativos con p-valores de 0.000238 y 2e-16, respectivamente, hecho que se confirma al observar los intervarlos de confianza, ya que ninguno contiene el cero. A continuación, calculamos la significancia del modelo logístico completo (p-valor del modelo), a través del test Likelihood ratio.

Evaluación del Modelo

# Diferencia de residuos
dif_residuo <- ModeloRed1$null.deviance -ModeloRed1$deviance
# Grados libertad
d <-ModeloRed1$df.null - ModeloRed1$df.residual
# p-value
p_value <- pchisq(q = dif_residuo,df = d, lower.tail = FALSE)
Analisis=data.frame("Diferencia de residuos"=round(dif_residuo, 4),"Grados de libertad"=d,"p-value"=round(p_value, 4));Analisis

##   Diferencia.de.residuos Grados.de.libertad p.value
## 1                13.7656                  1   2e-04

En este caso, el estadístico con distribución chi-cuadrado corresponde a 13.7656, con un grado de libertad; y un p-valor de 2e-04, que comparado con el nivel de significancia del 0.05 implica que el modelo Poisson con una variable regresora P1 es significativo.

Modelo Completo para P2

ModeloPoiss2 <- glm(P2 ~A+B+C+D+E+F+G+H,family=poisson(link = "log"), data=cr)
summary(ModeloPoiss2)

## 
## Call:
## glm(formula = P2 ~ A + B + C + D + E + F + G + H, family = poisson(link = "log"), 
##     data = cr)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -1.3465  -0.2317   0.1161   0.2307   0.8413  
## 
## Coefficients:
##              Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  2.827089   0.057911  48.818   <2e-16 ***
## A            0.087202   0.070376   1.239   0.2153    
## B            0.015667   0.057325   0.273   0.7846    
## C           -0.072859   0.070199  -1.038   0.2993    
## D            0.004174   0.062461   0.067   0.9467    
## E            0.153593   0.074346   2.066   0.0388 *  
## F            0.098952   0.073310   1.350   0.1771    
## G            0.098265   0.070550   1.393   0.1637    
## H           -0.031929   0.070233  -0.455   0.6494    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 15.4574  on 17  degrees of freedom
## Residual deviance:  5.0995  on  9  degrees of freedom
## AIC: Inf
## 
## Number of Fisher Scoring iterations: 4

Se concluye que a un nivel de significancia del 5% el intercepto y el parámetro E son significativos con p-valores de 2e-16 y 0.0388, respectivamente. La variable E es la única que tiene efecto sobre la variable respuesta y1, entonces construimos nuevamente un modelo Poisson, pero empleando solamente E (temperatura de laminación) como variable regresora.

Modelo Reducido para P2

ModeloRed2<- glm(P2 ~E,family=poisson(link = "log"), data=cr)
summary(ModeloRed2 )

## 
## Call:
## glm(formula = P2 ~ E, family = poisson(link = "log"), data = cr)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -2.0267  -0.3570   0.0061   0.5997   1.3599  
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  2.83788    0.05723  49.584   <2e-16 ***
## E            0.14537    0.06997   2.078   0.0378 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 15.457  on 17  degrees of freedom
## Residual deviance: 11.118  on 16  degrees of freedom
## AIC: Inf
## 
## Number of Fisher Scoring iterations: 4

confint(ModeloRed2)

## Waiting for profiling to be done...

##                   2.5 %    97.5 %
## (Intercept) 2.723530151 2.9479641
## E           0.008574888 0.2830824

Del análisis se concluye que a un nivel de significancia del 5%, los parámetros E (temperatura de laminación) e intercepto son significativos con p-valores de 0.0378 y 2e-16, respectivamente, hecho que que se confirma al observar los intervarlos de confianza ya que ninguno contiene el cero. A continuación, calculamos la significancia del modelo logístico completo (p-valor del modelo), a través del test Likelihood ratio.

Evaluación del Modelo

# Diferencia de residuos
dif_residu <- ModeloRed2$null.deviance -ModeloRed2$deviance
# Grados libertad
d <-ModeloRed2$df.null - ModeloRed2$df.residual
# p-value
p_value <- pchisq(q = dif_residu,df = d, lower.tail = FALSE)
Analisis=data.frame("Diferencia de residuos"=round(dif_residu, 4),"Grados de libertad"=d,"p-value"=round(p_value, 4));Analisis

##   Diferencia.de.residuos Grados.de.libertad p.value
## 1                  4.339                  1  0.0372

En este caso, el estadístico con distribución chi-cuadrado corresponde a 4.339 con un grado de libertad y un p-valor de 0.0372, que comparado con el nivel de significancia del 0.05 implica que el modelo Poisson con una variable regresora P1 es significativo.