Modelo de regresión para razones y proporciones: variantes y aplicaciones socio-económicas y biosanitarias

Modelo de Regresión para razones y proporciones: Variantes y aplicaciones socio-económicas y biosanitarias Trabajo Fin de Grado Dept. de Estadística e Investigación Operativa FACULTAD DE MATEMÁTICAS Curso 2014-2015 Autor: Ana María Real Pérez Tutor: Juan Manuel Muñoz Pichardo . Abstract Many researchers use regression models to analyse data and study the relationship between some variables. The linear regression model, in particular, is commonly used in applications. This paper proposes a model based on the assumption that the response is beta distributed, which is very flexible for modelling proportions. For this reason, this model is appropriate for situations where the response is restricted to the interval (0, 1). The goal is to collect the modelling and inferential procedures for beta regression, including a description of betareg package which provides the class of beta regressions in the R system for statistical computing. Furthermore, some practical applications that employ real data are mentioned and described at the end. . Índice general 1. Distribución beta 7 1.1. Introducción............................ 7 1.2. Caracterización.......................... 10 1.3. Propiedades y relaciones . . . . . . . . . . . . . . . . . . . . . 12 1.4. Interpretación de parámetros. Otras parametrizaciones . . . . 13 1.5. Aplicaciones............................ 17 2. Regresión beta 19 2.1. Elmodelo............................. 19 2.2. Estimación de parámetros . . . . . . . . . . . . . . . . . . . . 21 2.3. Contrastes de hipótesis . . . . . . . . . . . . . . . . . . . . . . 25 2.4. Medidasdeajuste ........................ 27 2.5. Otras variantes del modelo . . . . . . . . . . . . . . . . . . . 28 3. Aplicaciones con R 29 3.1. El paquete Betareg ....................... 29 3.2. Aplicaciones prácticas . . . . . . . . . . . . . . . . . . . . . . 38 5 Capítulo 1 Distribución beta 1.1. Introducción La distribución beta es una familia de distribuciones de probabilidad continuas definidas en el intervalo [0,1] parametrizadas por dos parámetros, denotados αyβ, que controlan la distribución de la variable aleatoria. La distribución beta es usada para modelar el comportamiento de variables aleatorias limitadas por intervalos de longitud finita. En particular, es una distribución adecuada para porcentajes y proporciones. Un estudio realizado por Robert Kieschnick, de la Universidad de Texas, y BD McCullough, de la Universidad Drexel [1], compara siete modelos de regresión básicos para proporciones observadas en el intervalo abierto (0,1). Estos modelos son: 1. Distribución normal: función respuesta lineal En este modelo se asume que la distribución condicionada es una distribución normal y además la esperanza condicionada es una función lineal de unas variables predictoras x, es decir, f(y|x)∼N(k(x), σ2)con k(x) = x0β 2. Distribución normal logística aditiva Usando la transformación logit en la variable dependiente, ln(y/(1 − y)), tenemos el modelo ln y 1−y!=x0β+ 7 donde se asume que se distribuye según una normal estándar. Si la variable ysigue una distribución normal logística aditiva, se puede demostrar que la variable z= ln(y/(1 −y)) sigue una distribución normal. 3. Distribución normal censurada Modelo en el cual definimos la variable yi=       0, y∗ i≤0 y∗ i,0< y∗ i<1 1, y∗ i≥1 asumiendo que y∗ ise distribuye según una normal y que se puede escribir como y∗ i=x0 iβ+ui, i = 1,2, . . . , n donde uise consideran independientes e idénticamente distribuidas según una normal estándar. 4. Distribución normal: función respuesta no lineal Para el estudio de la función esperanza condicionada muchos investigadores recomiendan usar la función logística acumulada, que viene dada por: yi=1 1 + e−(α+βxi)+i donde isigue una distribución normal estándar. 5. Distribución beta Consideramos que la variable yestá distribuida de la siguiente forma: f(y) = 1 B(p, q)yp−1(1 −y)q−1 donde 0≤y≤1yB(p, q)es la función beta. Para el modelo de regresión beta podemos seguir el modelo de regresión lineal y asumir que la media es una función lineal de otras variables independientes x. Es decir, asumimos que la variable ysigue una distribución beta y que E(y|x) = p p+q=x0β 8 Otro modelo sería usar la función vínculo logit, que consiste en considerar E(yi|xi) = 1 1 + exp(−x0 iβ) 6. Distribución simple Otro modelo de regresión para proporciones se basa en la distribución simple, que viene dada por f(y;µ, σ2) = h2πσ2{y(1 −y)}3i−1/2exp (−1 2d(y;µ)) para 0< y < 1, donde d(y;µ) = (y−µ)2 y(1 −y)µ2(1 −µ)2 es la desviación, y 0<µ<1. El modelo lineal generalizado puede ser aplicado a modelos basados en la distribución simple, usando como función vínculo la función logit. 7. Modelo de regresión quasi-paramétrico Todos los casos anteriores eran modelos de regresión paramétricos, por lo que el último modelo es el modelo de regresión de quasi-verosimilitud para proporciones del intervalo (0,1) usando un número finito de observaciones límite, es decir, 0s y 1s: li(β) = yiln[G(xi)] + (1 −yi) ln[1 −G(xi)] con G(·)una función de las xique verifique 0< G(·)<1. En el estudio se utiliza la función logística en los casos prácticos. Una vez presentados los modelos de regresión más usados para el análisis de proporciones, éstos son aplicados a datos reales para poder comparar las ventajas y desventajas de cada modelo. A partir de este análisis comparativo, Kieschnick y McCullough concluyen su estudio recomendando a los investigadores usar en el futuro un modelo de regresión paramétrico basado en la distribución beta. Este capítulo se estructura en 5 secciones. En la segunda sección se presentan distintas caracterizaciones de la distribución beta , en la sección 1.3 se recogen algunas propiedades de interés de la distribución, en la siguiente 9 α≥1,β < 1 •Forma de J con cola a la izquierda •Coeficiente de asimetría negativo •Distribución estrictamente creciente y convexa •moda = 1 •1/2< mediana < 1 •0< var(X)<−11+5√5 2 α= 1,β > 1 •Coeficiente de asimetría positivo •Distribución estrictamente decreciente •Simétrica a la distribución función potencia, f(x) = xken el intervalo [0,1] •moda = 0 •Si α= 1 y1< β < 2la función es cóncava con 1−1/√2< mediana < 1/2y1/18 < var(X)<1/12 •Si α= 1 yβ= 2 es una recta de pendiente -2 con mediana = 1−1/√2yvar(X)=1/18 •Si α= 1 yβ > 2la distribución tiene forma de J con cola a la derecha, es convexa con 0< mediana < 1−1/√2y0< var(X)< 1/18 α > 1,β= 1 •Coeficiente de asimetría negativo •Distribución estrictamente creciente igual a la distribución función potencia en [0,1] •moda = 1 •Si 1<< 2yβ= 1 la función es cóncava con 1/2< mediana < 1/√2y1/18 < var(X)<1/12 •Si α= 2 yβ= 1 es una recta de pendiente +2 con mediana = 1/√2yvar(X)=1/18 •Si α > 2yβ= 1 obtiene una forma de J con cola a la izquierda y convexa con 1/√2< mediana < 1y0< var(X)<1/18 16 Figura 1.2: Figura superior izquierda: distribución beta con β= 2,5yα variando de 0 a 9. Figura superior derecha: distribución beta con β= 5,5 yαvariando de 0 a 9. Figura inferior: distribución beta con β= 8 yα variando de 0 a 10. 1.5. Aplicaciones La distribución beta tiene importantes aplicaciones como en la teoría de los estadísticos ordenados. El k-ésimo estadístico ordenado de una muestra de tamaño nde una distribución continua uniforme sigue una distribución beta, es decir, X(k)∼Beta(k, n + 1 −k). Esta distribución también tiene uso en la inferencia Bayesiana, ya que produce una familia de distribuciones de probabilidad a priori conjugadas, como por ejemplo Beta(1,1) (distribución uniforme en el intervalo (0,1)), Beta(0,0), (distribución Bernoulli cuando α, β →0) y Beta(1/2,1/2). Finalmente, como se recoge en la introducción de este capítulo y se ilustra en el capítulo siguiente, la distribución beta juega un papel fundamental en la modelización de razones y proporciones. 17 . 18 Capítulo 2 Regresión beta 2.1. El modelo El modelo de regresión comúnmente tiene por objetivo modelar y analizar la relación de una variable objetivo con otras variables explicativas a través de un conjunto de datos observados. En particular, el modelo de regresión lineal se usa con frecuencia, sin embargo no es apropiado para situaciones en la que la variable respuesta está restringida al intervalo (0,1). Una posible solución es transformar la variable dependiente para que tome valores en la recta real, y entonces modelar la media de la variable respuesta transformada. Esta aproximación tiene algunos inconvenientes como el hecho de que los parámetros del modelo no pueden ser fácilmente interpretados en términos de la variable original. El objetivo es proponer un modelo de regresión apto para situaciones donde la variable dependiente yes continua en el intervalo unidad, es decir, 0< y < 1. Este capítulo de estructura en 5 secciones. En la primera describimos la estructura del modelo, la inferencia correspondiente se desarrolla en la segunda sección, en la siguiente se presentan los contrastes de hipótesis y las medidas de ajuste en la sección 1.4. Por último, se presenta una posible extensión del modelo clásico. La bibliografía utilizada para el desarrollo de este capítulo es: Ferrari yCribari-Neto (2004). “Beta Regression for Modelling Rates and Proportions”[4] Cribari-Neto yQueiroz (2014). “On testing inference in beta regressions”[5] Hipótesis distribucional 19 El modelo se basa en asumir que la variable objetivo yse distribuye según la distribución beta, la cual es apropiada para modelar proporciones, como se ha recogido previamente en el Capítulo 1. Tomamos una nueva parametrización. Sea µ=α α+βyφ=α+βse tiene que α=µφ yβ= (1 −µ)φ. En consecuencia, la esperanza y la varianza de la variable objetivo se pueden escribir en función de estos nuevos parámetros de la siguiente forma: E[y] = µ V ar[y] = µ(1 −µ) 1 + φ donde µes la media de la variable respuesta y φpuede ser interpretado como el parámetro de precisión, en el sentido de, fijado µ, a mayor valor de φ, menor es la varianza de la variable. La función de densidad puede ser escrita con la nueva parametrización como: f(y, µ, φ) = Γ(φ) Γ(µφ)Γ((1 −µ)φ)yµφ−1(1 −y)(1−µ)φ−1,0< y < 1 donde 0< µ < 1yφ > 0. El modelo propuesto es también útil para situaciones donde la variable respuesta se restringe a intervalos genéricos (a,b), donde aybson escalares a < b. En este caso se modela con y−a b−aen vez de modelar con ydirectamente. Hipótesis estructural Sean y1, . . . , ynvariables aleatorias independientes, donde para cada yt, t= 1, . . . , n, se asume la hipótesis distribucional anterior con media µty parámetro de precisión φdesconocido. El modelo se obtiene asumiendo que la media de cada yt,µt, es vinculado al predictor lineal, ηt, como sigue: g(µt) = Pk i=1 xtiβi=x> tβ=ηt donde β= (β1, . . . , βk)>es un vector desconocido de parámetros de regresión (β∈R)y(xt1, . . . , xtk)es un vector de kcomponentes (k < n). Finalmente, g: (0,1) 7−→ Res una función vínculo estrictamente creciente y dos veces diferenciable. La varianza de ytes una función de µty, en consecuencia, de los valores covariantes xt. Por tanto, la varianza de la variable respuesta puede ser incorporada al modelo. Algunas posibles elecciones para la función vínculo g(·)son: 20 Función logit:g(µ) = log (µ 1−µ) Función probit:g(µ)=Φ−1(µ), donde Φ(·)es la función de distribución de una variable aleatoria normal estándar. Función log-log:g(µ) = −log {−log(1 −µ)} Función log-log complementaria:g(µ) = log {−log(1 −µ)} En particular, la función vínculo logit se puede escribir: x> tβ= log (µ 1−µ)⇐⇒ µt=ex> tβ 1 + ex> tβ donde x> t= (xt1, . . . , xtk),t= 1, . . . , n. Aquí, los parámetros de regresión tienen una interpretación interesante. Supongamos que el valor de la i-ésima variable, xti, es incrementado en cunidades, manteniendo las demás variables con su valor original, y sea µ∗la media de ybajo los nuevos valores de las variables explicativas y µla media de ycon los valores originales. Entonces, es fácil comprobar que ecβi=µ∗/(1 −µ∗) µ/(1 −µ) de manera que ecβies igual a la odds ratio. Para el desarrollo de la siguiente sección, correspondiente a la estimación de parámetros, se ha elegido la función vínculo logit. 2.2. Estimación de parámetros La estimación de parámetros suele resolverse por el método de máxima verosimilitud. La función de log-verosimilitud basada en una muestra de n variables respuesta independientes distribuidas según una beta es: l(β, φ) = Pn t=1 lt(µt, φ) donde lt(µt, φ) = log Γ(φ)−log Γ(µtφ)−log Γ((1 −µt)φ)+(µtφ−1) log yt+ ((1 − µt)φ−1) log(1 −yt) La función score, U= (Uβ(β, φ), Uφ(β, φ))>, es obtenida derivando la función de log-versimilitud con respecto a los parámetros desconocidos del modelo (β, φ). 21 ∂l(β, φ) ∂βi = n X t=1 ∂lt(µt, φ) ∂µt dµt dηt ∂ηt ∂βi Notamos que dµt/dηt= 1/g0(µt)y que ∂lt(µt, φ) ∂µt =φ"log yt 1−yt−{ψ(µtφ)−ψ((1 −µt)φ)}# donde ψ(·)es la función digamma, recogida en el Capítulo 1. Si llamamos y∗ t= log {yt/(1 −yt)}yµ∗ t=ψ(µtφ)−ψ((1−µt)φ)se puede escribir la derivada de la función de log-verosimilitud de la siguiente forma: ∂l(β, φ) ∂βi =φ n X t=1 (y∗ t−µ∗ t)1 g0(µt)xti de la cual obtenemos la expresión matricial de la función score para β: Uβ(β, φ) = φX>T(y∗−µ∗) donde Xes una matriz n×kcuya t-ésima componente es x> t,y∗= (y∗ 1, . . . , y∗ n)>, µ∗= (µ∗ 1, . . . , µ∗ n)>yTes una matriz diagonal cuya expresión es T=diag 1/g0(µ1),...,1/g0(µn) Se puede demostrar de manera similar que la función score para φpuede ser escrita de la siguiente forma: Uφ(β, φ) = n X t=1 {µt(y∗ t−µ∗ t) + log(1 −yt)−ψ((1 −µt)φ) + ψ(φ)} El siguiente paso es obtener la matriz de información de Fisher. La segunda derivada de l(β, φ)con respecto a los parámetros βs viene dada por ∂2l(β, φ) ∂βi∂βj =Pn t=1 ∂ ∂µt ∂lt(µt, φ) ∂µt dµt dηt!dµt dηt ∂ηt ∂βj xti =Pn t=1 ∂2lt(µt, φ) ∂µ2 t dµt dηt +∂lt(µt, φ) ∂µt ∂ ∂µt dµt dηt!dµt dηt xtixtj Como E(∂lt(µt, φ)/∂µt)=0, se tiene E ∂2l(β, φ) ∂βi∂βj!= n X t=1 E ∂2lt(µt, φ) ∂µ2 t! dµt dηt!2 xtixtj Ahora, derivando la expresión anterior ∂lt(µt, φ)/∂µtse obtiene 22 ∂2lt(µt, φ) ∂µ2 t =−φ2ψ0(µtφ) + ψ0((1 −µt)φ) y por tanto, denotando wt=φ{ψ0(µtφ) + ψ0((1 −µt)φ)}1 (g0(µt))2, E ∂2l(β, φ) ∂βi∂βj!=−φ n X t=1 wtxtixtj En forma matricial, se escribe E ∂2l(β, φ) ∂β∂β>!=−φX>WX siendo W=diag {w1, . . . , wn}. La segunda derivada de l(β, φ)con respecto a algún βiyφpuede ser escrita como ∂2l(β, φ) ∂βi∂φ = n X t=1 "(y∗ t−µ∗ t)−φ∂µ∗ t ∂φ #1 g0(µt)xti Si E(y∗ t) = µ∗ ty∂µ∗ t/∂φ =ψ0(µtφ)µt−ψ0((1 −µt)φ)(1 −µt), se llega a E ∂2l(β, φ) ∂βi∂φ !=− n X t=1 ct 1 g0(µt)xti donde ct=φ{ψ0(µtφ)µt−ψ0((1 −µt)φ)(1 −µt)}. Sea c= (c1, . . . , cn)>, en notación matricial se escribe E ∂2l(β, φ) ∂βi∂φ !=−X>T c Por último, se obtiene ∂2l(β, φ)/∂φ2derivando la expresión de Uφ(β, φ) respecto de φ. Sea dt=ψ0(µtφ)µ2 t+ψ0((1−µt)φ)(1 −µt)2−ψ0(φ), se llega a E ∂2l(β, φ) ∂φ2!=− n X t=1 dt que en forma matricial se puede escribir como E ∂2l(β, φ) ∂φ2!=−tr(D) siendo D=diag {d1, . . . , dn}. Entonces la matriz de información de Fisher viene dada por 23 K=K(β, φ) = Kββ Kβφ Kφβ Kφφ! donde Kββ =φX>W X Kβφ =K> φβ =X>Tc Kφφ =tr(D) = Pn t=1 ψ0(µtφ)µ2 t+ψ0((1 −µt)φ)(1 −µt)2−ψ0(φ) Se verifica que bajo condiciones de regularidad, cuando el tamaño de muestra es suficientemente grande, se tiene el siguiente comportamiento asintótico de los estimadores de máxima verosimilitud ˆ β, ˆ φ: ˆ β ˆ φ!∼Nk+1 β φ!, K−1! siendo K−1la inversa de la matriz de información de Fisher K−1(β, φ) = Kββ Kβφ Kφβ Kφφ! donde Kββ =1 φ(X>WX)−1(Ik+X>Tcc>T>X(X>wX)−1 γφ ) Kβφ = (Kφβ)>=−1 γφ(X>WX)−1X>T c Kφφ =γ−1 donde Ikes la matriz identidad k×kyγ=tr(D)−φ−1c>T>X(X>WX)−1X>Tc. Los estimadores de máxima versimilitud de βyφ, es decir ˆ βyˆ φ, son solución del sistema de ecuaciones (Uβ(β, φ)=0 Uφ(β, φ)=0 Sin embargo, estas ecuaciones no son lineales, por lo que es necesario resolverlas por maximización numérica de la función de log-verosimilitud usando un algoritmo de optimización no lineal, como métodos Newton o Quasi-Newton; para más detalles, ver Nocedal y Wright (1999)[6]. Para ello, se requiere especificar un valor inicial para las iteraciones, como el estimador mínimo-cuadrático para βobtenido de una regresión lineal de las variables respuesta tranformadas g(y1), . . . , g(yn)en X, es decir, (X>X)−1X>z, donde z= (g(y1), . . . , g(yn))>. También se necesita un valor inicial para φ. A partir de la expresión de la varianza de las variables respuesta transformadas llegamos a la siguiente expresión de φ: 24 var(g(yt)) = µt(1 −µt) 1 + φ⇒φ=µt(1 −µt) var(g(yt)) −1 Observamos que var(g(yt)) ≈var {g(µt)+(yt−µt)g0(µt)}=var(yt){g0(µt)}2 y de ahí se tiene que var(yt)≈var(g(yt)) {g0(µt)}−2. Por tanto, un posible valor inicial para φes 1 nPn t=1 ˇµt(1 −ˇµt) ˇσ2 t−1 donde ˇµt=g−1(x> t(X>X)−1X>z) ˇσ2 t=ˇe>ˇe (n−k){g0(ˇµt)}2 con ˇe=z−X(X>X)−1X>zel vector de residuos de minimos cuadrados. 2.3. Contrastes de hipótesis Consideremos el contraste de hipótesis para el vector de parámetros β= (β1, . . . , βm)>, con m<k, (H0:β=β(0) H1:β6=β(0) donde β(0) = (β(0) 1, . . . , β(0) m)>es un vector de dimensión mdado. El estadístico basado en la log-verosimilitud viene dado por: w1=−2l(˜ β, ˜ φ)−l(ˆ β, ˆ φ) donde l(β, φ)es la función de log-verosimilituud y (˜ β>,˜ φ)>es el estimador de máxima verosimilitud de (β>, φ)>bajo H0. Valores grandes de este estadístico significaría una gran diferencia entre la máxima log-versimilitud con el máximo obtenido asumiendo como cierta la hipítesis nula, con lo que iría a favor del rechazo de H0. Bajo condiciones de regularidad y bajo H0se verifica que w1sigue asintóticamente una distribución chi-cuadrado con mgrados de libertad. Otros estadísticos usados para el análisis de este contraste son: 25 converged función lógica que indica si optim converge. call función original call. formula fórmula original. terms lista con los elementos "mean", "precision 2 "full"que contienen las condiciones de los respectivos modelos. levels lista con los elementos "mean", "precision" y "full"que contienen los niveles de las varibles categóricas. contrasts lista con los elementos "mean", "precision"que contienen los contrastes correspondientes a los niveles de los respectivos modelos. model estructura del modelo en caso de model=TRUE. yvector de la variable respuesta en caso de y=TRUE xlista con los elementos "mean", "precision" que contienen las matrices de los respectivos modelos, en caso de x=TRUE Comando beta.control Todos los parámetros en betareg son estimados por máxima versimilitud usando optim con opciones de control en betareg.control. Los valores iniciales pueden ser elegidos via start o estimados por lm.wfit. betareg.control(phi = TRUE, method = "BFGS", maxit = 5000, hessian = FALSE, trace = FALSE, start = NULL, fsmaxit = 200, fstol = 1e-8, ...) Argumentos: phi indica si el parámetro de precisión φdebe ser tratado como un parámetro del modelo (TRUE por defecto). method lista de caracteres que especifica el argumento method de optim 32 maxit número entero que especifica el argumento maxit (máximo número de iteraciones) de optim. trace argumento lógico o número entero que controla si la información durante el progreso de la optimización debe ser mostrada. hessian argumento lógico que controla si la matriz hessiana resultante de optim debe ser usada para la estimaciónde la matriz de covarianzas. Por defecto, es empleada. start vector opcional con valores iniciales para todos los parámetros, incluido el parámetro φ. fsmaxit número máximo de iteraciones en el método iterativo marcador de Fisher. fstol tolerancia numérica para la convergencia del método iterativo de Fisher. ... argumentos para optim. Plot en betareg El comando plot para objetos betareg produce varios tipos de gráficos de diagnóstico, incluyendo varios tipos de residuos, medidas de influencia, etc. plot(x, which = 1:4, caption = c(“Residuals vs indices of obs.", “Cook’s distance plot", "Generalized leverage vs predicted values", “Residuals vs linear predictor", "Half-normal plot of residuals", "Predicted vs observed values"), sub.caption = paste(deparse(x$call), collapse = " @"), main = @, ask = prod(par("mfcol")) <length(which) && dev.interactive(), ..., type = "sweighted2", nsim = 100, level = 0.9) Argumentos: xmodelo objeto de clase "betareg" which argumento numérico. Si sólo es necesario un subconjunto de gráficos, especifica un subconjunto de los números 1:6. caption caracteres. Comentario al pie de la gráfica. sub.caption caracteres. Título común bajo las figuras en caso de haber múltiples. main caracteres. Título de cada gráfica además del pie de imagen. ask argumento lógico. Si TRUE, se pregunta al usuario antes de cada gráfica. ... otros parámetros para representar funciones. type carácter que indica el tipo de residuo que se usa, ver residuals.betareg. nsim numérico. Número de simulaciones en half-normal plots. level numérico. Nivel de confianza en half-normal plots. 33 Predicciones en betareg Las predicciones en un modelo de regresión beta pueden darse en valores del intervalo (0,1) o en una escala del predictor lineal. En R, usamos la función predict. predict(object, newdata = NULL, type = c(response", "link", "precision", "variance", "quantile"), na.action = na.pass, at = 0.5, ...) Argumentos: object modelo objeto de clase "betareg" newdata opcionalmente, es una base de datos para buscar variables para predecir. Si se omite, se usa las observaciones originales. type carácter que indica el tipo de predicciones: correspondiente a la media de las variables respuesta ("response"), al predictor lineal ("link"), al parámetro de precisión ("precision"), a la varianza de las variables respuesta ("variance"), o a los cuantiles de la distribición de las variables respuesta ("quantile") na.action función que determina lo que hacer con los valores perdidos en newdata. at vector numérico que indica los cuantiles para predecir en caso de type = "quantile". Por defecto, at = 0.5. Base de datos GasolineYield. A continuación, se usa el paquete betareg para aplicar el modelo de regresión beta a una base de datos llamada GasolineYield, que contiene las proporciones de petróleo convertido en gasolina tras la destilación fraccionada. Este conjunto de datos esta formado por 32 observaciones de 6 variables: yield proporción de petróleo convertido el gasolina tras la destilación fraccionada. gravity gravedad del petróleo (en grados API). pressure presión del petróleo (lbf/in2). temp10 temperatura (en grados F) en la cual el 10 % del petróleo se ha vaporizado. temp temperatura (en grados F) en la cual todo el petróleo se ha vaporizado. batch factor que indica el conjunto de condiciones gravity, pressure ytemp10. 34 La variable dependiente es la proporción de petróleo convertido en gasolina, siendo el resto las variables explicativas del modelo. Diez subconjuntos de valores de las tres primeras variables explicativas corresponden a diez tipos de petróleo distintos, los cuales fueron sujetos en la destilación a condiciones controladas experimentalmente. Estas condiciones se recogen en la variable batch. Código: >data("GasolineYield", package = "betareg") >gy <- betareg(yield gravity + pressure + temp10 + temp, data = GasolineYield) >summary(gy) Los resultados del código anterior son los siguientes: Call: betareg(formula = yield ∼gravity + pressure + temp10 + temp, data = GasolineYield) Standardized weighted residuals 2: Min 1Q Median 3Q Max -2.1189 -0.6985 -0.0088 0.6306 2.1572 Coefficients (mean model with logit link): Estimate Std. Error z value Pr(>|z|) (Intercept) -2.6949422 0.7625693 -3.534 0.000409 *** gravity 0.0045412 0.0071419 0.636 0.524871 pressure 0.0304135 0.0281007 1.082 0.279117 temp10 -0.0110449 0.0022640 -4.879 1.07e-06 *** temp 0.0105650 0.0005154 20.499 <2e-16 *** Phi coefficients (precision model with identity link): Estimate Std. Error z value Pr(>|z|) (phi) 248.24 62.02 4.003 6.26e-05 *** - - - Signif. codes: 0 ’***’ 0.001 ’**’ 0.01 ’*’ 0.05 ’.’ 0.1 ’ ’ 1 Type of estimator: ML (maximum likelihood) Log-likelihood: 75.68 on 6 Df Pseudo R-squared: 0.9398 Number of iterations: 144 (BFGS) + 5 (Fisher scoring) Por tanto,los resultados del modelo indica la significación de las variables 35 temp10 ytemp en la modelización de la variable yield. Y sin embargo, la no significación de las variables gravity ypressure. La representación gráfica de los resultados se obtiene a través de las siguientes órdenes: >par(mfrow = c(3, 2)) >plot(gy, which = 1:6) >par(mfrow = c(1, 1)) 36 Figura 3.1: Conjunto de gráficas generadas con el comando plot 37 3.2. Aplicaciones prácticas En esta sección se mencionan y describen algunos estudios reales relacionados con el modelo de regresión beta: Breton y otros (2014). “Prenatal Tobacco Smoke Exposure Is Associated with Childhood DNA CpG Methylation”[8] Fumar durante el embarazo puede provocar múltiples consecuencias negativas en la salud del niño. El objetivo de este artículo es relacionar modificaciones epigenéticas en niños asmáticos con exposición prenatal al tabaco. Para ello, usa el modelo de regresión beta con variable objetivo el porcentaje de metilación del ADN en una muestra de sangre y como variables predictoras medida de exposición al tabaco en el útero, edad, sexo, clínica de atención prenatal y tipo de célula estimada. Los autores de esta investigación concluyen que fumar durante el embarazo está relacionado con cambios epigenéticos que persisten durante la infancia. Swearingen y otros (2011). “Application of Beta Regression to Analyze Ischemic Stroke Volume in NINDS rt-PA Clinical Trials”[9] La cantidad de lesión ictus isquémica es difícil de analizar debido a la forma de su distribución subyacente. El estudio consiste en ver si la regresión beta es un buen modelo para estudiar esta variable. Para ello comparan los resultados obtenidos en este modelo con un análisis publicado con anterioridad. Finalmente, todas las variables previamente identificadas de importancia fueron detectadas en el modelo construido concluyendo así que el modelo de regresión beta es una alternativa legítima para analizar la cantidad de lesión provocada por un ictus isquémico. Hunger y otros (2011). “Analysis of SF-6D Index Data: Is Beta Regression Appropriate?”[10] El pino de corteza blanca (Pinus albicaulis) depende del cascanueces americano (Nucifraga columbiana) para la dispersión de sus semillas. Mientras el pino de corteza blanca sucumbe a una enfermedad exótica y a los escarabajos del pino de montaña (Dendroctonus ponderosae), disminuye la producción de piña, por lo que los cascanueces americanos los visitan con menos frecuencia, disminuyendo la probabilidad de dispersar sus semillas. Este estudio de investigación demuestra que 38 el modelo de regresión beta es el método estadístico más apropiado para modelar la relación entre la densidad de piña y la proporción del tiempo de estudio en el que se observan cascanueces americanos. 39 Bibliografía [1] Kieschnick yMcCullough (2003). “Regression analysis of variates observed on (0, 1): percentages, proportions and fractions”, Statistical Modelling 3, 193–213. [2] Johnson y otros (1995). “Chapter 21: Beta Distributions”, Continuous Univariate Distributions Vol. 2 (2nd ed.). Wiley. ISBN 978-0- 471-58494-0. [3] Wikipedia contributors (2015). “Beta distribution”, Wikipedia, The Free Encyclopedia. Wikipedia, The Free Encyclopedia, 28 Aug. 2015. Web. 5 Sep. 2015. [4] Ferrari yCribari-Neto (2004). “Beta Regression for Modelling Rates and Proportions”, Journal of Applied StatisticsJournal of Applied Statistics 31:7, 799-815, DOI: 10.1080/0266476042000214501. [5] Cribari-Neto yQueiroz (2014). “On testing inference in beta regressions”, Journal of Statistical Computation and Simulation 84:1, 186-203, DOI: 10.1080/00949655.2012.700456. [6] Nocedal yWright (1999). “Numerical Optimization”, New York: Springer-Verlag [7] Zeileis y otros (2013). “Package ‘betareg”’, Beta Regression, CRAN, 2013-09-16 18:59:12 [8] Breton y otros (2014). “Prenatal Tobacco Smoke Exposure Is Associated with Childhood DNA CpG Methylation”, PLoS ONE 9(6): e99716. [9] Swearingen y otros (2011). “Application of Beta Regression to Analyze Ischemic Stroke Volume in NINDS rt-PA Clinical Trials”, Neuroepidemiology, 37:73–82. [10] Hunger y otros (2011). “Analysis of SF-6D Index Data: Is Beta Regression Appropriate?”, Value in Health, 759–767 40