Modelización logística II
Alex De la Puente
1/12/2020
REGRESIÓN ORDINAL
Los datos ordinales son datos que toman valores que tienen un rden natural, para los cuales los intervalos entre valores no tiene por que tener un significado. Un ejemplo de datos ordinales: estado de salud: Excelente, muy bueno, bueno, regular, lao; encuestas de satisfacción: Mucho Poco, Nada,etc.
En este caso tendríamos una variable respuesta con una distribucón ordinal, y los modelos que tendremos que ajustar van a ser los modelos de regresión ordinal.
Cuando disponemos de datos ordinales podemostomar varias decisiones: + recodificar los datos, pasar a una variable dicotómica y utilizar la regresión logística + Analizar cada una de las categorías
En cualquiera de los dos casos tendríamos modelos faciles de interpretar pero inficientes estadísticamente hablando, por que estamos perdiendo potencia, estámos perdiendo tamaño muestral dado que no utilizan toda la información disponible en los datos.
Los modelos de regresión ordinal lo combinan todo en un unico modelo.
Modelos de Odds proporcionales
Estos modelos asumen que el riesgo (Odds) es el mismo cuando cambiamos de categoría. con lo cual a nivel de estimación de parametros es mucho ams sencillo que el modelo multinomial, porque el modelo mutinomial estimabamos diferentes parametros para cada una de nuestras funciones logit. En este caso no. Por que este modelo asume que tenemos esa proporcinalidad con respecto a los odds.
Lo que el modelo hace es transformar la escala ordinal a varios puntos de corte binarios. el numero de puntos de corte siempre va a ser el numero de categorias que tenemos en nuestra variable respuesta -1.
Podemos pensar que estos puntos de corte son los umbrales que necesitamos cruzar para pasar de una categoría a la siguiente (más alta).
En este modelo no estimamos la probabilidad de algo, sino la probabilidad de observar un resultado o algo menor.
Sea Y una variable respuesta que toma valores \(k = 1, . . . ,K\) Definimos \(p˘k = P(Y ≤ k|X)\), para \(k = 1, . . . ,K − 1\). Las \(p˘k\) son las probabilidades acumuladas y \(k\) los puntos de corte.
Los modelos que estimamos son los siguientes Para cada uno de los \(p_K\) obtenemos una función logit.
\[logit(\breve{p}_1)=log(\frac{p_1}{1-p_1}) =\delta_1 - \beta_1X1 + ... + \beta_qX_q )\] \[logit(\breve{p}_2)=log(\frac{p_1+p_2}{1-p_1-p_2}) =\delta_2 - \beta_1X1 + ... + \beta_qX_q )\] \[logit(\breve{p}_1)=log(\frac{p_1+p_2+...+p_k}{1-p_1-p_2-...-p_k}) =\delta_k - \beta_1X1 + ... + \beta_qX_q )\] Y se cumple que \(1 = p_1 + p_2 + . . . + p_K + p_K+1\). El modelo asume que los coeficientes que acompañan a las covariables son iguales, y la única diferencia se da en la ordenada en el origen, δk .
Los betas aqui son todos iguales a excepcion del intercepto Y eso no nos pasaba en la regresion multinomial. Con lo cual el numero de parametros que vamos a estimar en este caso en mucho menor.
Las \(\breve{p}\) me dicen la probabilidad acumulada de estar en esas categorías.
El modelo se llama de odds-proporcionales ya que el odds para cualquier k es proporcional a \(e^{−β1X1+...−βqXq}\)
Si una cierta combinación de variables explicativas duplica el odds de estar en la categoría 1, también dobla el odds de estar en la categoría 2, etc. Los odds solo se diferencian en el punto de partida (cuando todas las covariables son cero), es decir, en \(e^{δ_k}\).
EJEMPLO
En este ejemplo tenemos datos de un estudio cuyo objetivo es buscar los factores que influyen en un estudiante a la hora de decidir si va a realizar estudios de postgrado o no. Se dipone de una muestra de 400 estudiantes de grado a los que se les preguntó si era: nada probable, algo problable, o muy probable que realizaran estudios de postgrado. Se recogieron datos sobre el nivel de educación de sus padres (pared 0=bajo, 1=alto), si el centro en el que estaban estudiando era público o privado (public), y su nota media del expediente académico (gpa).
## Cargamos librerías
library(MASS)## Warning: package 'MASS' was built under R version 4.0.3library(foreign)
library(VGAM)## Warning: package 'VGAM' was built under R version 4.0.3## Loading required package: stats4## Loading required package: splineslibrary(effects)## Warning: package 'effects' was built under R version 4.0.3## Loading required package: carData## Warning: package 'carData' was built under R version 4.0.3## lattice theme set by effectsTheme()
## See ?effectsTheme for details.## Lectura de los datos
setwd("C:/Users/Alex/Desktop/MÁSTER MATEMÁTICAS/MODELIZACIÓN ESTADÍSTICA/IRANZU BARRIO")
postgrado=read.dta(file="postgrado.dta")
postgrado$pared=factor(postgrado$pared)# le digo a R que son variables categóricas
postgrado$public=factor(postgrado$public)# le digo a R que son variables categóricas
summary(postgrado)## apply pared public gpa
## unlikely :220 0:337 0:343 Min. :1.900
## somewhat likely:140 1: 63 1: 57 1st Qu.:2.720
## very likely : 40 Median :2.990
## Mean :2.999
## 3rd Qu.:3.270
## Max. :4.000attach(postgrado)La funcion que vamos a utilizar para ajustar los modelos de regresión ordinal es la función polr.
polr es una función de la librería MASS (Documentación) que ajusta modelos de regresión para variables respuesta con distribución ordinal, con funciones enlace logit o probit.
Es importante tener en cuenta que esta función usa una parametrización que implica que los odds ratio son de estar en una categoría superior.
Los pasos para ajustar el modelo son los mismos que en la regresión multinomial.
ord1=polr(apply~pared,postgrado)
summary(ord1)##
## Re-fitting to get Hessian## Call:
## polr(formula = apply ~ pared, data = postgrado)
##
## Coefficients:
## Value Std. Error t value
## pared1 1.127 0.2634 4.28
##
## Intercepts:
## Value Std. Error t value
## unlikely|somewhat likely 0.3768 0.1103 3.4152
## somewhat likely|very likely 2.4519 0.1826 13.4302
##
## Residual Deviance: 722.7903
## AIC: 728.7903Por un lado tengo un coeficiente que me dice el beta asociado a la variable pared (1.127). La primera cosa que nos tenemos que fijar es cual es la categoría asociada para la cual estamos estimando este estimador.En este caso como me dice que es pared1, mi categoría de referencia es la 0, “tienen estudios bajos”. Y el beta que estamos estimando aqui es para “tienen estudios altos”. Este es el beta asociado a esa variable.
¿Por qué me sale solo un beta? por que es el mismo en las dos funciones logit. Pero tengo dos interceptos diferentes \(\delta_1=0.3768\) y \(\delta_2=2.4519\)
No obtenemos p-valores.En este caso me daría un poco igual saber si ese beta es significativo o si mi variable es significativa, por que solo tengo un beta, por lo que la pregunta es la misma.
Hacemos lo mismo para la variable public:
ord2=polr(apply~public,postgrado)
summary(ord2)##
## Re-fitting to get Hessian## Call:
## polr(formula = apply ~ public, data = postgrado)
##
## Coefficients:
## Value Std. Error t value
## public1 0.1833 0.2838 0.6459
##
## Intercepts:
## Value Std. Error t value
## unlikely|somewhat likely 0.2252 0.1076 2.0924
## somewhat likely|very likely 2.2231 0.1718 12.9429
##
## Residual Deviance: 740.7912
## AIC: 746.7912Y gpa:
ord3=polr(apply~gpa,postgrado)
summary(ord3)##
## Re-fitting to get Hessian## Call:
## polr(formula = apply ~ gpa, data = postgrado)
##
## Coefficients:
## Value Std. Error t value
## gpa 0.7249 0.2493 2.908
##
## Intercepts:
## Value Std. Error t value
## unlikely|somewhat likely 2.3748 0.7570 3.1372
## somewhat likely|very likely 4.3998 0.7815 5.6301
##
## Residual Deviance: 732.60
## AIC: 738.60podemos usar el test de wall para ver si los modelos univariantes son significativos o no. ¿Que hacemos ahora?
¿Que pasaría si tuviera una variable categórica con 4 categorías?, tendría 3 betas. Estudio si las betas son significativas con el test de wall, imaginaros que me da que una es significativa y dos no lo son.–> necesito un test de la razon de verosimilitud.
tengo que comparar la verosimilitud de ese modelo univariante que he ajustado, con la verosimilitud del modelo nuevo. Si soy capaz de reducir mucho esa verosimilitud querra decir que me merece la pena introducir esa variable en el modelo.
# Significación estadística . Test razòn de verosimilitud
ord0=polr(apply~1,postgrado)# ajusto el modelo nulo, le pongo un 1 después del gusanito por que no tengo ninguna covariable.
anova(ord0,ord1)# lo comparo con los demás ## Likelihood ratio tests of ordinal regression models
##
## Response: apply
## Model Resid. df Resid. Dev Test Df LR stat. Pr(Chi)
## 1 1 398 741.2053
## 2 pared 397 722.7903 1 vs 2 1 18.41498 1.776559e-05Y obtengo un p-valor de 1.776559e-05<0.05–> Rechazo la hipótesis nula.(betas=0) Me quedo con el modelo grande.
Quiere decir que la variable es significativa, me quedo con ese modelo.
anova(ord0,ord2)## Likelihood ratio tests of ordinal regression models
##
## Response: apply
## Model Resid. df Resid. Dev Test Df LR stat. Pr(Chi)
## 1 1 398 741.2053
## 2 public 397 740.7912 1 vs 2 1 0.4140527 0.5199196p-valor de 0.5199196, variable no significativa.
anova(ord0,ord3)## Likelihood ratio tests of ordinal regression models
##
## Response: apply
## Model Resid. df Resid. Dev Test Df LR stat. Pr(Chi)
## 1 1 398 741.2053
## 2 gpa 397 732.6000 1 vs 2 1 8.605321 0.0033518220.003351822<0.05, variable significativa.
Tenemos que influyen a la hora de hacer un potgrado o no, la nota media del expediente y la nota media de los padres. En cambio el tipo de escuela no influye.
Ahora contruímos un modelo multivariante con las dos variables significativas.
####################################
#AJUSTE DEL MODELO MULTIPLE
###################################
ord4=polr(apply~pared+gpa,postgrado)
summary(ord4)##
## Re-fitting to get Hessian## Call:
## polr(formula = apply ~ pared + gpa, data = postgrado)
##
## Coefficients:
## Value Std. Error t value
## pared1 1.0457 0.2656 3.937
## gpa 0.6042 0.2539 2.379
##
## Intercepts:
## Value Std. Error t value
## unlikely|somewhat likely 2.1763 0.7671 2.8370
## somewhat likely|very likely 4.2716 0.7922 5.3924
##
## Residual Deviance: 717.0638
## AIC: 725.0638Ahora tengo un beta más. \[logit(\breve{p}_1) =\delta_1 - \beta_1X1=2.1763+1.0457X_1 + 0.6042X_2 )\] \[logit(\breve{p}_2) =\delta_2 - \beta_1X1=4.2716 +1.0457X_1+ 0.6042X_2)\] Ahora otra vez yo voy a querer saber si esas dos variables son significativas en ese modelo.
Ahora puedo comparar el modelo multivariante con el modelo que tiene solo la variable pared. Con esto veo si la variable gpa es significativa dentro del modelo multivariante.
anova(ord1,ord4)## Likelihood ratio tests of ordinal regression models
##
## Response: apply
## Model Resid. df Resid. Dev Test Df LR stat. Pr(Chi)
## 1 pared 397 722.7903
## 2 pared + gpa 396 717.0638 1 vs 2 1 5.726514 0.016710620.01671062<0.05–>TRUE: Rechazo H0, la variable es significativa y me quedo con el modelo multivariante
Por otro lado, como tengo solo un beta asociido a cada una de las variables, puedo mirar su significación con un test de wall o equivalentemente con el intervalo de confianza (mirando si el cero esta dentro del intervalo del beta, o si el 1 no está dentro del intervalo de confianza del odds ratio)
# Intervalos de confianza
ci=confint(ord4)## Waiting for profiling to be done...##
## Re-fitting to get Hessianexp(cbind(OR=coef(ord4),ci))## OR 2.5 % 97.5 %
## pared1 2.845412 1.693056 4.806454
## gpa 1.829873 1.115180 3.022320En este caso, como el 1 no esta dentro de ninguno de los intervalos de confianza para los dos modelos.
La posibilidad de pasar de nada probable a algo o muy probable, es 2.8 veces mayor para los alumnos cuyos padres tienen estudios que para los que no los tienen (también sería 2.8 veces más probable para de nada o algo probable a muy probable).
La posibilidad de pasar de nada probable a algo o muy probable, es 1.8 veces mayor por cada unidad de incremento en la nota del expediente (también sería 1.8 veces más probable para de nada o algo probable a muy probable).
Hemos asumido que el cambio en odds se mantiene constante. Es decir, que yo cuando paso de nada a algo, probable o muy probable, tengo un odds que se me mantiene cuando paso de nada o algo, a muy probable.
Tenemos que comprobar que eso que hemos asumido, realmente es cierto.
Debemos comprobar que se satisfacen las condiciones del modelo de odds proporcionales. Es decir, si creamos dos variable dicotómicas correspodientes a las dos transiciones (dos variables respuesta) y ajustáramos dos modelos de regresión logística, los coeficientes en ambos casos serían muy similares.
Para comprobarlo podemos utilizar la librería VGAM:
Para ello existe la funcion vglm de la libreria vgam que nos permite ajustar el modelo asumiendo que tengo riesgos proporcionales o que no los tengo.
# Comprobación Odds Proporcionales
apply2 = ordered(postgrado$apply)
m0 <- vglm(apply2~pared+gpa,family=cumulative(parallel=T),
data=postgrado)#si quiero que los odds ratio sean proporcionales--> parallel=TRUE
m1 <- vglm(apply2~pared+gpa,family=cumulative(parallel=F),
data=postgrado)ahora tenemos que calcular el valor del estadístico
# logLik(m1) devuelve el valor de la verosimilitud para el perimer modelo
# logLik(m2) devuelve el valor de la verosimilitud para el perimer modelo
test.po <- 2*logLik(m1)-2*logLik(m0)
df.po <- length(coef(m1))-length(coef(m0))# necesito conocer los grados de libertad estoy teniendo en cuenta. Si yo miro la longitud de mis coeficientes en un modelo y otro y los resto, lo que obtengo esla diferencia de parametros estimados que tengo entre un modelo y otro modelo.
test.po # valor del estadístico## [1] 0.8357332df.po # grados de libertad## [1] 21-pchisq(test.po,df=df.po) #sigue una distribución chi cuadrado.## [1] 0.6584501¿Por qué hago 1-pchisq? Porque esta probabilidad asume lowertail=T. Es decir, me calcula que probabiidad hay de una chi cuadrado tome un valor menor o igual que G.
pchisq(test.po,df=df.po,lower.tail = F)## [1] 0.6584501Da lo mismo.
REGRESIÓN DE POISSON
Vamos a tener modelos de Poisson cuando tengamos una variable respuesta que tenga esa distribución. ¿Y eso cuando va a pasar?
Gran cantidad de datos son recogidos por cientícos, médicos, o empresas como datos de conteo. En general, este tipo de datos aparecen en cuatro formatos distintos:
Frecuencias, donde contamos cuantas veces ocurre algo, pero no sabemos cuantas veces no ocurre.
Proporciones, donde ambos, el número de veces que ocurre algo, y el tamaño total del grupo es conocido.
Datos por categorías, donde la variable cuenta cuántos individuos hay en cada nivel de una variable categórica.
Datos binarios, donde recogemos la presencia o ausencia de una característica.
Hay cuatro razones por las que es erroneo utlizar un modelo de regresión lineal (asumimos normalidad de la variable respuesta) para datos de conteo:
Puede dar lugar a predicciones negativas.
La varianza de la variable respuesta no es independiente de la media.
Los errores no siguen una distribución Normal.
Los ceros que aparecen en la variable respuesta dan problemas a la hora de transformar la variable.
Por todas esas razones, cuando tenemos frecuencias y cuando tenemos conteos, el modelo de interés es el modelo de la regresión de Poisson.
La principal diferencia entre la distribución de Poisson y la Binomial, es que, aunque ambas cuentan el número de veces que ocurre algo, en la distribución de Poisson no sabemos cuantas veces no ocurrió, y en la Binomial sí lo sabemos. Sea Y una variable con distribución de Poisson de parámetro λ (Y : Pois(λ)) . Entonces su función de probabilidad viene definida por: $$f (y) = P(Y = y){ \[\begin{array}{lcc} \frac{λ^ye^{−λ}}{y!} &\space y ∈ {0, 1, 2, . . .} \\ \\ 0! & en\space otro\space caso\\ \end{array}\].$$ Si \(Y : Pois(λ)\), se cumple que \(µ = E(Y) = λ\) y \(σ^2 = Var(Y) = λ\), es decir, \(E(Y ) = Var(Y )\).
Relación entre la distribución de Poisson y la Binomial
La variable aleatoria de Poisson tiene una gran variedad de aplicaciones en diversas áreas ya que puede utilizarse como una aproximación de una variable aleatoria binomial con los parámetros (n, p) cuando n es grande y p es lo suficientemente pequeño para que λ = np sea moderada.
Por lo tanto, lo que estamos haciendo es una aproximación: λ = np para n grande y p pequeño. Nota: si utilizamos una regresión logística con datos agrupados obtenemos los mismos coeficientes que si utilizamos una regresión de Poisson, pero la interpretación de ambas sólo es similar cuando la prevalencia (p = P(Y = 1)) es baja.
Para más detalle sobre variables aleatorias y la relación entre ellas ver Ross (2014).
Planteamiento del modelo
Sea %Y_i%la variable respuesta que mide el número de veces que le ha ocurrido un determinado evento al individuo \(i\), \(i = 1, . . . , n\), y que deseamos explicarlo en función de un vector de variables predictoras \(X_i\).Si \(Yi: Pois(µi)\), debemos relacionar \(µi\) con el vector de covariables \(Xi\) . Dado que se debe cumplir que \(µ ≥ 0\), entonces la función link que se utiliza habitualmente es la función log.
\[log(E(Y )) = log(µ) = log(λ) = β0 + β1X1 + . . . + βqXq\] De modo que es una transformación de la media la que tiene una relación lineal con las covariables, y exp(βj) representa un efecto multiplicador del j-ésimo predictor sobre la media, ya que al incrementar Xj en una unidad, la media queda multiplicada por un factor exp(βj). Una de las ventajas de utilizar este tipo de modelo es que, en general, en el caso de datos de conteo, los predictores son multiplicativos y no aditivos, es decir, observamos conteos pequeños para efectos pequeños y conteos grandes para efectos grandes (o viceversa).
Estimación de parámetros
Dado que la variable respuesta Y sigue una distribución de Poisson, la función log-verosimilitud es la siguiente: \[l(\beta) = \sum_{i=1}^{n}(y_{ik}x_i^T\beta- exp(x_i^T\beta)-log(y_i!))\] El estimador \(\hat{β}\) del vector de parámetros \(β\), será por tanto aquel que maximice la log-verosimilitud. Si tomamos derivadas con respeto a los \(βj\), \(j = 1, . . . , q\) e igualamos a cero obtenemos: \(x^ty = x^t exp(x^tβ)\), donde \(x\) es la matriz del modelo.
ES EXACTAMENTE LO MISMO QUE LA REGRESIÓN LOGISTICA Y LINEAL ¿Qué es lo que cambia? Que mi definición de la función de verosimilitud cambia y por tanto de la log verosimilitud cambia. Porque mi función de probabilidad o mi función de densidad cambia en función que tiene la distribución que tiene mi variable respuesta.Pero todo lo demás es exactamente igual.
Al igual que ocurre con la estimación de los parámetros del modelo de regresión logística, no existe una fórmula explícita para el estimador \(\hat{β}\), por lo que es necesarioutilizar métodos numéricos para su resolución (iterative reweighted least squales).
Deviance
Para un modelo de regresión de Poisson la Deviance viene dada por: \[D = 2\sum_{i=1}^{n}(y_ilog(y_i/\hat{µ}_i)) − (y_i − \hat{µ}_i),\] Donde, \(\hat{µ}_i = exp(\hatβ_0 + \hatβ_1x_{1i} + . . . + \hatβ_qx_{qi})\).
La deviance es una medida de lo bien que el modelo se ajusta a los datos - si el modelo se ajusta bien, los valores observados se acercarán a su media estimada - causando que ambos términos en la expresión de la Deviance sean pequeños, y por lo tanto que la deviance misma sea pequeña.
La fórmula para la deviance anterior se deriva como el cociente de verosimilitud comparando el modelo estimado con el llamado modelo saturado. El modelo saturado es el modelo para el cual los valores estimados del modelo coinciden exactamente con los valores observados (de la misma forma que en la regresión logística).
Comparación de modelos anidados
En el modelo de regresión de Poisson se puede utilizar la misma inferencia asintótica que para el modelo logístico. Podemos evaluar la bondad de ajuste del modelo ajustado comprobando la deviance del modelo frente a una distribución \(χ^2\) con grados de libertad iguales a los del modelo.
Por otro lado, podemos comparar modelos anidados tomando las diferencias de las deviance de ambos y comparándolo con una distribución \(χ^2\) con grados de libertad iguales a la diferencia en el número de parámetros de los dos modelos. De esta forma, podremos estudiar la significación estadística de cada una de las variables del modelo.
Notad que de la teoría de la probabilidad obtenemos que a medida que el tamaño de la muestra se hace grande \((n → ∞)\), la diferencia de las deviance sigue una distribución \(χ^2\) bajo la hipótesis nula de que el modelo más simple está correctamente especificado.
###Sobre el test de la razón de verosimilitud
Pawitan afirma en su libro In All Likelihood que el test de razón de verosimilitud (o diferencia en deviance) para medir la bondad de ajuste es aceptable para datos de Poisson siempre que las medias no sean demasiado pequeñas.
Cuando la media es grande, una distribución de Poisson se aproxima a la distribución normal, y por lo tanto la función logaritmo como función link es aproximadamente lineal.
En la práctica, consideraremos que este test es válido siempre que el número de casos en la mayoría de estratos sea > 5.
Si no, puedo tener problemas al asumir que mi distribución es Gi cuadrado
Bondad de Ajuste y residuos del modelo
Para ver si el modelo se ha ajustado correctamente, utilizaremos los residuos estandarizados: \[ri =\frac{y_i − \hat{µ}_i}{\sqrt{\hat{µ}_i}}\] Si el modelo está ajustado correctamente, estos residuos deben seguir una distribución N(0, 1) por lo tanto la mayoría de ellos deben estar en el intervalo (−2, 5; 2, 5).
EJEMPLO
Un experimento agrícola consta de 90 parcelas de pradera, de 25m2 cada una, que se diferencian entre sí, en biomasa, pH del suelo, y la riqueza en especies.
Es sabido que la riqueza de especies decrece cuando aumenta la biomasa, pero la cuestión de interés aquí es si esa disminución cambia con el pH del suelo.
Las parcelas se clasificaron de acuerdo a tres niveles de pH: alto (= 2), medio (= 1) y bajo (= 0). La variable respuesta es el número de especies por parcela, y las variables predictoras son la biomasa media medida en el mes de Junio, y la variable categórica correspondiente al pH del suelo.
# Lectura de los datos
setwd("C:/Users/Alex/Desktop/MÁSTER MATEMÁTICAS/MODELIZACIÓN ESTADÍSTICA/IRANZU BARRIO")
especies <- read.table("especies.txt",header=TRUE)
especies$pH<-as.factor(especies$pH)
names(especies)## [1] "pH" "Biomasa" "Especies"summary(especies)## pH Biomasa Especies
## 0:30 Min. :0.05017 Min. : 2.00
## 1:30 1st Qu.:1.44132 1st Qu.:12.25
## 2:30 Median :3.10836 Median :18.50
## Mean :3.55701 Mean :19.46
## 3rd Qu.:5.08570 3rd Qu.:25.75
## Max. :9.98177 Max. :44.00attach(especies)# Análisis exploratorio
plot(Biomasa,Especies,type="n")
points(Biomasa[pH==0],Especies[pH==0])
points(Biomasa[pH==1],Especies[pH==1],pch=2,col=2)
points(Biomasa[pH==2],Especies[pH==2],pch=3,col=3)
legend(4,45,legend=c("Bajo","Medio","Alto"),pch=c(1,2,3),col=c(1,2,3))
abline(lm(Especies[pH==0]~Biomasa[pH==0]), col=1)
abline(lm(Especies[pH==1]~Biomasa[pH==1]), col=2)
abline(lm(Especies[pH==2]~Biomasa[pH==2]), col=3)
Ahi esta dibujado cual sería si ajustase un modelo de regresión lineal en el que yo obtendría las especies por la biomasa en función de los phs.
No es correcto hacer un modelo de regresión lineal puesto que si yo extiendo la linea, obtendría un valor negativo de especies a medida que aumente el valor de la biomasa, y eso no tiene sentido.
# Ajustamos los Modelos univariantes
m0=glm(Especies~Biomasa,family=poisson)# le tenemos que decir a que familia pertenece nuestra variable respuesta, por que sino va a ajustar un modelo lieneal
summary(m0)$coefficients## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 3.18409377 0.039159103 81.311714 0.000000e+00
## Biomasa -0.06444054 0.009837528 -6.550481 5.735197e-11La función glm nos vale para la regresión logística, nos vale para la regresión lineal.
m1=glm(Especies~pH,family=poisson)
summary(m1)$coefficients## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 2.4481274 0.05368221 45.604076 0.000000e+00
## pH1 0.5476049 0.06744215 8.119624 4.676311e-16
## pH2 0.8402745 0.06423050 13.082173 4.163671e-39Si quiero hacer el test de la razón de verosimilitud, de nuevo uso anova
anova(m1,test='Chi')## Analysis of Deviance Table
##
## Model: poisson, link: log
##
## Response: Especies
##
## Terms added sequentially (first to last)
##
##
## Df Deviance Resid. Df Resid. Dev Pr(>Chi)
## NULL 89 452.35
## pH 2 187.23 87 265.12 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Lo siguiente es ajustar un modelo con las dos variables, donde le introducimos Biomasa y ph. Y ajustamos un modelo adicional donde le introduzco la interacción entre la biomasa y el ph.¿Por qué? Porque hemos visto que sabemos que la biomasa influye en el numero de especies. Lo que yo quiero saber es si esa inluencia de la biomasa en el número de especies es distinta en función del pH en el que estoy. Es decir, si esa influencia de la masa cambia por 1 cuando el ph en neutro, cambia por 2 cuando el ph es alto o cambia por 3 cuando el p es bajo. O si por el contrario me da igual el ph que tiene mi suelo, que la influencia de la biomasa en igual
# Ajuste del modelo multiple
especies1=glm(Especies~Biomasa+pH,family=poisson)
especies2=update(especies1,~.+Biomasa:pH)
summary(especies2)$coefficients## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 2.9525524 0.08239758 35.833000 3.383847e-281
## Biomasa -0.2621580 0.03803334 -6.892846 5.468701e-12
## pH1 0.4841087 0.10723006 4.514673 6.341455e-06
## pH2 0.8155712 0.10283673 7.930739 2.178459e-15
## Biomasa:pH1 0.1231362 0.04269689 2.883961 3.927073e-03
## Biomasa:pH2 0.1550282 0.04003169 3.872636 1.076643e-04anova(especies2,especies1,test='Chi')## Analysis of Deviance Table
##
## Model 1: Especies ~ Biomasa + pH + Biomasa:pH
## Model 2: Especies ~ Biomasa + pH
## Resid. Df Resid. Dev Df Deviance Pr(>Chi)
## 1 84 83.201
## 2 86 99.242 -2 -16.04 0.0003288 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Trás ajustar los modelos univariantes tanto para la biomasa como para el pH, vemos que en ambos casos tenemos modelos significativos. En un modelo univariante la biomasa es significativa, y en otro modelo univariante, el pH es significativo.
Tenemos 3 categorías para el pH, por tanto debemos obtener 2 betas para esa variable. Una vez que vemos que las dos variables son estadísticamente significativas, nos plantemos la necesidad de incluirlas en un modelo multivariante, donde ajustamos un modelo aditivo con la biomasa y el pH (especies1). Y una vez vemos que ese modelo es significativo nos planteamos la necesidad de introducir una interacción.
¿Qué implica en un modelo de este tipo que una interacción sea significativa?
En un modelo aditivo (especies1) yo lo que asumo es que el cambio que yo tengo al pasar de un pH bajo a un pH medio es el mismo independiente del nivel de Biomasa que yo tenga. O equivalentemente, el efecto que yo tengo al añadir en una parcela una unidad más de biomasa es el mismo en parcelas con pH bajo, con pH medio y en parcelas con pH alto.
Si yo tengo un modelo aditivo en el que no tenemos interacción donde no hay interacción. Yo tengo la variable Biomasa, y mi variable respuesta es el logaritmo del número de especies, yo lo que estoy buscando es el beta que me marca la pendiente de la biomasa.
si tengo un modelo aditivo, ese beta1 que me indica la pendiente de ese modelo lineal, es la misma, cuando estoy en el pH bajo, cuando estoy en el pH neutro, y cuando estoy en el pH alto. ¿Que cambia? La ordenada en el origen.
Sin embargo, cuando meto una interacción en el modelo, esa pendiente me va a cambiar en cada uno de los pHs. Es decir, el comportamiento de mi biomasa con el numero de especies es diferente con cada uno de los pHs.
Eso quiere decir que en un modelo con interacción voy a tener que estimar muchos más parámetros.
En un modelo aditivo yo hubiese tenido log(especies) = b0 + b1(biomasa) + b2(pHmedio) + b3(phalto). 4 betas. Sin emabrgo, en el modelocon interacción tengo 2 betas mas. Esos 2 betas más corresponden a la interacción de la biomasa con el phmedio, y ala interacción de la biomasa con el pHalto.
¿como tengo que reescribir este modelo en el caso que tenga la interacción? log(especies) = b0 + b1(biomasa) + b2(phmedio) + b3(biomasa:pHalto) +b4(biomasa:phmedio) + b5(biomasa:phalto)
¿Cómo interpretámos eso?
Si yo estoy en el ph bajo, por una unidad que aumen
El efecto de la biomasa cuando estoy en el pH medio ya no es solo b1, b1 +b4, por que un incremento de una unidad en la biomasa me aumenta (b1+b4) unidades el logaritmo de la exponencial. ¿Qué pasa si voy al pH alto? Sigo teniendo que el efecto de la biomasa (b1), pero cuando voy a ph alto, asociado a la biomasa tengo un b5, y si yo quiero analizar que efecto tiene la biomasa en el logaritmo de las especies, ese efecto es b1+b5.
Ahora, a mi no me interesa interpretar en terminos de el logaritmo del numero de especies, que es muy poco intuitivo. Me interesa hablar sobre el número de especies. Qué pasa con el número de especies.
Me voy al caso del modelo más sencillo: el modelo sin interacción. ¿Cómo vamos a interpretar el efecto de la biomasa?
Yo voy a considerar este modelo, para un valor de x: \(Biomasa=x\) y añado una unidad : \(Biomasa=x+1\)
¿cuál sería el logaritmo de las especies, si mi Biomasa toma el valor x? Voy a tener b0 + b1x unidades de biomasa + b2( si estoy en el ph medio) + b3( si estoy en el ph alto).
¿cual es la ecuación de mi modelo si quiero analizar la biomasa para el valor x+1?: Voy a tener b0 + b1(x+1) unidades de biomasa + b2(si estoy en el ph medio) + b3(si estoy en el ph alto).
Y hago la resta de estas dos ecuaciones. Es decir, yo lo que quiero es ver cómo cambia el logaritmo del número medio de especies cuando tengo una biomasa de x+1 y cuando tengo una biomasa de x.
Entonces si hago log(numero de especies medio cuando la biomasa = x+1) - log(numero de especies medio cuando la biomasa = x), si os fijais, tengo todo exactamente igual a excepcion de este termino de aqui, y me va a quedar b1(x+1)-b1x, con lo cual, me va a quedar b1.
Esto lo puedo reescribir cómo el (logaritmo del número de especies cuando la biomasa es x+1)/log(numero de especies medio cuando la biomasa = x)
o equivalente: num especies medio cuando biomasa x+1= e^b1+log(numero de especies cuando la biomasa es x)
eso quiere decir que cuando yo añado una unidad a la cantidad de biomasa, el numero de especies que espero se multiplica por e elevado a mi parámetro estimado en el modelo. es decir, el numero de especies que yo espero tener por una unidad mas de biomasa es e^b1 veces mas mayorque el numero de especies medio que tengo cuando mi biomasa es x.
Hemos dicho que tenemos una interacción significativa, lo cual quiere decir que el efecto de la biomasa va a ser diferente en cada uno de los niveles de ph. Es decir, cómo cambie el numero de especies en función del cambio en la biomasa, va a ser distinto si esty en un ph alto, si estoy en un ph bajo, o en un ph medio.
\[log(Especies) = 2.95 − 0.262Biomasa + 0.484pH_{medio} + 0.815pH_{alto} +0.123pH_{medio}Biomasa + 0.155pH_{alto}Biomasa\]
Si estoy en un caso de ph bajo (mi categoria de referencia): \[log(Especies) = 2.95 − 0.262Biomasa\] Lo que quiere decir que una unidad de incremento en la Biomasa me multiplica el numero de especies en \(e^{\beta_1}\).Es decir, una unidada adicional de Biomasa supondría un número de especies \(e^{-0.262} = 0.769\) unidades mayor en suelos con pH bajo. Equivalentemente, una unidad menor de Biomasa supondría un \(\frac{1}{0.769}=1.30\) menos de especies en suelos con pH bajo.
pH medio : \(log(Especies) = 2.95 − 0.262Biomasa + 0.484pH_{medio} +0.123pH_{medio}Biomasa\) Una unidad de Biomasa adicional supondría un número de especies \(e^{-0.262+0.123}=0.87\) unidades mayor en suelos con pH medio. Equivalentemente, una unidad menor de Biomasa supondría un \(\frac{1}{0.87}=1.15\) menos de especies en suelos con pH medio.
pH alto : \(log(Especies) = 2.95 − 0.262Biomasa + 0.815pH_{alto} + 0.155pH_{alto}Biomasa\) Una unidad de Biomasa adicional supondría un número de especies \(e^{-0.262+0.155}=0.898\) unidades mayor en suelos con pH alto. Equivalentemente, una unidad menor de Biomasa supondría un \(\frac{1}{0.87}=1.112\) menos de especies en suelos con pH alto.
Con esto no habriamos acabado de interpretar el modelo. Nos falta interpretar que pasa cuando pasamos de un ph bajo al medio, y de un ph bajo al ph alto.
–se queda pendiente.–ç
# Gráfico del modelo estimado
plot(Biomasa,Especies,type="n")
points(Biomasa[pH==0],Especies[pH==0])
points(Biomasa[pH==1],Especies[pH==1],pch=2,col=2)
points(Biomasa[pH==2],Especies[pH==2],pch=3,col=3)
legend(4,45,legend=c("Bajo","Medio","Alto"),pch=c(1,2,3),col=c(1,2,3))
legend(6,45,legend=c("Modelo con interacción","Modelo sin interacción"),lty=c(1,2),col=c(1,1))
ajustados1=fitted(especies1) #modelo sin interacción
ajustados2=fitted(especies2) #modelo con interacción
x=Biomasa[pH==0]
or=order(x)
y1=ajustados1[pH==0]
y2=ajustados2[pH==0]
lines(x[or],y1[or],lty=2)
lines(x[or],y2[or])
x=Biomasa[pH==1]
or=order(x)
y1=ajustados1[pH==1]
y2=ajustados2[pH==1]
lines(x[or],y1[or],col=2,lty=2)
lines(x[or],y2[or],col=2)
x=Biomasa[pH==2]
or=order(x)
y1=ajustados1[pH==2]
y2=ajustados2[pH==2]
lines(x[or],y1[or],col=3,lty=2)
lines(x[or],y2[or],col=3)
Las rayas discontinuas corresponden al modelo sin interacción, que correspondería a las rectas paralelas que tenemos si lo que ponemos en y, es el logaritmo de las especies.
Y las de linea continua son las del modelo con interacción. estas 3 curvas, si las dibujamos en la escala logaritmica, van a tender a cruzarse. No tienen por que cruzarse en el rango que observamos, pero si que vamos a ver que no son paralelas.
Observaciones
En el caso de datos de Poisson, la media y la varianza son iguales. Si eso no se cumple, no es cierto que la variable respuesta que estamos modelizando sigue una distribución de Poisson.Si los datos no cumplen esta hipótesis los errores estándar de los estimadores de losparámetros serán incorrectos.
Entonces lo que podemos hacer es analizar cual es esa varianza que estamos estimando con ese modelo que estamos ajustando.Para ello, podemos introducir un parámetro de dispersión, que haga que gracias a el la varianza y la media sean iguales. \(φ\), demodo que \(Var(Y ) = φµ\) (si \(φ = 1\) estamos en el caso Poisson). Dicho parametro se puede estimar como: \[\hat{φ} =\frac{Deviance}{n − q}\] ¿Cómo se estima ese parámetro? Se calcula la Deviance, y se divide por los grados de libertad de ese modelo que es el tamaño muestral que tenemos menos el numero de parametros que estamos estimando.
En el ejemplo anterior, Deaviance = 83.2 y n − q = 90 − 6 = 84, por lo que el parámetro es muy próximo a 1. Podemos utilizar en la función glm, la familia quasipoisson, que automáticamente estima el parámetro de dispersión. La estimación puntual de los parámetros será la misma pero cambiarán los errores estándar.
Podemos considerar que estamos ajustando bien un modelo de Poisson siempre y cuando ese parametro de dispersión estimado sea proximo a 1.
IMPORTANTE: La estimación puntual de los parametros si tenemos varianza igual a la media o si no la tenemos, va a ser muy parecidas, pero lo que nos va a cambiar es la variabilidad que estamos estimando para esos parámetros, es decir, los errores estándar asociados a esos parámetros.
Si no tenemos la misma varianza y la misma media, en lugar de decirle a R que tenemos una distribución de Poisson, le podemos decir que tenemos una distribución de quasi-poisson.
Podemos utilizar en la función glm, la familia quasipoisson, que automáticamente estima el parámetro de dispersión. La estimación puntual de los parámetros será la misma pero cambiarán los errores estándar.
Regersión de Poisson para tasas de incidencia
Hasta ahora hemos estado hablando del número de veces que ocurre algo, pero muchas veces no nos interesa el numero bruto como tal, sino la tasa. La tasa del numero de veces que ocurre algo pero en un intervalo de tiempo.
Podemos usar la regresión logística para estimar la prevalencia (proporción), pero no podemos usarla para estimar la incidencia (incidencia=probabilidad de que ocurra algo en un intervalo de timepo), ya que en este caso necesitamos utilizar el tiempo de exposición al riesgo, ya que en la tasas de incidencia medimos el número de sucesos en la unidad de tiempo.
Por que con la regresión logística no tenemos información del tiempo. Tenemos informacion de cuantas veces ocurre algo, y de cuantas veces no ocurre algo.
Si necesitámos analizar el tiempo de exposicion. (numero de accidentes en un año, numero de positivos en una semana…) entonces ya estoy hablando de tasas, porque estoy teniendo en cuenta esa exposición, el tiempo de exposición en el que yo estoy haciendo ese conteo.
Supongamos que la variable respuesta Y representa el número de eventos que aparecen con una tasa λ por unidad de tiempo de exposición y que la exposición es E (E se suele dar en términos de personas por año o 1000 personas por año) y va a determinar las unidades en las que se expresa la tasa.
(la tasa sería el número de eventos en un tiempo determinado)
Dado que nuestro interés es la tasa, λ y que esta puede depender de distintas covariables, podríamos ajustar el modelo: \(log(λ) = β0 + β1X1 + . . . + βqXq\)
PERO, CUIDADO! es Y quien sigue una distribución de Poisson y no λ. . . . Y nuestro objetivo siempre es estimar la experanza matemática (media) de Y , es decir E(Y). ¿Cómo podemos replantear el modelo cuando lo que a nosotros nos interesa es ajustar un modelo de Poisson para la tasa, pero sabiendo que esa tasa no sigue una distribución de poisson?
Dado que , Y : Pois(µ) y el número medio de eventos E(Y ) = µ = λ ∗ E (tasa x el tiempo de exposición), y de nuevo usando la función logaritmo como función link: \[log(µ) = log(λ ∗ E) = log(λ) + log(E) = β_0 + β_1X_1 + . . . + β_qX_q\]
Al log(E) se le llama offset y es una cantidad fija y no hay que estimar ningún parámetro para ella. Eso quiere decir que me puedo olvidar de esa exposición cuando ajusto el modelo. Me puedo olvidar relativamente: lo tengo que tener en cuenta, tengo que saber que tengo esa exposición ahi y que lo que yo estoy ajustando es una tasa, no es el número de eventos,pero no tengo que cambiar nada más en mi modelo.
Lo que le vamos a tener que decir al modelo es que nosotros estamos ajustando una tasa, y le tenemos que decir cual es ese tiempo de exposición, para el que nosotros estámos calculando esa tasa. Eso se lo metemos cómo un valor fijo, una constante fija y ajustamos el modelo teniendo en cuenta que estamos estimando para una tasa.
Es mas común ver modelos de Poisson para estimar tasas que para estimar numero de casos.
INTERPRETACIÓN
La interpretación de los parámetros va unida al concepto de riesgo relativo. Por ejemplo, supongamos que tenemos una sola covariable X, entonces: \[log(µ) = β_0 + β_1X\] o lo que es lo mismo, \[log(λ) = β_0 + β_1X − log(E) = β_0 + β_1X − offset\]
¿Cómo interpretamos el modelo cuando lo ajustamos para tasas y no para el número de casos? Sean λ0 la tasa de incidencia cuando X toma el valor x0, y λ1 la tasa de incidencia cuando X toma el valor x0 + 1. Entonces: \[log(λ_0) = β_0 + β_1x_0 − offset\] \[log(λ_1) = β_0 + β_1(x_0 + 1) − offset\] \[log (\frac{λ_1}{λ_0})= log(λ_1) − log(λ_0) = β1 ⇒ RR =\frac{λ_1}{λ_0}= e^{β_1}\]
Es decir, e β1 es el riesgo relativo entre las poblaciones donde X = x0 y X = x0 + 1, y e β0 es la tasa de incidencia cuando todas las covariables son cero. Lectura sobre RR y OR (https://data.library.virginia.edu/comparing-proportions-with-relative-risk-and-odds-ratios/)
EJEMPLO REGRESIÓN POISSON PARA TASAS DE INCIDENCIA
# Lectura de los datos
setwd("C:/Users/Alex/Desktop/MÁSTER MATEMÁTICAS/MODELIZACIÓN ESTADÍSTICA/IRANZU BARRIO")
diabetes <- read.table("diabetes.txt",header=TRUE)
names(diabetes)## [1] "edad" "sexo" "mes" "periodo" "casos" "poblacion"
## [7] "estacion"En este estudio se han planteado los siguientes objetivos:
Ver si hay cambios en la incidencia en los 7 años de registro globalmente.
Explorar si la incidencia varía en función de los grupos de edad y/o sexo.
Tenemos que calcular cual es esa exposicion sobre la que vamos a calcular la tasa.
diabetes$edad=factor(diabetes$edad)
diabetes$periodo=factor(diabetes$periodo)
diabetes$estacion=factor(diabetes$estacion)
diabetes$sexo=factor(diabetes$sexo)
diabetes$mes=factor(diabetes$mes)
diabetes[1:10,]## edad sexo mes periodo casos poblacion estacion
## 1 1 2 12 1997 2 110324 1
## 2 1 2 8 1997 1 110324 2
## 3 2 2 9 1997 0 120910 2
## 4 3 2 1 1997 1 146712 1
## 5 1 1 11 1997 1 116134 1
## 6 2 2 10 1997 3 120910 2
## 7 3 1 5 1997 3 154741 2
## 8 1 2 6 1997 0 110324 2
## 9 1 2 10 1997 0 110324 2
## 10 3 2 3 1997 3 146712 1attach(diabetes)
tasa=casos/poblacion
tasa2=tasa*100000
tasa2[edad==1 & mes==12 & sexo==2 & periodo==1997]## [1] 1.812842λ = y/E, es decir, λ = casos/poblacion. Entonces, por ejemplo, para niñas, entre 0 y 4 años, en el mes de Diciembre del año 1997, y = 2 y E = 110324, por lo tanto,λ = 2/110324 = 0, 0000181 o 1.81 cada 100000 personas.
Si calculao la tasa y me interesa por la población, lo que le tengo que meter al modelo, el offset, es el logaritmo de la exposición. Por lo tanto, la exposición en este caso serñia mi población, y tengo que calcular el logaritmo de la población, que es lo que le tengo que meter como offset al modelo.
logpobla <- log(poblacion)
diabetes1 <- glm(casos~periodo+offset(logpobla),family=poisson) # Con la opcion offset le estoy indicando que es un modelo para tasas.
summary(diabetes1)##
## Call:
## glm(formula = casos ~ periodo + offset(logpobla), family = poisson)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.0439 -0.7709 -0.3837 0.7293 2.7097
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -11.21951 0.08944 -125.438 <2e-16 ***
## periodo1998 0.02386 0.12649 0.189 0.850
## periodo1999 0.04955 0.12599 0.393 0.694
## periodo2000 0.13387 0.12391 1.080 0.280
## periodo2001 -0.15710 0.13238 -1.187 0.235
## periodo2002 0.08106 0.12369 0.655 0.512
## periodo2003 0.03106 0.12527 0.248 0.804
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 686.27 on 503 degrees of freedom
## Residual deviance: 680.44 on 497 degrees of freedom
## AIC: 1716.9
##
## Number of Fisher Scoring iterations: 5Da igual que ponga el offset en postivo o en negativo por que no deja de ser una cantidad constante.Ajusto mi variable respuesta. Aqui no tengo que hacer nada, le tengo que decir que son los casos, por que ya le estoy diciendo al meterle eloffser, que lo que estoy estimando es la tasa y no el numero de casos.Y de nuevo le digo que la familia es poisson.
Ajusto el modelo en este caso introduciendo solo cómo variable explicativa el periodo. Y voy a obtener una estimación beta ( mi modelo es el logatirmo de la tasa b198, b299, b3*2000…) hemos cogido cómo referencia un año que será el 1997, y lo que yo quiero ver es si puedo decir que todos esos betas son iguales a cero, o si existen algunos betas que son disitntos de cero( uno o más) , de forma que hagan que la tasa de incidencia de diabetes mellitus no sea la misma para todos los años.
Es decir, lo de siempre, lo unico que ahora estámos hablando de tasas de incidencia. ¿Que puedo hacer? Por un lado me fijo en los p-valores correspondientes al test de wall, y en este caso lo que veo es que ninguno es distinto con respecto al 97. Lo que no llego a ver ahi es si hay diferencias entre algunos que están ahi dentro. Como los p-valores son mayores a 0.05, lo que puedo decir claramente es que no hay diferencias entre el 98 y el 97, no hay diferencias entre 99 y el 97, entre el 2000 y el 97…¿pero las habría entre el 2000 y el 2001? Con esta salida yo no lo sé.
¿Qué tengo que hacer? El test de la razon de verosimilitud.
anova(diabetes1, test="Chisq")## Analysis of Deviance Table
##
## Model: poisson, link: log
##
## Response: casos
##
## Terms added sequentially (first to last)
##
##
## Df Deviance Resid. Df Resid. Dev Pr(>Chi)
## NULL 503 686.27
## periodo 6 5.8389 497 680.44 0.4415Obtengo un p-valor mayor a 0.05, por tanto no puedo decir que haya diferencias en la tasa de incidencia de diabetes mellitus a lo largo de los años.Esa tasa de incidenciapuedo asumir que permanece constante, no varía a lo largo de los años.
¿Qué ocurre si considero la variable sexo? si son niñas o son niños.
diabetes2 <- glm(casos~sexo+offset(logpobla),family=poisson) # Con la opcion offset le estoy indicando que es un modelo para tasas.
summary(diabetes2)$coefficients## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -11.17502599 0.04657464 -239.9380198 0.0000000
## sexo2 -0.03774752 0.06728807 -0.5609839 0.5748085Lo que veo es que no hay evidencias estadísticamente significativas en las tasas de incidencia de diabetes mellitus entre niños uy niñas. Es decir, que se puede dar la misma tasa de incidencia tanto en niños que en niñas. Y lo miro para la edad. Y de nuevo como veis en todos los modelos lo que estoy haciendo es incluir el offset del logaritmo de la población.
diabetes3 <- glm(casos~edad+offset(logpobla),family=poisson) # Con la opcion offset le estoy indicando que es un modelo para tasas.
summary(diabetes3)$coefficients## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -11.4517478 0.06804131 -168.305812 0.000000e+00
## edad2 0.3897737 0.08762647 4.448127 8.662219e-06
## edad3 0.3293682 0.08696074 3.787551 1.521397e-04Ajusto el modelo para la edad.Recordamos que es una variable categorica que toma 3 categorías con lo cual cogemos a los más jovenes cómo referencia.
Y en este caso lo que vemos es que si hay evidencias estadisticamente significativas. Tanto el grupo de 5-9 años como el de 10-15 años son estadisticamente diferentes con respecto a los mas jovenes de 0-4 años.
Cómo tenemos betas positivos, (recordad que para la interpretación lo que tenemos que hacer calcular el exponencial de esos parámetros estimados ) vamos a tener riesgos relativos superiores a 1. Cómo estamos modelizando las tasas de incidencia,la interpretación la vamos a hacer en terminos de riesgos relativos, entonces el riesgo relativo en los niños entre 5-9 años va a ser e^0.389 veces mayor que el de los niños más jovenes.
Es decir, a medida que auemtna la edad, el riesgo de que la tasa de incidencia aumente es mayor.
Como veis, yo ahi lo unico que se es que hay diferencias entre los niños de 5-9 y 0-4 y los niños de 10-15 y 0-4, pero no sé si hay diferencias entre los niños de 10-15 y 5-9. Si quisiese ver eso tendría que cambiar la categoría de referencia y comparar esas dos categorías que me interesan.
entonces: + el numero de año no influye, es decir, la tasa de incidencia se mantiene constante a lo largo de los años. + si son niños o niñas no influye, + lo que influe es la edad.
Podemos plantear, aunque hemos visto que ser niño o niña no influye, podríamos ver si el echo de introducir el sexo junto con la edad, influiría.
Muchas veces las variables que en un modelo univariante tienen un p-valor ni significativo pero que está por debajo de 0.2, se meten en el modelo multivariante, ya que los parámetros que estimamos varian en función del numero de variables que yo meto en el modelo.
hay veces que si mi variable no es significativa pero tampocoo está muy lejos de serlo, puede pasar que cuando la meto en el modelo en conbinación con otras variables, puede resultar que sea significativa. Es decir, que si yo tengo en cuenta la edad, el sexo si sea sifnificativo.
diabetes4<-glm(casos~edad+sexo+offset(logpobla), family=poisson)
anova(diabetes3, diabetes4,test='Chisq')## Analysis of Deviance Table
##
## Model 1: casos ~ edad + offset(logpobla)
## Model 2: casos ~ edad + sexo + offset(logpobla)
## Resid. Df Resid. Dev Df Deviance Pr(>Chi)
## 1 501 663.43
## 2 500 663.10 1 0.32205 0.5704En este caso no lo es. Comparamos el modelo que solo tiene la edad con el modelo que tiene la edad y el sexo, y lo que vemos es que no rechazamos la hipotesis nula, es decir, nos quedamos con el modelo mas pequeño, el modelo más simple, que en este caso, tiene la edad.
Entonces ahora lo que nos interesaría es estimar la tasa de incidencia en base a este modelo para mis individuos en la población. como hemos dicho que la variable sexo no es significativa, la tasa de incidencia estimada para niños y niñas va a ser la misma
diabetes.ajustados=predict(diabetes3,type="response") ## Numero medio de niños con diabetes diferente por cada grupo de edad.
tasa.ajustada=diabetes.ajustados/poblacion ## Tasa de diabetes
tapply(tasa.ajustada,sexo,mean) ## tasa de diabetes diferente por cada sexo.## 1 2
## 1.370227e-05 1.370227e-05Sin embargo, como el modelo si tiene en cuenta la edad, la tasa de incidencia ajustada es diferente para los 3 grupos de edad. pero como veis, la mayor diferencia se da para los 1 y 2, para los 2 y 3 las tasas de incidencia son muy parecidas.
tapply(tasa.ajustada,edad,mean) ## tasa de diabetes diferente por cada grupo de edad.## 1 2 3
## 1.063088e-05 1.569805e-05 1.477787e-05O que podrñiamos imaginar, y habria que hacerlo, es comprobar si entre esos dos grupos de edad hay diferencias estadisticamente significativas.
Con la funcion predict, considrando mi base de datos, le digo que quiero los valores estimados en termino de la variable respuesta. Entonces, mi variable respuesta sigue siendo el numero de casos.
con lo cual, si quiero ajustar la tasa, ( como estoy haciendo esto, ya he metido en el modelo eloffset, esa constante la esta teniedno en cuenta. ) lo que pasa es que como estoy ajustando el nuemro de casos de diabetes mellitus, si quiero ajustar la tasa, tengo que hacer el numero de casos dividido por la población total.