Regresión xeralizada aplicada

Traballo Fin de Grao Regresión Xeralizada Aplicada Alba Camino Enríquez Xullo, 2022 UNIVERSIDADE DE SANTIAGO DE COMPOSTELA GRAO DE MATEMÁTICAS Traballo Fin de Grao Regresión Xeralizada Aplicada Alba Camino Enríquez Xullo, 2022 UNIVERSIDADE DE SANTIAGO DE COMPOSTELA iii iv Traballo proposto Área de Coñecemento: Estatística e Investigación Operativa Título: Regresión Xeralizada Aplicada Breve descrición do contido Os modelos de regresión serven para explicar e modelar a relación que existe entre unha variable resposta e unha ou máis variables explicativas. Tomando como base o modelo de regresión lineal simple clásico, o obxectivo deste traballo será o de presentar distintas extensións que permitan xeralizar dito modelo e analizar o seu funcionamento. Este traballo de fin de grado terá un forte carácter aplicado, polo tanto, os modelos revisados serán seleccionados e empregados en función da natureza dos datos que se dispoñan. Recomendacións Faraway, J.J. (2006). Extending the Linear Model with R: Generalized Linear, Mixed Effects and Nonparametric Regression Models. Chapman and Hall. Sheather, S.J. (2009). A modern approach to regression with R. Springer. Outras observacións O tratamento dos datos realizarase empregando o software estadístico de uso libre R (https://www.r-project.org/) Índice Resumo viii Introdución xi 1. Modelo de Regresión Loxística 1 1.1. Estimación dos parámetros do modelo . . . . . . . . . . . . . . . . . . . . . . . . 3 1.1.1. Métodos iterativos para o cálculo das estimacións . . . . . . . . . . . . . . 5 1.2. Inferencia sobre os parámetros do modelo . . . . . . . . . . . . . . . . . . . . . . 7 1.2.1. Contraste de modelos mediante deviance ................... 7 1.2.2. Intervalos de confianza para os parámetros de regresión . . . . . . . . . . 8 1.3. SeleccióndoModelo .................................. 9 1.3.1. Detección de datos atípicos . . . . . . . . . . . . . . . . . . . . . . . . . . 10 1.4. Conclusión........................................ 11 2. Modelo de Poisson 13 2.1. Plantexando o Modelo de Poisson . . . . . . . . . . . . . . . . . . . . . . . . . . 14 2.2. Estimación dos parámetros do modelo . . . . . . . . . . . . . . . . . . . . . . . . 15 2.2.1. Métodos iterativos para o cálculo das estimacións . . . . . . . . . . . . . . 16 2.3. Inferencia sobre os parámetros do modelo . . . . . . . . . . . . . . . . . . . . . . 17 2.3.1. Contraste de modelos mediante deviance ................... 17 2.3.2. Intervalos de confianza para os parámetros de regresión . . . . . . . . . . 19 v vi ÍNDICE 2.4. SeleccióndoModelo .................................. 20 2.4.1. Detección de datos atípicos . . . . . . . . . . . . . . . . . . . . . . . . . . 20 2.5. Sobredispersión..................................... 20 2.6. Conclusión ....................................... 21 3. Aplicación do Modelo de Regresión Loxística 23 3.1. Diagnóstico do cancro de mama . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23 3.1.1. Descrición dos datos . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23 3.1.2. Aplicación do modelo de regresión loxística. . . . . . . . . . . . . . . . . . 25 3.1.3. Inferencia sobre o modelo. . . . . . . . . . . . . . . . . . . . . . . . . . . . 28 3.1.4. Diagnose sobre o modelo. . . . . . . . . . . . . . . . . . . . . . . . . . . . 31 3.1.5. Conclusión. ................................... 34 4. Aplicación do Modelo de Poisson 39 4.1. Demanda de atención médica . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39 4.1.1. Descrición dos datos . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39 4.1.2. Aplicación do modelo de Poisson. . . . . . . . . . . . . . . . . . . . . . . . 40 4.1.3. Inferencia sobre o modelo. . . . . . . . . . . . . . . . . . . . . . . . . . . . 44 4.1.4. Diagnose sobre o modelo. . . . . . . . . . . . . . . . . . . . . . . . . . . . 45 4.1.5. Sobredispersión do modelo. . . . . . . . . . . . . . . . . . . . . . . . . . . 47 4.1.6. Conclusión ................................... 47 Bibliografía 51 Appendices 53 2 1. Modelo de Regresión Loxística verifique 0≤g−1(η)≤1para calquera valor η. Ainda que podemos pensar en distintas funcións gque satisfagan estas propiedades, a elección máis popular é a función loxística ou función logit dada pola seguinte expresión: η=g(p) = ln p 1−p, ou equivalentemente p=expη 1 + expη.(1.2) A combinación do emprego desta función logit xunto cun predictor lineal é o que recibe o nome de Regresión Loxística. Figura 1.1: Relación entre a probabilidade da resposta, p, e o predictor lineal η. Na Figura 1.1 podemos apreciar que a curva loxística é case lineal no seu rango medio. Isto quere dicir, que para o modelaxe de respostas con probabilidades próximas a 0,5o comportamento da regresión loxística e lineal non vai ser moi diferente. Ademais pódese apreciar que nos extremos a curva achégase a cero e a un pero sen chegar nunca a tales límites, o que quere dicir que a regresión loxística non vai predicir algo inevitable ou imposible. Recapitulando, recordemos que pnas expresións anteriores fai referencia a probabilidade de éxito fronte a 1−pque fai referencia a de fracaso. A función loxística consiste en aplicar un logaritmo ao cociente de ambas probabilidades. O cociente desas probabilidades indicadas é o que coñecemos como odds. As odds defínense como a probabilidade de éxito entre a de fracaso odds =p 1−p 1.1. Estimación dos parámetros do modelo 3 ou de xeito equivalente p=odds 1 + odds. Estas son unha escala alternativa a probabilidade para representar o azar, pódense ver como unha forma de expresar os pagos das apostas. Así podemos dicir que a probabilidade de que o FC Barcelona gañe o clásico fronte ao Real Madrid é de 2/3ou o que é o mesmo que a odd vale 2. Falando en términos de apostas diríamos que as apostas están 2 a 1 a favor do FC Barcelona. Unha importante vantaxe matemática das odds é que poden tomar calquera valor real positivo, mentres que a probabilidade de éxito só podía tomar valores no intervalo [0,1]. Dado que obviamente ao ser un cociente de cantidades positivas as odds van a ser positivas podemos aplicarlle un logaritmo e así transformalas nunha cantidade real calquera. Empregando o concepto de odds e combinándoo con (1.1), o modelo de regresión loxística é: η= ln odds = ln p 1−p=β0+β1X1+β2X2+... +βqXq ou odds = expβ0·expβ1X1·expβ2X2·... ·expβqXq. De aquí podemos interpretar β1como segue: un aumento unitario en X1con X2, ... , Xqfixadas aumenta os ln(odds)de éxito de β1ou equivalentemente aumenta as probabilidades de éxito dun factor expβ1. Polo que en certo modo a interpretación dos coeficientes expoñenciais pode ser máis práctica. 1.1. Estimación dos parámetros do modelo A continuación imos estimar os parámetros do modelo de regresión loxística, para iso usaremos o xa coñecido método de máxima verosimilitude. O estimador de máxima verosimilitude é o valor ou valores dos parámetros que maximiza a función masa de probabilidade ou densidade da mostra nas observacións. Definindo a función de verosimilitude para unha distribución con función masa de probabilidade Pθque depende dun parámetro θ, como: α(θ) = Pθ(X1=x1, ... , Xn=xn) = n Y i=1 Pθ(X=xi), o estimador de máxima verosimilitude é aquel valor ˆ θpara o cal temos a seguinte igualdade: α(ˆ θ) = m´ax θα(θ). Como a función logaritmo é monótona crecente e podemos supoñer que α(θ)é positiva, ˆ θé o estimador de máxima verosimilitude se e só se é un máximo de: l(θ) = ln (α(θ)) .(1.3) 4 1. Modelo de Regresión Loxística A obtención dos candidatos a máximos pódese conseguir derivando (1.3) con respecto de θ. Céntrandonos agora, no caso que nos ocupa, sabemos que no modelo loxístico a variable resposta segue unha distribución de Bernoulli con parámetro pi, a súa función masa de probabilidade aplicada a unha observación yi, onde i∈ {1, ... , n}, é : P(Y=yi;pi) = pyi i(1 −pi)1−yi, e a súa función de verosimilitude ten a seguinte expresión: α(y,p) = n Y i=1 pyi i(1 −pi)1−yi= n Y i=1 pi 1−piyi (1 −pi). Aplicando as funcións expoñencial e logarítmo e empregando as propiedades destes, chegamos a seguinte expresión equivalente: α(y,p) = (n Y i=1 (1 −pi))(n Y i=1 exp ln pi 1−piyi). Tendo en conta que ln pi 1−piyi=yiln pi 1−pi=yiηiepi=exp ηi 1+exp ηie a continuación remplazando ηi= ln pi 1−pie1−pi=1 1+exp ηiobtemos a seguinte expresión simplificada: α(η) = (n Y i=1 1 1 + exp ηi)exp "n X i=1 yiηi#.(1.4) Empregando agora o logaritmo para (1.4) e a propiedade do logaritmo dun cociente obténse: l(η) = n X i=1 yiηi−ln (1 + exp ηi).(1.5) Partindo de (1.5) e tendo en conta que ηi=β0+β1xi1+... +βqxiq ou o que é o mesmo tomando ηi=x′ iβ, sendo x′ io trasposto do vector de prediccións xi. Teremos así, a seguinte expresión en función de β: l(β) = n X i=1 yix′ iβ−ln(1 + exp(x′ iβ)).(1.6) Plantexamos a continuación a expresión da derivada da log-verosimilitude. Deste xeito, igualándoas a cero, obteremos o parámetro que será candidato a ser o máximo, e polo tanto o estimador de β ∂l(β) ∂β= n X i=1 yix′ i−exp(x′ iβ) 1 + exp(x′ iβ)x′ i= 0.(1.7) En resumo, o método de máxima verosimilitude permítenos obter un estimador, para iso o que se fai é derivar a log-verosimilitude e igualar a cero, chegando a (1.7). Polo que, finalmente para obter o estimador desexado debemos resolver as ecuacións (1.7) e así chegar a unha expresión 1.1. Estimación dos parámetros do modelo 5 para βj. No caso do modelo lineal, as estimacións dos parámetros obtíñanse a partir do método de mínimos cadrados co cal acadabamos un sistema de necuacións e nincógnitas, coñecidas como ecuacións normais de regresión con solucións explícitas. Sen embargo, neste caso a ausencia de solución explícita, fai que teñamos que recurrir a métodos iterativos co fin de achar o estimador buscado. Para resolver numéricamente as ecuacións (1.7) é preciso recorrer a matriz hessiana, dada pola expresión: ∂2l(β) ∂β∂β′=− n X i=1 xi exp(x′ iβ)x′ i 1 + exp(x′ iβ)2!.(1.8) A partir desta expresión (1.8) e tendo en conta que, ηi=x′ iβ,pi=exp ηi 1+exp ηie1−pi=1 1+exp ηi podemos reescribir a expresión para a matriz hessiana como segue: ∂2l(β) ∂β∂β′=− n X i=1 [xix′ ipi(1 −pi)].(1.9) Dado que precisamos resolver as ecuacións plantexadas pode ser de gran utilidade escribir matricialmente a expresión obtida para a matriz hessiana. Tendo en conta a seguinte notación: X=    1x11 ... x1q . . . 1xn1... xnq     ,(1.10) onde esta expresión (1.10) non é máis que a matriz que contén os valores das variables explicativas. E denotando por V a matriz diagonal, con valores pi(1 −pi), temos: V=    p1(1 −p1) 0 ... 0 ... 0... 0pn(1 −pn)     , finalmente a matriz hessiana pódese escribir: H=∂2l(β) ∂β∂β′=− n X i=1 [xix′ ipi(1 −pi)] = −X′V X. (1.11) No caso que nos ocupa, a partir da expresión para a matriz hessiana (1.11) vemos que é preciso empregar a teoría estándar para resolver un sistema de n+1 ecuacións e n+1 incógnitas e así obter aproximacións dos erros estándares. Como non se coñece unha solución explícita, temos que recorrer ao uso de métodos iterativos. 1.1.1. Métodos iterativos para o cálculo das estimacións Para resolver as ecuacións de verosimilitude do modelo loxístico existen diversos métodos iterativos. O primeiro e máis coñecido para nós é o método de Newton-Raphson. Este método 6 1. Modelo de Regresión Loxística parte dun iterante inicial, β0, e obtén sucesivamente o valor do seguinte parámetro mediante a seguinte expresión (Hastie, Tibshirani, e Friedman, 2009): βk+1 =βk−∂2l(β) ∂β∂β′|β=βk−1∂l(β) ∂β|β=βk,(1.12) sendo ko número de iteración. Dada a expresión xeral para o método de Newton-Raphson simplemente queda sustituir en (1.12) para o caso particular que nos ocupa, é dicir, o relativo a regresión loxística. Para iso debemos recordar a expresión matricial para a matriz hessiana descrita en (1.11) e as relacións dadas en (1.1) e en (1.2) así como a expresión da matriz de variables explicativas (1.10). Neste punto podemos expresar as ecuacións de verosimilitude dadas en (1.7) con notación matricial.: ∂l(β) ∂β=X′(y−p)=0,(1.13) onde yé o vector y= (y1, ... , yn), é dicir é a mostra de Ye análogamente p= (p1, ... , pn). Tendo en conta estas expresións , volvemos a escritura xeral do método (1.12) e sustituimos polas expresións relativas a matriz hessiana e as ecuacións de verosimilitude ambas en forma matricial, obtendo: βk+1 =βk+ (X′VkX)−1X′(y−pk),(1.14) onde Vké a matriz diagonal Vcon valores pi(1 −pi)epké o vector pmudando βpolo valor de βna iteración anterior, é dicir por βk. A partir de (1.14) podemos comezar o proceso iterativo, para o cal, en primeiro lugar debemos partir dun iterante inicial β0e establecer un umbral ϵ. Este valor ϵpermitirános saber se o noso algoritmo converxe ou non. Principalmente o proceso iterativo a seguir consistirá no cálculo de pusando p=exp(Xβk) 1+exp(Xβk)e no cálculo de Vcuxa diagonal vén dada por Vii =pi(1 −pi). Estos termos calculados iránse introducindo progresivamente na expresión (1.14) para cada iteración. Repetíndose o proceso ata que ∥βk−βk+1∥< ϵ, cando isto ocorra se é que sucede, usaremos o valor obtido de βkcomo estimador do vector de parámetros βe rematamos o procedemento. Se isto último, non se da ao cabo dunha cantidade determinada de iteracións, podemos dicir que non hai converxencia. Sen embargo, a hora de levar a cabo a implementación deste método no software de R (R Core Team, 2021) debemos desenvolver manualmente o proceso descrito, o que pode resultar tedioso a hora de escribir o código. Para realizar a regresión loxística, o software R (R Core Team, 2021), na función glm(...,family="binomial" ), non emprega este método por defecto se non que utiliza o método de mínimos cadrados reponderados iterativamente IRLS polas súas siglas en inglés (iteratively reweighted least squares). O método IRLS está estreitamente relacionado co método de Newton-Raphson (Agresti, 1990). Para comprender esta relación compre obter a matriz de información de Fisher (Colaboradores 1.2. Inferencia sobre os parámetros do modelo 7 de Wikipedia, 2022), que para o noso caso concreto, esta coincide co oposto da matriz hessiana dada en (1.11): J=X′V X =−∂2l(β) ∂β∂β′.(1.15) A continuación volvendo a expresión xeral do método de Newton-Raphson (1.12) e neste caso escribindoo en términos da matriz de información de Fisher teremos: βk+1 =βk+J−1X′(y−p).(1.16) Agora multiplicando a ambos lados de (1.16) por J=X′V X , chegamos a: (X′V X)βk+1 = (X′V X)βk+ (X′V X)(X′V X)−1X′(y−p) X′VhXβk+V−1(y−p)i=X′Vzk, deste xeito temos que: βk+1 = (X′V X)−1X′Vzk,(1.17) onde, zk=Xβk+V−1(y−p).(1.18) Chegados a isto, xa temos todo o necesario para iniciar o proceso iterativo partindo de (1.17) e avanzando nos valores de Vezkdebido a que en cada iteración o valor de pvaría acorde aos valores de βcalculados na iteración anterior. Este algoritmo recibe o nome de mínimos cadrados reponderados iterativamente xa que en cada iteración resolvemos un problema de mínimos cadrados ponderados pola matriz V: arg m´ın βk+1=β (z−Xβ)′V(z−Xβ).(1.19) 1.2. Inferencia sobre os parámetros do modelo A continuación levaremos a cabo certos procesos de inferencia sobre o noso modelo de regresión loxística, como: realizar un contraste para ver se se pode asumir que os parámetros βjson iguais a cero e construir intervalos de confianza para ditos parámetros. 1.2.1. Contraste de modelos mediante deviance O obxectivo desta sección será obter o mellor axuste posible para o noso modelo. Con este fin, empregaremos o test de razón de verosimilitudes, este método basease en comparar a verosimilitude baixo unha condición, H0coa verosimilitude baixo o modelo con outra condición alternativa H1. Para o caso que nos ocupa e partindo da función de verosimilitude descrita en (1.4), podemos 8 1. Modelo de Regresión Loxística escribir esta última en términos dos βsimplemente tendo en conta que pi=exp(xi ′β) 1+exp(xi ′β). Obtendo así, a función de verosimilitude en términos das βcomo: α(β|(x1, y1)... (xn, yn)) = n Y i=1 exp (xi ′β) 1 + exp (xi ′β)yi1 1 + exp (xi ′β)1−yi .(1.20) Recorrendo ao método citado, dados dous modelos con distinto número de variables explicativas, un modelo con lvariables explicativas e función de verosimilitude αl, obtida empregando o estimador de βna expresión (1.20), e un modelo con svariables explicativas e verosimilitude αs. Onde o modelo de menor tamaño representa un subconxunto do grande, é dicir, son modelos aniñados (s > l). A seguinte relación estadística pódese empregar como un estadístico válido para comparar ambos modelos: −2 log αl αs .(1.21) A distribución deste estadístico pódese aproximar mediante unha ji-cadrado con tantos graos de liberdade como parámetros se perden ao pasar do modelo longo ao máis curto (s−l). Relacionado co estadístico (1.21) aparece o concepto de deviance, que mide a desviación do modelo loxístico axustado respecto dun modelo perfecto. Este modelo ideal coñécese como modelo saturado e é un modelo que proporciona predicións da variable resposta que coinciden cos propios valores da mostra. No caso de modelos con variable resposta discreta, como o que nos ocupa, é máis sinxelo obter un modelo saturado xa que a variable resposta ten menos valores posibles. Dado o modelo axustado e o modelo saturado a deviance coincide coa razón de verosimilitudes (1.21) cando se teñen en conta o modelo axustado e o modelo saturado: DModelo =−2 log Verosimilitude Modelo Verosimilitude Modelo Saturado.(1.22) Empregando por unha parte a expresión para a verosimilitude (1.20) e por outra a definición para a deviance do modelo (1.22) e tendo en conta que a verosimilitude do modelo saturado é 1, adeviance para o noso modelo loxístico é: D=−2 n X i=1 [yilog( ˆpi) + (1 −yi) log(1 −ˆpi)],(1.23) onde ˆpi=exp ˆηi 1+exp ˆηison os valores do modelo axustado, con ˆηi=ˆ β0+ˆ β1xi1+... +ˆ βqxiq. 1.2.2. Intervalos de confianza para os parámetros de regresión Consideraremos o intervalo de confianza de Wald (Cepeda-Cuervo et al., 2008) para o noso estimador ˆ βiobtido a partir do método de máxima verosimilitude descrito na Sección 1.1. Neste 1.3. Selección do Modelo 9 punto, acudindo a teoría asintótica sobre os estimadores de máxima verosimilitude ou a un teorema do límite central podemos deducir que a distribución de ˆ βconverxe a N(β,(X′V X)−1). É dicir, ˆ βsegue unha distribución normal de media βe matriz de varianzas-covarianzas J−1= (X′V X)−1. Facendo uso do resultado anterior, o intervalo de confianza para βié o seguinte: ˆ βi±zα/2ET (ˆ βi),(1.24) onde zα/2é o cuantil 1−α/2dunha distribución normal estándar e no tamaño mostral. Con isto, simplemente bastaría obter unha fórmula correcta para o erro típico de ˆ βi. O erro típico áchase partindo de que ˆ βisegue unha distribución normal cuxa matriz de varianzas-covarianzas coincide precisamente coa matriz de información de Fisher J, definida en (1.15). A continuación, estimando os parámetros βe tendo presente tamén que J=X′V X , en lugar de Vintroduciremos ˆ Vsustituíndo en Vos parámetros polas súas estimacións. Definindo agora unha matriz ˆ J=X′ˆ V X unha forma de obter o erro típico asociado a un βiconcreto consistirá en facer simplemente a raíz cadrada do elemento (i, i)da matriz ˆ J−1. Deste xeito, dado calquera estimador ˆ βi, poderemos obter o seu erro típico e por conseguinte a forma do intervalo de Wald do mesmo. Teóricamente este intervalo ten, para valores grandes de n, un nivel de confianza aproximado de 100(1 −α) % 1.3. Selección do Modelo Podemos descubrir que non todas as variables predictoras obtidas son útiles para explicar a resposta polo que tentaremos identificar un subconxunto destas que modele o mellor posible a resposta. Para isto podemos empregar un método de axuste de fácil implementación como o método de eliminación back-ward. Para a implementación deste, partimos do modelo completo, con todas as variables predictoras dispoñibles, comparamos secuencialmente este modelo con todos os modelos resultantes de eliminar cada unha destas variables predictoras. Para isto, podemos obter os p-valores correspondentes a realizar o contraste H0:βj= 0 fronte a H1:βj= 0. Unha estratexia habitual é eliminar aquela variable predictora cuxo p-valor asociado tome o valor máis elevado. Repetimos o proceso ata que non se poidan eliminar máis variables predictoras sen unha pérdida do axuste estadísticamente significativa. A significación de cada variable pódemola ver realizando no software de R (R Core Team, 2021) un summary do modelo glm(... , family="binomial") aplicado anteriormente e ollando o valor obtido do p-valor asociado a cada variable. Con family="binomial" especifícase que función de probabilidade utilizamos e glm por defecto emprega a función logit como función de enlace. A continuación o summary emprega a normalidade asintótica introducida na sección (1.2.2) para realizar o contraste mencionado. Deste xeito, se o p-valor obtido no summary é moi baixo indicanos que debemos rexeitar a hipótese nula e polo tanto considerar que os nosos coeficientes son distintos de 0. Análogamente tamén 10 1. Modelo de Regresión Loxística podemos empregar o método forward, neste caso en lugar de suprimir o término menos significativo, engadiremos aquela variable predictora, que ao introducila, o modelo resultante sexa máis significativo, é dicir, teña p-valores máis pequenos asociados as súas variables. A pesar de que os algoritmos descritos son de fácil implementación, non son os mellores para identificar o modelo desexado. Se en lugar de observar a significación de cada coeficiente queremos construír unha medida global para o modelo, podemos empregar un criterio moi popular, o criterio de información de Akaike, AIC. Este criterio para un modelo con verosimilitude lqe número de parámetros q defínese por: AIC =−2 log(lq)+2q. (1.25) Outra opción perfectamente válida para a mesma labor é o emprego do criterio de información bayesiano ou BIC. Este último está estreitamente ligado co criterio anterior e vén dado por: BIC =−2 log (lq) + qlog (n),(1.26) onde né o tamaño da mostra. Explorando modelos con distinto número de parámetros (de variables predictoras), finalmente seleccionaremos aquel modelo cun valor máis pequeno tanto de AIC como de BIC. Xa que se o valor do AIC ou BIC é pequeno significa que contamos cun modelo con gran verosimilitude e poucos parámetros. A idea principal é atopar aquel modelo que incorpore variables realmente útiles para así incrementar a verosimilitude. 1.3.1. Detección de datos atípicos Plantexamonos facer unha diagnose sobre o modelo de regresión loxística para encontrar algún punto inusual. Como no caso dos modelos lineais, os residuos son o máis importante para determinar que tan bos axustes temos para o modelo e onde sería aconsellable unha modificación ou mellora. Podemos calcular os residuos como a diferenza entre os valores observados e axustados. No caso do modelo lineal os residuos presentaban a mesma varianza e a suma residual de cadrados Pn i=1 ˆϵi2, era de gran utilidade para a diagnose do modelo. Pola contra, no modelo loxístico os residuos brutos non teñen a mesma varianza e como cantidade equivalente a suma residual de cadrados pódese empregar a deviance residual, Pn i=1 r2 i. Debido a isto, podemos recorrer aos residuos estandarizados, estos son o valor do residuo dividido entre unha estimación da súa desviación estándar. No caso do modelo de regresión loxística os residuos estandarizados son da forma: ri=yi−ˆpi pˆpi(1 −ˆpi).(1.27) Así, residuos estandarizados demasiado grandes en valor absoluto indican que a observación correspondente pode ser anómala ou atípica. Para o caso dos modelos lineais pódese considerar 1.4. Conclusión 11 como critero para que un dato sexa atípico aquél cun residuo estandarizado maior que 2 ou menor que -2. Ademais, se o modelo axustado é acertado, estos residuos seguirán unha distribución normal estándar. No caso da regresión loxística, a expresión para a deviance residual, Pn i=1 r2 i, ten a mesma distribución asintótica que a deviance, é dicir, unha chi-cadrado cos mesmos graos de liberdade. A partir dos residuos, tamén podemos elaborar gráficos en R que contrasten os residuos calculados fronte as predicións e interpretalos de forma similar aos modelos lineais. Pondo especial atención na detección de puntos inusuais, igual que para os modelos lineais, podemos examinar os leverages ou apalancamentos. Estos calcúlanse a partir da matriz H= ˆ V1/2X(X′ˆ V X)−1X′ˆ V1/2, onde ˆ Vé a matriz obtida a partir de V, sustituíndo os parámetros polos seus estimadores. O valor do leverage para a observación i-ésima vén dado polo elemento i-ésimo da diagonal principal da matriz H. Canto máis grande sexa o valor do leverage máis pequena será a varianza do residuo. Esta varianza pequena interprétase negativamente, pois as observacións asociadas poderán convertirse en observacións demasiado influíntes no modelo. 1.4. Conclusión Nun modelo de regresión loxística podemos tratar as mesmas cuestións que para o caso dos modelos lineais. Algunhas destas cuestións resólvense dun xeito moi similar en ambos modelos, mentres que noutros aspectos para a o modelo loxístico surxen máis dificultades. Dado que este último ten unha variable resposta discreta e a interpretación das variables predictoras obtidas pode resultar máis complexa. Deste xeito, a diferenza principal entre ambos modelos é que no caso do modelo de regresión loxística a variable resposta segue unha distribución discreta en concreto unha distribución de Bernoulli en lugar dunha distribución continua como é a normal para o caso lineal. Isto refléxase na saída prevista para a esperanza da resposta, que no caso do modelo loxístico débese atopar no intervalo [0,1]. 18 2. Modelo de Poisson con outra condición H1. Neste caso partindo da función de verosimilitude (2.6) e escribindoa en términos das β, empregando ηi=x′ iβtemos: α(β|(x1, y1)... (xn, yn)) = n Y i=1 exp(−exp(x′ iβ)) exp(x′ iβyi) yi!.(2.19) De novo como para o modelo loxístico, dados dous modelos con distinto número de variables explicativas, un con lvariables explicativas e outro con s, sendo s>l, e funcións de verosimilitude αleαsrespectivamente. Podemos de novo suxerir o seguinte estadístico para comparar os dous modelos aniñados: −2 log αl αs , o cal se poderá aproximar por unha ji-cadrado con s−lgraos de liberdade. A partir deste estadístico, xorde de novo o concepto de deviance definido en (1.22). Empregando a expresión para a verosimilitude (2.19) xunto coa definición (1.22) e tendo en conta que a verosimilitude do modelo saturado é 1, podemos expresar a deviance da regresión de Poisson como segue: D= 2 n X i=1 yilog yi ˆµi−(yi−ˆµi),(2.20) empregando ademais que ˆµi= exp(x′ iˆ β)e onde ˆ βé o estimador de máxima verosimilitude de β. Un valor elevado deste estadístico pode indicar un axuste pobre para o modelo. Mediante a diferenza das deviance e comparandoas con unha distribución χ2con tantos graos de liberdade como a diferenza entre o número de parámetros dos dous modelos, podemos comparar tamén dous modelos aniñados. Bondade de axuste Agora podemos testar a bondade de axuste do modelo proposto comparando a deviance do modelo fronte unha distribución χ2con tantos grados de liberdade como presente o noso modelo (McCullagh e Nelder, 1989). A maiores tamén se pode testar a significación individual dos predictores e construir intervalos de confianza para β, usando o erro estándar. Unha alternativa aχ2é efectuar un contraste. No cal a hipótese nula, H0, é que o modelo imposto é o correcto, fronte a hipótese alternativa, H1, que non o sexa. Desta forma, compre calcular un estadístico de contraste a partir da nosa mostra, neste caso o estadístico de Pearson X2: X2= n X i=1 (yi−ˆµi)2 ˆµi ,(2.21) que dado un nivel de significación α, imposto previamente, este estadístico (2.21) permitiranos efectuar o contraste. Dado que, baixo certas condicións de regularidade, a distribución de (2.21) 2.3. Inferencia sobre os parámetros do modelo 19 pódese aproximar por unha ji-cadrado con tantos graos de liberdade como a diferenza entre o número de observacións e a cantidade de parámetros do modelo e tendo isto en conta, poderemos obter un p-valor asociado a este estadístico. Co cal concluiremos o noso contraste, se o p-valor obtido é menor que o criterio de significación αpreestablecido, rexeitaremos a hipótese nula; no caso contrario acaptarémola. Cando o número de observacións é o suficientemente grande, a deviance e o X2son equivalentes. Unha alternativa, para testar o axuste do modelo é a interpretación do pseudo-R2de McFadden. Este vén dado pola seguinte expresión D2=NullDeviance−ResidualDeviance NullDeviance , onde ResidualDeviance é a diferenza entre a deviance do modelo que non depende de ningunha variable menos a do modelo que inclúe as variables explicativas, mentres que a NullDeviance é a deviance para un modelo que non depende de ningunha variable. Os valores deste coeficiente, D2, comprendidos entre 0,2e0,4indican según McFadden, un bo axuste para o modelo (Hensher e Stopher, 2021). 2.3.2. Intervalos de confianza para os parámetros de regresión Os intervalos de confianza para o modelo de Poisson pódense escribir de xeito análogo ao modelo de regresión loxística do Capítulo 1. Polo que consideraremos o intervalo de confianza de Wald (Cepeda-Cuervo et al., 2008) para o noso parámetro βj. Ademais , en virtude do teorema do límite central podemos deducir que a distribución de ˆ βconverxe a N(β,(X′V X)−1)(Cameron e Trivedi, 1998), sendo Va matriz definida en (2.13). Con isto, o intervalo de confianza de Wald para βjvén dado como: ˆ βj±zα/2ET (ˆ βj),(2.22) onde zα/2é o cuantil (1 −α/2) dunha distribución normal estándar. Véxase no Capítulo 1, que simplemente queda obter a fórmula para o erro típico de ˆ βj, como indicamos. Denotamos por ˆ Ja matriz ˆ J=X′ˆ V X, onde ˆ Vé a matriz obtida a partir da expresión de Vsustituíndo os parámetros polas súas estimacións. Deste xeito, o erro típico asociado a un ˆ βjconcreto vén dado como segue: ET (ˆ βj) = hˆ J−1 jj i1/2.(2.23) Deste modo, para calquera estimador ˆ βjpódese obter o seu erro típico e escribir o intervalo de Wald do parámetro βj. Este intervalo ten, para valores grandes de n, un nivel de confianza aproximado de 100(1 −α) %. Sen embargo para modelos con poucas observacións os coeficientes xeralmente non se achegan normalidade asumida. 20 2. Modelo de Poisson 2.4. Selección do Modelo En primeiro lugar podemonos plantexar identificar aquelas variables predictoras que non son útiles para modelar a nosa resposta. Para esta labor, podemos acudir aos métodos back-ward e forward detallados no Capítulo 1, e finalmente obter o modelo cuxas variables predictoras sexan máis significativas. Co fin de detectar de forma global o conxunto de predicións que mellor modelen a nosa resposta tamén podemos recorrer ao criterio de información de Akaike, AIC, ou ao criterio de información bayesiano, BIC, descritos polas expresións (1.25) e (1.26) respectivamente. Desta forma ao comparar o axuste de modelos, os valores máis baixos tanto de AIC como de BIC indicarán mellores axustes para o noso modelo. Dado que un valor baixo do AIC ou BIC implica maior verosimilitude e menos parámetros. 2.4.1. Detección de datos atípicos Análogamente ao modelo regresión loxística e ao visto para modelos lineais, se nos plantexamos facer unha diagnose sobre o modelo de Poisson podemos empregar os residuos, que determinarán que tan bos axustes temos para o noso modelo e se sería necesario algunha modificación. No modelo de Poisson, os residuos non teñen a mesma varianza polo que unha forma de medir a diferenza entre os valores observados e axustados neste modelo é a deviance residual. Establecemos a deviance residual de tal forma que se teña Pn i=1 r2 i=deviance, onde para este modelo concreto, os residuos da deviance veñen dados por: ri=sign(yi−ˆµi) [2 (yilog (yi/ˆµi)−yi+ ˆµi)]1/2,(2.24) onde sign a función é a función que obtén o signo do que tomemos como entrada. Desta forma, se o valor absoluto destes residuos é demasiado alto pode ser debido a que a observación correspondente é atípica. No caso da regresión lineal, podíamos considerar como candidato a dato atípico aquél cun residuo da deviance superior a 2 ou inferior a -2. 2.5. Sobredispersión Como sabemos, unha das características principais da regresión de Poisson é que a media e a varianza coinciden. Sen embargo, cando axustamos un modelo de Poisson a un conxunto de datos determinado, pode darse que estos valores difiran entre si. Isto é o que se coñece como sobredispersión do modelo de Poisson e ten consecuencias cualitativas similares ao incumprimento da hipótese de homocedasticidade no modelo de regresión lineal (Hilbe, 2014). 2.6. Conclusión 21 No paquete AER do software de R (R Core Team, 2021) hai un test que nos permitirá contrastar a sobredispersión do modelo de Poisson. O comando necesario para este test é dispersiontest. Dito test, realiza un contraste onde a hipótese nula H0é que o modelo de Poisson é equidisperso, fronte a alternativa H1de sobredispersión ou subdispersión, esta última menos frecuente (Hilbe, 2014). Observando o p-valor asociado ao test podemos rexeitar a hipótese nula se este valor é pequeno ou no caso contrario aceptala. No caso en que teñamos sobredispersión para o modelo, pódese empregar como alternativa a regresión de Poisson a regresión binomial negativa (McCullagh e Nelder, 1989). 2.6. Conclusión O modelo de regresión de Poisson presenta diferenzas significativas respecto do modelo lineal estándar e o modelo de regresión loxística. A diferenza principal é que para o modelo de Poisson a nosa variable resposta segue unha distribución discreta en concreto unha distribución de Poisson en lugar dunha distribución continua como é a distribución normal para o caso lineal. Recordemos que para o modelo loxístico a distribución de probabilidade da nosa variable resposta tamén era discreta pero seguía unha distribución de Bernoulli. Outra diferenza significativa é a varianza do modelo, no caso do modelo lineal esta é constante, mentres que no caso do modelo de regresión loxística e de Poisson non o é, ademais para este último a varianza coincide coa media. 22 2. Modelo de Poisson Capítulo 3 Aplicación do Modelo de Regresión Loxística Chegados a este punto, e recorrendo aos conceptos teóricos explicados no Capítulo 1, aplicaremos o modelo de regresión loxística a un conxunto de datos predeterminado. É dicir, analizaremos un conxunto de datos sobre o cancro de mama en Wisconsin, extraídos do repositorio de aprendizaxe automático UCI que en inglés é coñecido como UCI Machine Learning Repository (Dua e Graff, 2017). 3.1. Diagnóstico do cancro de mama 3.1.1. Descrición dos datos Dispoñemos dun conxunto de datos relativos ao diagnóstico do cancro de mama en Wisconsin. Este conxunto de datos foi doado por Nick Street en 1995 e os seus creadores son: Dr. William H. Wolberg, do departamento de Ciruxía Xeral, W. Nick Street e Olvi L. Mangasarian, do departamento de Ciencias da Computación, todos eles pertencentes a universidade de Wisconsin (Dua e Graff, 2017). A base de datos recolle un rexistro de 569 casos clínicos particulares, fichados cada un deles polo número de identificación do paciente. O noso obxetivo de traballo será partindo dunha serie de características, determinar se un tumor é maligno (M) ou benigno (B). As características das que dispomos, calculáronse a partir dunha imaxe dixitalizada dunha aspiración con agulla fina (FNA) de cada masa mamaria e describen principalmente particularidades dos núcleos celulares presentes na imaxe. 23 24 3. Aplicación do Modelo de Regresión Loxística Centraremonos nas variables indicadas a continuación. Dentro deste conxunto de variables, a variable resposta do noso modelo será diagnosis, mentres que as demais serán as variables explicativas que pretenden modelar a resposta: Media do radio (radius_mean): media das distancias dende o centro do núcleo celular aos puntos do perímetro do mesmo. Tomará valores como: 18, 20.6, 19.7, 11.4, 20.3, ... Media da textura (texture_mean): variable de valores númericos (como 10.4, 17.8, 21.2, 20.4, 14.3, ...) referente a media das desviacións estándar dos valores da escala de grises. Media do perímetro (perimeter_mean): media do perímetro dos núcleos celulares, os valores posibles para esta variable poden ser 122.8, 132.9, 130, 77.6, 135.1, ... Media da área (area_mean): media da área dos núcleos celulares. Podemos observar os seguintes valores asociados a esta variable: 1001, 1326, 1203, 386, 1297, ... Media da uniformidade (smoothness_mean): media da uniformidade, é dicir, media das variacións locais das lonxitudes dos radios, presenta valores como: 0.1184, 0.0847, 0.1096, 0.1425, 0.1003, ... Media da compacidade (compactness_mean): media da compacidade que se calcula como perimetro2 area −1. Esto da lugar aos seguintes posibles valores para dita variable: 0.2776, 0.0786, 0.1599, 0.2839, 0.1328, ... Media da concavidade (concavity_mean): media da concavidade, é dicir, severidade das porcións cóncavas do contorno. Contamos cos seguintes valores asociados a variable: 0.3001, 0.0869, 0.1974, 0.2414, 0.198, ... Media dos puntos cóncavos (concave.points_mean): media dos puntos cóncavos, é dicir, número de porcións cóncavas do contorno. Pode tomar os seguintes valores: 0.1471, 0.0702, 0.1279, 0.1052, 0.1043, ... Media da simetría (symmetry_mean): media das simetrías. Variable númerica con valores como: 0.242, 0.181, 0.207, 0.26, 0.181, ... Media da dimensión do fractal (fractal_dimension_mean): media das dimensións dos fractais. A variable pode tomar os seguintes valores: 0.0787, 0.0567, 0.06, 0.0974, 0.0588, ... Dado este conxunto de variables que describen características asociadas a cada masa mamaria, o obxetivo é propor o mellor modelo posible para o diagnóstico dun tumor. Polo tanto, podemos plantexar un modelo de regresión loxística que nos permita modelar a nosa resposta dicotómica, de se certo tumor resulta ser maligno ou benigno. 3.1. Diagnóstico do cancro de mama 25 3.1.2. Aplicación do modelo de regresión loxística. Escritura e selección do modelo. A regresión loxística permite estimar a probabilidade dunha variable cualitativa binaria en función de certas variables cuantitativas. A saída prevista para a esperanza da nosa resposta atoparase dentro do intervalo [0,1], mentres que se empregamos un axuste lineal para modelar a nosa resposta é posible que os valores pronosticados para a nosa resposta non se encadren en tal intervalo. Observando a variable resposta, diagnosis, vemos que do total de diagnósticos a maioría resultarán ser tumores malignos, 357, fronte a 212 que serán benignos. Podemos aplicar aos nosos datos un modelo de regresión loxística, para ver que variables explicativas das descritas inflúen no diagnóstico do cancro. Para iso comezaremos escribindo o modelo formado pola resposta diagnosis e todas as variables explicativas indicadas e faremos un summary do mesmo. Call: glm(formula = diagnosis ~ radius_mean + texture_mean + perimeter_mean + area_mean + smoothness_mean + compactness_mean + concavity_mean + concave.points_mean + symmetry_mean + fractal_dimension_mean, family = binomial, data = datos) Deviance Residuals: Min 1Q Median 3Q Max -1.95590 -0.14839 -0.03943 0.00429 2.91690 Coefficients: Estimate Std. Error z value Pr(>|z|) (Intercept) -7.35952 12.85259 -0.573 0.5669 radius_mean -2.04930 3.71588 -0.551 0.5813 texture_mean 0.38473 0.06454 5.961 2.5e-09 *** perimeter_mean -0.07151 0.50516 -0.142 0.8874 area_mean 0.03980 0.01674 2.377 0.0174 * smoothness_mean 76.43227 31.95492 2.392 0.0168 * compactness_mean -1.46242 20.34249 -0.072 0.9427 concavity_mean 8.46870 8.12003 1.043 0.2970 concave.points_mean 66.82176 28.52910 2.342 0.0192 * symmetry_mean 16.27824 10.63059 1.531 0.1257 fractal_dimension_mean -68.33703 85.55666 -0.799 0.4244 26 3. Aplicación do Modelo de Regresión Loxística --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 (Dispersion parameter for binomial family taken to be 1) Null deviance: 751.44 on 568 degrees of freedom Residual deviance: 146.13 on 558 degrees of freedom AIC: 168.13 Number of Fisher Scoring iterations: 9 Comezaremos ollando os p-valores obtidos para cada unha das variables. O contraste que se plantexa ten como hipótese nula, H0, que o coeficiente asociado a certa variable sexa cero ou equivalentemente que certa variable non explica a nosa resposta, fronte a hipótese alternativa, H1, que indica o contrario. A interpretación destes p-valores é moi similar a do modelo lineal. É dicir, as variables de interese serán aquelas que teñan uns p-valores máis pequenos que un nivel dado, xa que os p-valores menores indicarán maior significación no modelo da variable asociada. En concreto, no caso que nos ocupa vemos que as únicas variables que teñen p-valores baixos, é dicir, son significativas para un nivel de significación α= 0,05, son texture_mean, area_mean,smoothness_mean econcave.points_mean. Con isto, descubrimos que non todas as variables preditoras indicadas son de utilidade para explicar a nosa resposta, polo que tentaremos seleccionar un subconxunto de variables que modele o mellor posible a resposta. Para iso, empregaremos o método backware descrito na Sección 1.3. A estratexia habitual é eliminar aquela variable predictora cuxo p-valor asociado tome o valor máis elevado. Volvendo a ollar o summary, vemos que a variable menos significativa é compactness_mean (cun p-valor asociado igual a 0.9427). Se volvemos axustar o modelo de regresión loxística prescindindo desta variable explicativa, obteríamos que perimeter_mean é a variable menos significativa, cun p-valor asociado igual a 0.7750. Repetindo iterativamente este proceso ata que os p-valores asociados a todas as variables aleatorias sexan menores que 0.05, vanse eliminando sucesivamente as variables concavity_mean,fractal_dimension_mean esymmetry_mean. Obtendo como resultado, un modelo final no que as únicas variables explicativas serían: radius_mean,texture_mean,area_mean, smoothness_mean econcave.points_mean. Neste último modelo plantexado, todas as variables predictoras son significativas polo que pomos fin o algoritmo do método backware. Sen embargo, podemos observar que non temos significación para o intercepto xa que o p-valor asociado para este é igual a 0.55983. Como explicamos no Capítulo 1, na Sección 1.3, se en lugar de observar a significación de 3.1. Diagnóstico do cancro de mama 27 cada coeficiente o que queremos é construír unha medida global para a selección do modelo, podemos empregar o criterio de información de Akaike, AIC. Para isto, simplemente debemos comparar os valores de AIC asociados a cada un dos modelos plantexados e quedarnos con aquel que presente un valor de AIC menor, xa que ese será o modelo máis verosímil. O AIC asociado ao primeiro modelo é 168.13, mentres que os asociados aos seguintes son 166.14, 164.22, 163.28, 162.38 e 162.68 respectivamente. Tendo en conta estes valores, o modelo máis plausible acorde a dito criterio é o penúltimo, é dicir, aquel cuxas variables explicativas son: radius_mean, texture_mean,area_mean,smoothness_mean,concave.points_mean,symmetry_mean. Interpretación dos coeficientes. Partindo do modelo seleccionado anteriormente, é dicir, aquel cuxas variables explicativas son radius_mean,texture_mean,area_mean,smoothness_mean,concave.points_mean esymmetry_mean podemos por atención nos coeficientes obtidos. A interpretación destes coeficientes difire da do modelo lineal. Como explicamos no Capítulo 1, o modelo de regresión loxística emprega unha función enlace coñecida como función logit e dada pola expresión (1.2). En consecuencia, para levar a cabo a interpretación desexada debemos recorrer ao concepto de odds definido no Capítulo 1, que para o modelo loxístico é igual a odds = expβ0·expβ1X1·expβ2X2·... ·expβqXq. Desta forma, en lugar de interpretar directamente, ˆ β, o vector de coeficientes axustados asociados a cada variable, que no noso caso particular ten a seguinte expresión: ˆ β= (−8,6108,−2,7251,0,3852,0,0430,58,7854,73,7015,15,5621),(3.1) e onde a primeira entrada deste vector correspode ao coeficiente asociado ao intercepto, a segunda ao coeficiente asociado a variable radius_mean e así sucesivamente ata a última entrada correspondente a variable symmetry_mean, pode ser máis útil a interpretación dos valores resultantes de aplicarlle a expoñencial a cada coeficiente. Deste xeito, os coeficientes do modelo de regresión loxística interprétanse como o logarítmo das odds ratio. É dicir, agora β1pódese interpretar como segue: un aumento unitario en X1con X2, ... , Xqfixados aumenta o logaritmo das odds asociado a β1ou o que é o mesmo, aumenta as probabilidades de éxito dun factor exp β1. Se nos fixamos no coeficiente axustado asociado a variable texture_mean, 0.3852, este indica que o logaritmo das odds ratio de diagnosis, aumenta en 0.3852 unidades por cada unidade que aumenta texture_mean. Ou equivalentemente, o aumento dunha unidade da variable texture_mean, permanecendo as demais constantes, implicará que a probabilidade de ter un tumor maligno se multiplique por 1.469941. Por outra banda, o coeficiente relativo a variable radius_mean, -2.72, indícanos que o logaritmo das odds ratio de diagnosis diminúe en 2.72 unidades por cada unidade que aumenta radius_mean. Ou o que é o mesmo que por cada unidade que aumente radius_mean a probabilidade de ter un tumor maligno multiplícase por 0.065. Deste modo, a interpretación do valor de ˆ βjasociado a cada variable difire según o signo deste. Se o valor de ˆ βjé positivo indica 34 3. Aplicación do Modelo de Regresión Loxística Figura 3.2: Residuos do modelo con valor elevado. head(hat.valores.ma) 113 492 538 505 77 529 0.18483481 0.10553580 0.10463376 0.09305295 0.06958510 0.06948025 Polo tanto no caso que nos ocupa, obtemos en xeral valores pequenos para os leverages, o que indica que a varianza do residuo non é pequena. Esto é positivo, dado que unha varianza do residuo pequena pode dar lugar a observacións demasiado influíntes no modelo. 3.1.5. Conclusión. Partindo do modelo con mellor significación para as variables, é dicir aquel que inclúe como variables explicativas as seguintes: texture_mean,area_mean,smoothness_mean econcave.points_mean, será de gran interese comparar as predicións obtidas coas observacións. Desta maneira, veremos que tan bo resulta o noso modelo para predicir a nosa resposta. Para levar a cabo este análise, asumiremos que o modelo clasifica a un paciente cun tumor maligno se ˆpi≥0,5. 3.1. Diagnóstico do cancro de mama 35 Figura 3.3: Gráficos de residuos empregando o paquete ’DHARMa’. Ademais, empregaremos o programa R (R Core Team, 2021), para elaborar unha matriz que reflexe a compación mencionada. Ademais, para maior claridade, tamén podemos representar dita matriz mediante un mapa de calor, como mostra a Figura 3.4 . Deste xeito, os valores en verde, é dicir os da diagonal principal da matriz, correspóndense cos valores estimados correctamente polo modelo. Mentres que a outra diagonal, indicada na Figura 3.4 en vermello, representa os casos onde o modelo se equivoca. predicciones observaciones 0 1 0 343 14 1 20 192 Desta forma o modelo é capaz de clasificar correctamente 343+192 343+14+20+192 = 0,94024, aproximadamente un 94 % das observacións ou o que é o mesmo a tasa de erro do modelo é aproximadamente un 6 % (Faraway, 2016). A fracción de pacientes que se prevé que non sexan diagnosticados cun tumor benigno é 343 343+14 = 0,96078, o que resulta un valor elevado. Este valor é o que se coñece como especificidad do test (Faraway, 2016). En contraste, a proporción que será diagnosticada dun tumor maligno é o que se coñece como sensibilidade (Faraway, 2016), é será 192 20+192 = 0,90566. Polo tanto vemos que é probable que o noso proceso preditivo detecte a presenza dun tumor maligno a partir das variables predictoras das que dispoñemos. Unha forma alternativa para determinar se o noso modelo obtén as predicións correctas 36 3. Aplicación do Modelo de Regresión Loxística Figura 3.4: Matriz de confusión en forma de mapa de calor. consiste en elaborar un gráfico que represente a probabilidade dun tumor predita polo modelo, ˆpi, fronte aos valores observados. Se o modelo é acertado, este gráfico mostrará valores preto do (0,0) e do (1,1). Na Figura 3.5, podemos observar dita representación para o noso modelo en cuestión. Vemos que obtemos unha gran cantidade de valores preto do (0,0) e do (1,1), co que podemos concluir que o modelo indicado é correcto para predecir a nosa variable resposta. 3.1. Diagnóstico do cancro de mama 37 Figura 3.5: Probabilidades preditas fronte a valores observados polo modelo. 38 3. Aplicación do Modelo de Regresión Loxística Capítulo 4 Aplicación do Modelo de Poisson Neste capítulo, aplicaremos o modelo de Poisson a un conxunto de datos predeterminado. Para elo, empregaremos as cuestións teóricas explicadas no Capítulo 2, que nos permitirán analizar os nosos datos sobre a demanda de atención médica. 4.1. Demanda de atención médica 4.1.1. Descrición dos datos Partimos dun conxunto de datos presente no paquete ’AER’ do software de R (R Core Team, 2021) e en concreto procedentes da librería ’MASS’. A base de datos recolle información transversal procedente da enquisa de saúde australiana dos anos 1977 e 1978. En concreto, o marco de datos contén 5190 observacións sobre 12 variables. Para o noso estudo, tomaremos como variable resposta a variable visits, que fai referencia ao número de visitas ao doctor nas últimas dúas semanas. O obxetivo do noso modelo será modelar esta demanda de atención médica en función das seguintes variables explicativas: Idade (age): idade en anos dividida entre 100. Ingresos (income): ingresos anuais en decenas de miles de dólares. Con estas variables, que describen particularidades de cada paciente, o obxetivo é plantexar un modelo para explicar a demanda de atención médica en función das nosas variables explicativas. Dado que a variable resposta é discreta e ademais fai referencia a un número de feitos que ocorreron nun tempo determinado, podemos plantexar un modelo de regresión de Poisson para modelar dita resposta. 39 40 4. Aplicación do Modelo de Poisson 4.1.2. Aplicación do modelo de Poisson. Escritura e selección do modelo. Como indicamos no Capítulo 2, o modelo de Poisson é un dos modelos que permiten modelar datos de conteo. Para aplicarlle este modelo a nosa resposta, esta debe ser discreta e seguir unha distribución de Poisson. A propiedade principal desta distribución é que a media e a varianza coinciden. Ademais, ao contrario que para a regresión loxística, a resposta para a regresión de Poisson non se limita aos valores {0,1}, senón que pode tomar calquer valor enteiro positivo. Observando a variable dependente, visits, vemos que segue aparentemente unha distribución de Poisson como podemos observar na Figura 4.1, cuxa forma se asemella bastante as primeiras gráficas presentes na Figura 2.1. Desta forma, podemos plantexar un modelo de Poisson que Figura 4.1: Gráfico de barras coa frecuencia absoluta da demanda de atención médica. pretenda explicar o número de visitas ao consultorio médico en función das variables explicativas indicadas. Polo que escribimos o modelo formado pola resposta visits e as variables explicativas age eincome e facemos un summary do mesmo. Call: glm(formula = visits ~ age + income, family = poisson, data = DoctorVisits) 4.1. Demanda de atención médica 41 Deviance Residuals: Min 1Q Median 3Q Max -1.0203 -0.7986 -0.6691 -0.6176 6.3802 Coefficients: Estimate Std. Error z value Pr(>|z|) (Intercept) -1.60661 0.08596 -18.690 < 2e-16 *** age 1.35292 0.12666 10.681 < 2e-16 *** income -0.34165 0.07770 -4.397 1.1e-05 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 (Dispersion parameter for poisson family taken to be 1) Null deviance: 5634.8 on 5189 degrees of freedom Residual deviance: 5449.9 on 5187 degrees of freedom AIC: 7787.5 Number of Fisher Scoring iterations: 6 Comezamos igual que no Capítulo 3, ollando os p-valores obtidos para cada unha das nosas variables. De novo, a interpretación destes valores é bastante semellante a do modelo lineal e a do modelo de regresión loxística. Se este p-valor é baixo indicaranos que a variable asociada é significativa para modelar a resposta. Neste caso particular, vemos que todas as nosas variables resultan ser significativas, para un nivel de significación α= 0,05. Polo que todas elas son de utilidade para explicar a demanda de atención médica. Se en lugar de observar a significación de cada coeficiente o que queremos é construír unha medida global para a selección do modelo, podemos empregar o criterio de información de Akaike, AIC. Para isto, simplemente debemos comparar os valores de AIC asociados a cada un dos modelos plantexados e quedarnos con aquel que presente un valor de AIC menor, xa que ese será o modelo máis verosímil. Para isto podemos ter en conta o modelo plantexado anteriormente e composto pola variable resposta visitas e as varibles explicativas idade e ingresos, outro modelo formado só pola variable resposta e como explicativa a variable idade e un último modelo composto pola resposta e a variable explicativa ingresos. Call: glm(formula = visits ~ age, family = poisson, data = DoctorVisits) 42 4. Aplicación do Modelo de Poisson Deviance Residuals: Min 1Q Median 3Q Max -0.9650 -0.8266 -0.6552 -0.6402 6.5608 Coefficients: Estimate Std. Error z value Pr(>|z|) (Intercept) -1.87944 0.06248 -30.08 <2e-16 *** age 1.54878 0.12060 12.84 <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: 5634.8 on 5189 degrees of freedom Residual deviance: 5470.0 on 5188 degrees of freedom AIC: 7805.6 Number of Fisher Scoring iterations: 6 Call: glm(formula = visits ~ income, family = poisson, data = DoctorVisits) Deviance Residuals: Min 1Q Median 3Q Max -0.9147 -0.8235 -0.7526 -0.6193 6.7789 Coefficients: Estimate Std. Error z value Pr(>|z|) (Intercept) -0.87152 0.04567 -19.082 < 2e-16 *** income -0.59991 0.07486 -8.014 1.11e-15 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 (Dispersion parameter for poisson family taken to be 1) Null deviance: 5634.8 on 5189 degrees of freedom Residual deviance: 5566.5 on 5188 degrees of freedom 4.1. Demanda de atención médica 43 AIC: 7902 Number of Fisher Scoring iterations: 6 Agora, simplemente ollando no summary de cada modelo o valor do AIC obtido, temos que o modelo máis plausible tendo en conta o criterio do AIC é o composto pola resposta e ambas variables explicativas, idade e ingresos. Dado que o AIC asociado a este modelo é 7787,5, mentres que o asociado aos outros modelos indicados é 7805,6e7902 respectivamente. Interpretación dos coeficientes. Tal e como explicamos na Sección 2.1, o modelo de Poisson emprega unha función enlace logarítmica. Deste xeito poderemos expresar a media da nosa resposta como: log(µi) = ηi=x′ iβ ou equivalentemente µi= exp(x′ iβ). Como consecuencia destas expresións descritas, poderemos interpretar os nosos parámetros da forma: se certa variable explicativa xjaumenta en nunidades e as demais permanecen fixas, a media para a variable de Poisson multiplicase pola potencia nésima de exp(βj). Para o caso que nos ocupa, a estimación do vector de coeficientes asociado ao noso modelo vén dado pola seguinte expresión: ˆ β= (−1,60,1,35,−0,34),(4.1) onde a primeira entrada do vector ˆ β, corresponde ao coeficiente asociado ao intercepto e as demais aos coeficientes asociados as respectivas variables explicativas, age eincome. Agora, partindo dos coeficientes obtidos para o noso modelo (4.1) e tendo en conta a función enlace empregada sábese que a interpretación da expoñencial destos coeficientes resulta máis interesante. exp ˆ β= (0,20,3,86,0,71),(4.2) con isto, podemos interpretar que o aumento dunha unidade da variable idade, permanecendo as restantes fixas, implicará que o número medio de visitas ao consultorio se multiplique por 3,86. Por outra banda, no caso da variable income, referente ao ingreso familiar anual, o seu aumento diminuirá o número de visitas ao consultorio. Máis concretamente, un aumento unitario da variable ingresos, mantendo as demais variables explicativas constantes, significa que o número de visitas ao consultorio médico se multiplican por 0,71. En conclusión, teremos que se o valor da expoñencial do coeficiente é maior que 1, dita variable aumentará o valor da media da nosa resposta, como é o caso da variable age, mentres que se o valor da expoñencial do coeficiente é menor que 1, como no caso da variable income, diminuirá o valor da mesma. No caso de que fose xustamente 1, a variable non terá ningunha influencia para o modelaxe da nosa resposta. 50 4. Aplicación do Modelo de Poisson Bibliografía Agresti, A. G. (1990). Categorical data analysis. New York: John Wiley and Sons. Cameron, A. C., e Trivedi, P. K. (1998). Regression analysis of count data. Cambridge: Cambridge University Press. Cepeda-Cuervo, E., Aguilar, W., Cervantes, V., Corrales, M., Díaz, I., e Rodríguez, D. (2008). Intervalos de confianza e intervalos de credibilidad para una proporción. Revista Colombiana de Estadística,31(2), 211-228. Colaboradores de Wikipedia. (2022). Información de Fisher — Wikipedia, la enciclopedia libre. https://en.wikipedia.org/w/index.php?title=Fisher_information&oldid= 1079985740. ([En línea; consultado el 2 de mayo de 2022]) Dua, D., e Graff, C. (2017). UCI Machine Learning Repository. http://archive.ics.uci.edu/ ml. ([En línea; consultado el 12 de mayo de 2022]) Faraway, J. J. (2016). Extending the linear model with R: Generalized Linear, Mixed Effects and Nonparametric Regression Models. Boca Raton: CRC Press. Fox, J. (2008). Applied regression analysis and generalized linear models (2nd ed.). Los Angeles: Sage. Hartig, F. (2022). DHARMa: Residual Diagnostics for Hierarchical (Multi-Level / Mixed) Regression Models [Manual]. Obtido de:https://CRAN.R-project.org/package=DHARMa (R package version 0.4.5) Hastie, T. J., Tibshirani, R. J., e Friedman, J. (2009). The elements of statistical learning: data mining, inference, and prediction (2nd ed ed.). New York: Springer. Hensher, D., e Stopher, P. (2021). Behavioural travel modelling. Taylor & Francis. Obtido de:https://books.google.com.do/books?id=Z_UlEAAAQBAJ Hilbe, J. M. (2014). Modeling count data. New York: Cambridge University Press. McCullagh, P., e Nelder, J. A. (1989). Generalized linear models (2nd ed.). London: Chapman and Hall. R Core Team. (2021). R: A Language and Environment for Statistical Computing [Manual]. Vienna, Austria. Obtido de:https://www.R-project.org/ Rivera, J. I. y. G. (n.d.). Chapter 8 regresión de poisson | modelos lineales generaliza51 52 Bibliografía dos con r. Obtido de:https://bookdown.org/jaimeisaacp/bookglm/regresi%C3%B3n-de -poisson.html Venables, W., e Ripley, B. (2002). Modern applied statistics with S (4th ed.). Berlin: Springer. Anexos 53 55 O análise de datos desenrolado ao longo dos Capítulos 3 e 4, realizouse mediante o software estatístico de R. Para maior precisión, adxunto o código empregado para replicar todos os resultados expostos. #Lectura dos datos library(readxl) library(dplyr) ## ## Attaching package: ’dplyr’ ## The following objects are masked from ’package:stats’: ## ## filter, lag ## The following objects are masked from ’package:base’: ## ## intersect, setdiff, setequal, union datos <- read.csv("data.csv") attach(datos) #Recodificación da variable diagnosis levels(datos$diagnosis) <- c("0","1") levels(datos$diagnosis) ## [1] "0" "1" datos$diagnosis<- factor(datos$diagnosis) table(datos$diagnosis) ## ## 0 1 ## 357 212 #Escribimos o modelo modelo<-glm(diagnosis~radius_mean+texture_mean+perimeter_mean+area_mean+ smoothness_mean+compactness_mean+concavity_mean+concave.points_mean+ 56 symmetry_mean+fractal_dimension_mean,data=datos,family=binomial) summary(modelo) ## ## Call: ## glm(formula = diagnosis ~ radius_mean + texture_mean + perimeter_mean + ## area_mean + smoothness_mean + compactness_mean + concavity_mean + ## concave.points_mean + symmetry_mean + fractal_dimension_mean, ## family = binomial, data = datos) ## ## Deviance Residuals: ## Min 1Q Median 3Q Max ## -1.95590 -0.14839 -0.03943 0.00429 2.91690 ## ## Coefficients: ## Estimate Std. Error z value Pr(>|z|) ## (Intercept) -7.35952 12.85259 -0.573 0.5669 ## radius_mean -2.04930 3.71588 -0.551 0.5813 ## texture_mean 0.38473 0.06454 5.961 2.5e-09 *** ## perimeter_mean -0.07151 0.50516 -0.142 0.8874 ## area_mean 0.03980 0.01674 2.377 0.0174 * ## smoothness_mean 76.43227 31.95492 2.392 0.0168 * ## compactness_mean -1.46242 20.34249 -0.072 0.9427 ## concavity_mean 8.46870 8.12003 1.043 0.2970 ## concave.points_mean 66.82176 28.52910 2.342 0.0192 * ## symmetry_mean 16.27824 10.63059 1.531 0.1257 ## fractal_dimension_mean -68.33703 85.55666 -0.799 0.4244 ## --- ## Signif. codes: 0 '***'0.001 '**'0.01 '*'0.05 '.'0.1 ' ' 1 ## ## (Dispersion parameter for binomial family taken to be 1) ## ## Null deviance: 751.44 on 568 degrees of freedom ## Residual deviance: 146.13 on 558 degrees of freedom ## AIC: 168.13 ## ## Number of Fisher Scoring iterations: 9 57 #Selección do modelo (método backware) modelo1<-glm(diagnosis~radius_mean+texture_mean+perimeter_mean+ area_mean+smoothness_mean+concavity_mean+ concave.points_mean+symmetry_mean+fractal_dimension_mean, data=datos,family=binomial) summary(modelo1) ## ## Call: ## glm(formula = diagnosis ~ radius_mean + texture_mean + perimeter_mean + ## area_mean + smoothness_mean + concavity_mean + concave.points_mean + ## symmetry_mean + fractal_dimension_mean, family = binomial, ## data = datos) ## ## Deviance Residuals: ## Min 1Q Median 3Q Max ## -1.97463 -0.14756 -0.03947 0.00428 2.92261 ## ## Coefficients: ## Estimate Std. Error z value Pr(>|z|) ## (Intercept) -7.04893 12.08068 -0.583 0.5596 ## radius_mean -1.89732 3.05194 -0.622 0.5342 ## texture_mean 0.38463 0.06455 5.959 2.54e-09 *** ## perimeter_mean -0.09813 0.34325 -0.286 0.7750 ## area_mean 0.03999 0.01649 2.424 0.0153 * ## smoothness_mean 76.08009 31.50588 2.415 0.0157 * ## concavity_mean 8.48311 8.10558 1.047 0.2953 ## concave.points_mean 66.77285 28.49046 2.344 0.0191 * ## symmetry_mean 16.16553 10.51539 1.537 0.1242 ## fractal_dimension_mean -71.99700 68.75680 -1.047 0.2950 ## --- ## Signif. codes: 0 '***'0.001 '**'0.01 '*'0.05 '.'0.1 ' ' 1 ## ## (Dispersion parameter for binomial family taken to be 1) ## ## Null deviance: 751.44 on 568 degrees of freedom ## Residual deviance: 146.14 on 559 degrees of freedom 58 ## AIC: 166.14 ## ## Number of Fisher Scoring iterations: 9 modelo2<-glm(diagnosis~radius_mean+texture_mean+area_mean+ smoothness_mean+concavity_mean+concave.points_mean+symmetry_mean+ fractal_dimension_mean,data=datos,family=binomial) summary(modelo2) ## ## Call: ## glm(formula = diagnosis ~ radius_mean + texture_mean + area_mean + ## smoothness_mean + concavity_mean + concave.points_mean + ## symmetry_mean + fractal_dimension_mean, family = binomial, ## data = datos) ## ## Deviance Residuals: ## Min 1Q Median 3Q Max ## -1.96847 -0.15195 -0.04024 0.00409 2.93549 ## ## Coefficients: ## Estimate Std. Error z value Pr(>|z|) ## (Intercept) -5.27847 10.31074 -0.512 0.60869 ## radius_mean -2.68473 1.32326 -2.029 0.04247 * ## texture_mean 0.38262 0.06413 5.966 2.42e-09 *** ## area_mean 0.04157 0.01554 2.675 0.00747 ** ## smoothness_mean 78.22119 30.57445 2.558 0.01052 * ## concavity_mean 8.25689 8.04476 1.026 0.30472 ## concave.points_mean 64.07659 26.75842 2.395 0.01664 * ## symmetry_mean 16.02120 10.47671 1.529 0.12621 ## fractal_dimension_mean -82.21451 58.85970 -1.397 0.16248 ## --- ## Signif. codes: 0 '***'0.001 '**'0.01 '*'0.05 '.'0.1 ' ' 1 ## ## (Dispersion parameter for binomial family taken to be 1) ## ## Null deviance: 751.44 on 568 degrees of freedom ## Residual deviance: 146.22 on 560 degrees of freedom 59 ## AIC: 164.22 ## ## Number of Fisher Scoring iterations: 9 modelo3<-glm(diagnosis ~radius_mean+texture_mean+area_mean+ smoothness_mean+concave.points_mean+ symmetry_mean+fractal_dimension_mean, data=datos,family=binomial) summary(modelo3) ## ## Call: ## glm(formula = diagnosis ~ radius_mean + texture_mean + area_mean + ## smoothness_mean + concave.points_mean + symmetry_mean + fractal_dimension_mean, ## family = binomial, data = datos) ## ## Deviance Residuals: ## Min 1Q Median 3Q Max ## -2.00824 -0.15647 -0.04323 0.00326 2.83457 ## ## Coefficients: ## Estimate Std. Error z value Pr(>|z|) ## (Intercept) -3.14334 9.86273 -0.319 0.74995 ## radius_mean -3.06595 1.24925 -2.454 0.01412 * ## texture_mean 0.38122 0.06374 5.981 2.22e-09 *** ## area_mean 0.04562 0.01485 3.072 0.00212 ** ## smoothness_mean 60.39599 23.93789 2.523 0.01163 * ## concave.points_mean 83.46572 19.16188 4.356 1.33e-05 *** ## symmetry_mean 17.04684 10.23216 1.666 0.09571 . ## fractal_dimension_mean -49.23528 47.54527 -1.036 0.30041 ## --- ## Signif. codes: 0 '***'0.001 '**'0.01 '*'0.05 '.'0.1 ' ' 1 ## ## (Dispersion parameter for binomial family taken to be 1) ## ## Null deviance: 751.44 on 568 degrees of freedom ## Residual deviance: 147.28 on 561 degrees of freedom ## AIC: 163.28 66 1-pchisq(G2, df =4) ## [1] 0.6889965 #Contraste empregando función anova library(car) ## Loading required package: carData ## ## Attaching package: ’car’ ## The following object is masked from ’package:dplyr’: ## ## recode Anova(modelo) ## Analysis of Deviance Table (Type II tests) ## ## Response: diagnosis ## LR Chisq Df Pr(>Chisq) ## radius_mean 0.306 1 0.57992 ## texture_mean 49.209 1 2.3e-12 *** ## perimeter_mean 0.020 1 0.88749 ## area_mean 5.496 1 0.01906 * ## smoothness_mean 6.289 1 0.01215 * ## compactness_mean 0.005 1 0.94270 ## concavity_mean 1.103 1 0.29370 ## concave.points_mean 5.795 1 0.01607 * ## symmetry_mean 2.307 1 0.12878 ## fractal_dimension_mean 0.647 1 0.42136 ## --- ## Signif. codes: ## 0 '***'0.001 '**'0.01 '*'0.05 '.'0.1 ' ' 1 modelo6<-glm(diagnosis~texture_mean+area_mean+smoothness_mean +concave.points_mean,data=datos,family = binomial) summary(modelo6) 67 ## ## Call: ## glm(formula = diagnosis ~ texture_mean + area_mean + smoothness_mean + ## concave.points_mean, family = binomial, data = datos) ## ## Deviance Residuals: ## Min 1Q Median 3Q Max ## -2.31798 -0.15623 -0.04212 0.01662 2.84201 ## ## Coefficients: ## Estimate Std. Error z value Pr(>|z|) ## (Intercept) -23.677816 3.882774 -6.098 1.07e-09 ## texture_mean 0.362687 0.060544 5.990 2.09e-09 ## area_mean 0.010342 0.002002 5.165 2.40e-07 ## smoothness_mean 59.471304 25.965153 2.290 0.022 ## concave.points_mean 76.571210 16.427864 4.661 3.15e-06 ## ## (Intercept) *** ## texture_mean *** ## area_mean *** ## smoothness_mean * ## concave.points_mean *** ## --- ## Signif. codes: ## 0 '***'0.001 '**'0.01 '*'0.05 '.'0.1 ' ' 1 ## ## (Dispersion parameter for binomial family taken to be 1) ## ## Null deviance: 751.44 on 568 degrees of freedom ## Residual deviance: 156.44 on 564 degrees of freedom ## AIC: 166.44 ## ## Number of Fisher Scoring iterations: 8 #Intervalos de confianza para os parámetros do modelo confint.default(modelo) ## 2.5 % 97.5 % 68 ## (Intercept) -32.55012751 17.83109229 ## radius_mean -9.33229662 5.23368682 ## texture_mean 0.25824448 0.51122420 ## perimeter_mean -1.06161528 0.91859445 ## area_mean 0.00698718 0.07260522 ## smoothness_mean 13.80178675 139.06276076 ## compactness_mean -41.33297960 38.40813510 ## concavity_mean -7.44627507 24.38367459 ## concave.points_mean 10.90575011 122.73776358 ## symmetry_mean -4.55732264 37.11380728 ## fractal_dimension_mean -236.02499401 99.35094023 #Diagnose sobre o modelo #Cálculo de residuos de Pearson res <- residuals(modelo6, type ="pearson") head(res) #residuos de Pearson dos 6 primeiros pacientes ##12345 ## 0.012625601 0.031857290 0.001676578 0.120231091 0.011750045 ## 6 ## 0.689101952 res.sig <- abs(res) >2#Residuos de Pearson significativos table(res.sig) ## res.sig ## FALSE TRUE ## 559 10 res.orde <- sort(abs(res[res.sig]), decreasing =TRUE)#mostramos os máis altos head(res.orde) ## 298 41 136 172 129 74 ## 7.466056 6.108532 4.130008 3.852646 3.698711 2.773438 #Probabilidades de diagnose de tumor maligno asociadas fitted.values(modelo6)[298] 69 ## 298 ## 0.01762363 fitted.values(modelo6)[41] ## 41 ## 0.02610001 fitted.values(modelo6)[136] ## 136 ## 0.05538029 fitted.values(modelo6)[172] ## 172 ## 0.06311982 fitted.values(modelo6)[129] ## 129 ## 0.9318822 fitted.values(modelo6)[74] ## 74 ## 0.1150489 #Cálculo de residuos estandarizados res.std <- rstandard(modelo6, type ="pearson") res.std.sig <- abs(res.std) >2#Residuos estandarizados significativos table(res.std.sig) ## res.std.sig ## FALSE TRUE ## 559 10 70 head(res.std[res.std.sig]) ## 32 41 74 129 136 172 ## 2.448381 6.124732 2.789193 -3.728307 4.144604 3.865401 res.ordeest <- sort(abs(res.std[res.std.sig]), decreasing =TRUE) # mostramos só os máis altos head(res.ordeest) ## 298 41 136 172 129 74 ## 7.479692 6.124732 4.144604 3.865401 3.728307 2.789193 #Gráficos para os residuos do modelo plot(res, cex =0.6) abline(h=c(-2,2), col ="red") 71 ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● 0 100 200 300 400 500 −4 −2 0 2 4 6 Index res #Residuos con valor absoluto elevado signif <- which(abs(res) >2) plot(res[signif], type ="n") text(1:length(signif), res[signif], label = signif, cex =0.4) #Gráficos de residuos empregando DHARMa #install.packages("DHARMa",repos = "http://cran.us.r-project.org") #library(DHARMa) #install.packages("glmmTMB",repos ="http://cran.us.r-project.org") #library(glmmTMB) 72 #resDhar <- simulateResiduals(modelo6) #plot(resDhar) #Observamos valores dos leverages hat.valores <- hatvalues(modelo6) head(hat.valores) ##1234 ## 2.745494e-04 1.365718e-03 6.419484e-06 1.171571e-02 ## 5 6 ## 2.180533e-04 6.057120e-02 #Conclusión predicciones <- ifelse(test = modelo6$fitted.values >0.5,yes =1,no =0) matriz_confusion <- table(modelo6$model$diagnosis, predicciones, dnn =c("observaciones","predicciones")) matriz_confusion ## predicciones ## observaciones 0 1 ## 0 343 14 ## 1 20 192 library(vcd) ## Loading required package: grid 73 2 4 6 8 10 −4 −2 0 2 4 6 Index res[signif] 32 41 74 129 136 172 206 264 298 561 mosaic(matriz_confusion, shade = T, colorize = T, gp =gpar(fill =matrix(c("green3","red2","red2","green3"), 2,2))) 74 predicciones observaciones 1 0 0 1 #CODIFICACIÓN 0 1 DA VARIABLE RESPOSTA datos$diagnosis <- as.character(datos$diagnosis) datos$diagnosis <- as.numeric(datos$diagnosis) ### plot(fitted.values(modelo6),datos$diagnosis) 75 ●●●●●● ●● ●●● ● ●● ● ●● ●● ●●● ● ●●●● ●●● ●● ●●●●● ● ●●● ● ●●● ● ● ● ● ●●●● ●● ● ●● ●●●● ● ● ● ● ●●●● ● ● ●● ● ● ● ●● ●● ● ●● ● ●● ● ● ●● ● ● ● ● ● ●●● ●● ●●●● ● ●● ● ●● ● ●●●●● ● ●● ● ●● ●●● ● ● ● ● ● ●● ● ●● ●● ● ●● ● ●●●● ● ● ●● ●● ●●●● ● ●●● ● ●● ● ● ●● ● ● ●● ● ● ●●●● ● ●● ●●● ● ● ● ● ●●● ● ●● ● ● ● ● ● ●● ● ●●● ● ● ● ● ●● ● ● ●●●● ●● ●● ● ●● ● ● ●● ●● ● ● ● ● ● ● ● ●● ● ● ●●● ● ● ●●●●● ● ● ●●●● ●●●●●● ●● ●● ● ●●●●● ● ● ● ●● ● ● ● ● ● ●● ●● ●● ●● ●●●●●●● ● ●● ● ● ● ●●●●●●●●●●●●●● ● ●●● ● ● ● ●●●● ●●● ●●●● ● ● ● ● ● ●●● ● ●●● ●●●● ●●● ● ● ●●●●● ●● ●● ●● ● ●●● ● ●● ● ●● ●● ● ●●●●● ● ●●● ● ●● ●● ●● ●●●● ● ●●●●● ●● ● ●●●● ● ● ●● ● ●●● ●● ●●●●●●● ● ● ●● ● ● ●●●●● ● ●● ● ● ● ● ● ● ● ● ● ●● ●●●●● ●● ●●●● ●● ● ●● ●●●●● ●●● ● ●● ●● ●●● ● ● ● ● ● ● ●● ●●● ●● ● ● ● ● ●●●● ● ● ●● ● ● ● ● ●● ●●● ● ● ●●● ●● ●● ●●● ● ● ●● ●●●● ●●●● ●●●●●●● ●● ●●●● ●● ●● ●●●●●● ● 0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0 fitted.values(modelo6) datos$diagnosis #Aplicación do modelo de regresión de Poisson #install.packages("AER", repos = "https://cran.r-project.org/package=AER") library("MASS") ## ## Attaching package: ’MASS’ ## The following object is masked from ’package:dplyr’: ## ## select 82 #library(DHARMa) #install.packages("glmmTMB") #library(glmmTMB) #resDhar <- simulateResiduals(modelo) #plot(resDhar) #Sobredispersión modelo library(AER) ## Loading required package: lmtest ## Loading required package: zoo ## ## Attaching package: ’zoo’ ## The following objects are masked from ’package:base’: ## ## as.Date, as.Date.numeric ## Loading required package: sandwich ## Loading required package: survival dispersiontest(modelo) ## ## Overdispersion test ## ## data: modelo ## z = 7.9502, p-value = 9.312e-16 ## alternative hypothesis: true dispersion is greater than 1 ## sample estimates: ## dispersion ## 2.027443 #Modelo de regresión binomial negativa modelo_bn<-glm.nb(visits~age+income,data = DoctorVisits) summary(modelo_bn) ## 83 ## Call: ## glm.nb(formula = visits ~ age + income, data = DoctorVisits, ## init.theta = 0.4239459284, link = log) ## ## Deviance Residuals: ## Min 1Q Median 3Q Max ## -0.8230 -0.6932 -0.5995 -0.5617 3.6708 ## ## Coefficients: ## Estimate Std. Error z value Pr(>|z|) ## (Intercept) -1.62405 0.10929 -14.860 < 2e-16 *** ## age 1.37043 0.16723 8.195 2.5e-16 *** ## income -0.32376 0.09872 -3.280 0.00104 ** ## --- ## Signif. codes: ## 0 '***'0.001 '**'0.01 '*'0.05 '.'0.1 ' ' 1 ## ## (Dispersion parameter for Negative Binomial(0.4239) family taken to be 1) ## ## Null deviance: 3099.2 on 5189 degrees of freedom ## Residual deviance: 2992.9 on 5187 degrees of freedom ## AIC: 7076.4 ## ## Number of Fisher Scoring iterations: 1 ## ## ## Theta: 0.4239 ## Std. Err.: 0.0309 ## ## 2 x log-likelihood: -7068.3970