Algunos modelos de regresión
Alejandro Quintela del Rio
2021-02-17
library(tidyverse)
library(vcd)
library(freqparcoord)
library(ggridges)
library(gridExtra)
library(kableExtra)
library(effects)
library(visreg)
library(DescTools)
library(readr)
library(pscl)
library("sandwich")
library("lmtest")
library("MASS")
library("car")
library(AER)
library(tidyverse)
library(stargazer) #for combining output
library(ggplot2) #for graphing
library("mpath")
library("zic")
library(finalfit)
library(reshape2) #convert wide to tall
library(summarytools) #freq
library(jtools)
library(mosaic)
library(sjPlot)
library(gtools)
library(haven)
library(performance)
library(radiant)
library(readxl)
library(correlation)
library(see) # for plotting
library(ggraph) # needs to be loaded
## convenience functions for visualization
clog <- function(x) log(x + 0.5)
cfac <- function(x, breaks = NULL) {
if(is.null(breaks)) breaks <- unique(quantile(x, 0:10/10))
x <- cut(x, breaks, include.lowest = TRUE, right = FALSE)
levels(x) <- paste(breaks[-length(breaks)], ifelse(diff(breaks) > 1,
c(paste("-", breaks[-c(1, length(breaks))] - 1, sep = ""), "+"), ""), sep = "")
return(x)
}
######## esto solo es para que no se escriban los warnings
oldw <- getOption("warn")
options(warn = -1)
########################################
# leemos archivo de datos
Bdatos <- read_excel("BaseDatos28FebreroSoloVariablesVisiblesxlsx.xlsx")
Bdatos <- Bdatos %>% mutate_if(is.character,as.factor)
Bdatos$Caidas=Bdatos$CAIDAS.ANO.ANT
Bdatos <-(Bdatos %>% filter(Bdatos$Caidas<38) )
Bdatos <- Bdatos %>%
mutate(trozos= case_when(
.$Caidas >= 0 & .$Caidas <4 ~ "0",
.$Caidas >= 4 ~ "1"))
Bdatos$trozos <- as.factor(Bdatos$trozos)
Bdatos <- Bdatos %>%
mutate(fac.caidas = case_when(
.$Caidas >= 0 & .$Caidas <1 ~ "0",
.$Caidas >= 1 & .$Caidas < 2 ~ "1",
.$Caidas >= 2 & .$Caidas < 4 ~ "2-3" ,
.$Caidas >= 4 & .$Caidas <= 5 ~ "4-5",
.$Caidas >= 6 & .$Caidas <= 10 ~ "6-10" ,
.$Caidas > 10 ~ "10+" )
)
ord_fac.caidas <- c("0","1","2-3","4-5","6-10","10+")
Bdatos <- Bdatos %>%
mutate(fac.caidas=factor(fac.caidas, levels=ord_fac.caidas))
Bdatos$caidas_log=clog(Bdatos$Caidas)
Bdatos$fac.edad=cfac(Bdatos$EDAD, breaks =c(20,40,60,80))
Bdatos$fac.medic=cfac(Bdatos$NUM.MEDIC, breaks= c(0,2,5,20))
Bdatos$fac.marcha=quantcut(Bdatos$SUMA.ALT.MARCHA, 3)
Bdatos$fac.obs=quantcut(Bdatos$SUMA.DIFICUL.OBSTACULOS, 3)
Bdatos$fac.tineti=cfac(Bdatos$TINETTI.TOTAL, breaks =c(10,19,24,28))
Bdatos$fac.grup.medic=quantcut(Bdatos$CAM.MED.12M, 3)
ord_GRADO.DI <- c("Leve","Moderado","Grave","Profundo")
Bdatos <- Bdatos %>%
mutate(GRADO.DI=factor(GRADO.DI, levels=ord_GRADO.DI))
ord_IMC <- c("Infrapeso","Normopeso","Sobrepeso","Obesidad I",
"Obesidad II", "Obesidad III")
Bdatos <- Bdatos %>%
mutate(IMC=factor(IMC, levels=ord_IMC))
attach(Bdatos)1
2 Regresión lineal simple
La regresión (de manera general, considerada en \(p+1\) dimensiones) estudia como influyen un conjunto de \(p\) variables \(X_1,X_2,..,X_p\) (en principio, independientes) en una variable dependiente \(Y\).
Cuando sólo se considera la relación entre \(Y\) y una única variable \(X\) suele hablarse de regresión lineal simple.
La ecuación de una regresión lineal simple (línea de regresión) es:
\[ y = b_0 + b_1x + \varepsilon \] El término \(\varepsilon\) es el error que aparece debido a que no todos los puntos de un diagrama de dispersión están sobre la recta (si todos los puntos están sobre la recta, la correlación lineal va a ser \(1\), no hay ningún error). Claramente, para un valor \(X\) concreto, \(Y\) puede tomar diferentes valores, y la regresión “promedia” esos valores.
De manera general, la regresión se utiliza cuando las variables \(X\) e \(Y\) son continuas (o discretas que tomen muchos valores), y los errores \(\varepsilon\) son errores que siguen una distribución normal (gaussiana) de varianza constante. ¿Esto qué quiere decir? Pues que, para un valor de \(X\) cualquiera, dado que \(Y\) no toma ún unico valor, los distintos puntos \(Y\) tienen que seguir una distribución normal, como se ve en esta gráfica:
Los errores son las distancias de los diferentes valores de \(Y\) al valor “promedio” que sería el punto que está en la recta.
Vamos a hacer un ejemplo con el conjunto de datos mlb, que corresponde a las estaturas y pesos de jugadores de la liga americana de beisbol.
data(mlb)
mlb$Estatura=mlb$Height # está en pulgadas
mlb$Peso=mlb$Weight # en libras
mlb$Edad=mlb$Age
ggplot(mlb, aes(x=Estatura, y=Peso))+
geom_point(alpha=.2)+
ggtitle("Estaturas y pesos de jugadores de la liga MLB") +
geom_smooth(method=lm, se=FALSE) 
Vamos a centrarnos en un valor particular para la estatura, por ejemplo \(75\) pulgadas:
ggplot(mlb, aes(x=Estatura, y=Peso))+
geom_point(alpha=.2)+
ggtitle("Estaturas y pesos de jugadores de la liga MLB") +
geom_smooth(method=lm, se=FALSE) +
geom_vline(xintercept = 75, col="red")
Analizamos el conjunto de datos que están solo en esa línea de color rojo:
x= mlb$Peso[mlb$Estatura==75]
summary(x)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 156.0 197.0 210.0 208.7 220.0 278.0op <- par(mfrow = c(1, 2)) # 2 gráficas en 1 linea, 2 columnas
boxplot(x, horizontal = T)
plot(density(x))
rug(x) 
par(op) # devolvemos los gráficos a 1 linea, 1 columnaAtendiendo al sumario de los datos, al boxplot y a la estimación de la densidad, parece que los datos de esa línea vertical cumplen la hipótesis de pertenecer a una distribución gaussiana. Lo único que podría perjudicar esa suposición serían un par de valores excesivamente altos que se observan a la derecha. En todo caso, si modificamos ligeramente el valor del párametro “bandwidth” (ancho de los intervalos) del estimador de la densidad
plot(density(x, bw= "nrd0",adjust=2) ) 
observamos que puede existir una buena aproximación a la campana de Gauss.
Es posible realizar una gráfica que refleje esta situación para cada valor del eje de las \(X\), mediante (viene a ser la gráfica de las estaturas y los pesos que está arriba, pero girada noventa grados) la siguiente figura. En ella se observa que, para cada valor de la estatura, los pesos de dicha estatura siguen una distribución gaussiana (aproximadamente, salvo casos extremos).
ggplot(mlb, aes(y=Estatura, x=Peso, group=Estatura))+
geom_density_ridges(jittered_points = TRUE,
position = position_points_jitter(width = 0.5, height = 0),
point_size = 1,
alpha = 0.7)+
ggtitle("Densidad de los Pesos de los jugadores para cada Estatura")## Picking joint bandwidth of 6.64
2.1 Supuestos teóricos
Los supuestos teóricos en los que se puede utilizar la regresión se basan en que todas estas distribuciones de valores (para cada valor de la variable \(X\)) sean normales (se aproximen a la campana de Gauss), de media cero y tengan varianza constante (la variabilidad sea siempre la misma, aproximadamente). Aquí, por ejemplo, observamos que el rango de valores en cada recta es casi siempre el mismo, excepto en las ocasiones en que aparecen valores atípicos (que quizá habría que considerar con más detalle, pero, si nos olvidamos de ellos, siempre habría una variabilidad más o menos equivalente).
Estas condiciones no son estrictas, es decir, en ocasiones se pueden usar las mismas técnicas aunque no se verifique la normalidad de los residuos, o aunque no tengan varianza constante.
2.2 Aplicación a los datos de caidas
2.2.1 Correlaciones
Con la base de datos de Caidas, vamos a estudiar la relación entre las variables Tinettiy Tug.media.
plot(TINETTI.TOTAL, TUG.MEDIA, main="Nube de puntos")
Calculamos la correlación lineal: \[r_{xy} = \frac{s_{xy}}{s_x s_y}\]
cor(TINETTI.TOTAL, TUG.MEDIA)## [1] -0.70600982.2.2 Cálculo de la línea de regresión
La función que se utiliza en R es lm
2.2.3 Modelo
y=TUG.MEDIA
x=TINETTI.TOTAL
linearModel <- lm(y ~ x, data=Bdatos)
linearModel ##
## Call:
## lm(formula = y ~ x, data = Bdatos)
##
## Coefficients:
## (Intercept) x
## 35.969 -1.021La recta es \(y=35.969-1.021\cdot x\).
2.2.4 Información del modelo
Se obtiene mediante la función summary.
summary(linearModel)##
## Call:
## lm(formula = y ~ x, data = Bdatos)
##
## Residuals:
## Min 1Q Median 3Q Max
## -14.729 -2.491 -1.107 1.675 17.766
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 35.96857 2.42517 14.83 <2e-16 ***
## x -1.02068 0.09807 -10.41 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 4.658 on 109 degrees of freedom
## Multiple R-squared: 0.4984, Adjusted R-squared: 0.4938
## F-statistic: 108.3 on 1 and 109 DF, p-value: < 2.2e-16El análisis detallado del sumario de un modelo lineal lo vemos en el apartado siguiente de la regresión múltiple.
2.2.5 Gráfica de la línea de regresión sobre la nube de puntos
Así de sencillo:
plot(x, y, xlab="Tinetti", ylab="Tug")
abline(linearModel)
2.2.6 Interpretación de los coeficientes de la recta
\(b_0\) es el término intercept(intersección de la recta con el eje vertical), puesto que es el valor que toma \(y\) cuando \(X\) vale cero. En este caso no da información alguna (recordemos que las conclusiones deben ceñirse al rango de datos en los que estemos trabajando. En este caso, por ejemplo, no hay valores de Tinetti menores que \(10\)).
\(b_1\) es la pendiente de la recta. Aquí es un valor negativo, lo que significa que la recta es decreciente. ¿Tiene algún significado el valor concreto? Supongamos que la variable \(X\) crece una unidad. \(X\) pasa de valer \(x_0\) a \(x_0 +1\). Entonces \(Y\) pasa de valer \(b_0 + b_1\cdot x_0\) a valer \(b_0 + b_1\cdot (x_0+1)\), es decir que \(Y\) crece \(b_0 + b_1\cdot (x_0+1) - \left( b_0 + b_1\cdot x_0 \right) = b_1\).
Vemos que el coeficiente \(b_1\) es lo que cambia la variable \(Y\) cuando \(X\) crece una unidad.
En nuestro caso, si Tinetti crece 1 unidad, Tug decrece \(1.02\) (en media, puesto que hay que considerar el error de los puntos que no están sobre la recta)
3 Regresión lineal múltiple
Es una extensión del modelo anterior, cuando consideramos que la variable \(Y\) depende de más variables \(X_1,X_2,...,X_p\).
\[Y = b_0 + b_1X_1 + b_2X_2 + ...b_pX_p + error\]
Volvamos a considerar el conjunto de datos mblde los jugadores de beisbol. El comando lm que hay que utilizar es el mismo de la regresión lineal simple, pero añadiendo más variables explicativas. Ejemplo:
regmul <- lm(Peso ~ Estatura + Edad, data=mlb)
summary(regmul)##
## Call:
## lm(formula = Peso ~ Estatura + Edad, data = mlb)
##
## Residuals:
## Min 1Q Median 3Q Max
## -50.535 -12.297 -0.297 10.824 74.300
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -187.6382 17.9447 -10.46 < 2e-16 ***
## Estatura 4.9236 0.2344 21.00 < 2e-16 ***
## Edad 0.9115 0.1257 7.25 8.25e-13 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 17.21 on 1012 degrees of freedom
## Multiple R-squared: 0.318, Adjusted R-squared: 0.3166
## F-statistic: 235.9 on 2 and 1012 DF, p-value: < 2.2e-16Analizamos ahora con detalle los resultados de summary:
El primer elemento (call) solo recuerda la fórmula que se ha utilizado (variable dependiente \(Y\), en este caso el Peso, frente a la Estatura y la Edad).
El sumario de los residuos muestra cómo se distribuyen los residuos (errores) del modelo. Los residuos son la diferencia entre los puntos reales y la predicción del modelo, y si los residuos se distribuyen simétricamente alrededor de 0 es una buena señal, porque la regresión lineal está asumiendo que estos errores siguen una distribución normal (como hemos visto en la regresión lineal simple). Podemos dibujarlos y también hacer una estimación de su densidad:
plot(regmul$residuals, main="Residuals Plot")
abline(h=0)
plot(density(regmul$residuals))
El primer gráfico muestra que los residuos son variables (positivos o negativos), con media en cero y varianza (desviación típica) constante (los valores más alejados del cero, tanto por arriba como por abajo no varían a lo largo de todo el gráfico. Si hubiera intervalos donde eso no se cumpliese, los residuos tendrían varianza no constante). Ejemplo:
Ejemplo de residuos con varianza no constante
Analicemos los componentes que aparecen en el término Coefficients:
El primero es la estimación de cada valor \(b_0\), \(b_1\) etc. de la ecuación de regresión.
En este caso, la estimación de \(b_1\) representa que, por cada unidad de aumento de la Estatura, el Peso crece en \(4.92\) unidades, siempre que la Edad permanezca constante (y si hubiera más variables consideradas, tendriamos que suponer que también permanecen constantes).
Si la estatura permanece constante, un aumento de la edad en una unidad hace que el peso crezca \(0.91\) unidades.
El error estándar estima la cantidad promedio que las estimaciones varían de la respuesta real (es la varianza de cada estimador \(b_i\) dividido entre el tamaño de la muestra). Intuitivamente, interesa que ese error sea lo más pequeño posible, porque nos estará diciendo que nuestra estimación es buena.
El \(t\ value\) de cada coeficiente representa cuán lejos (cuántas desviaciones estándar) está la estimación de \(0\). Para cada uno de los coeficientes, se está realizando un contraste de hipótesis, para ver si los parámetros \(b_0,b_1,...,b_p\) son cero o no lo son. El \(p-valor\) de esta prueba es \(Pr(>|t|)\), de manera que nos aparecerán asteriscos si el \(p-valor\) es más pequeño que \(0.05\). Como vemos, en este caso todos los términos estimados (\(b_0,b_1\) y \(b_2\)) son significativos.
El número \(R^2\) (multiple R-squared) mide la precisión del ajuste de la curva a los puntos (es la extensión del \(R^2\) de la regresión lineal simple). El número en el que nos debemos fijar es, en realidad, el Adjusted R-squared, que es un \(R^2\) corregido por el número de datos y variables. Este número, multiplicado por cien, nos da la precisión en tanto por ciento (\(31.66\) en este caso).
Es importante comentar que, siempre que se aumente el número de variables en un modelo de regresión lineal, el \(R^2\) queda igual o crece, nunca disminuye. Es por eso que deben seleccionarse bien las variables que se introducen en el modelo, puesto que utilizar el \(R^2\) para ello puede resultar muy engañoso.
Para los datos de caídas, vamos a plantear un modelo con el número de caidas en función de la edad, la medicación y el valor de Tinetti.total.
modelo1<-lm(CAIDAS.ANO.ANT ~ TINETTI.TOTAL + EDAD + NUM.MEDIC,
data=Bdatos)
summary(modelo1)##
## Call:
## lm(formula = CAIDAS.ANO.ANT ~ TINETTI.TOTAL + EDAD + NUM.MEDIC,
## data = Bdatos)
##
## Residuals:
## Min 1Q Median 3Q Max
## -8.435 -1.791 -0.411 1.042 16.007
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 11.17086 2.99500 3.730 0.000308 ***
## TINETTI.TOTAL -0.52221 0.09051 -5.770 7.81e-08 ***
## EDAD 0.07165 0.03877 1.848 0.067333 .
## NUM.MEDIC 0.30925 0.14151 2.185 0.031042 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 4.007 on 107 degrees of freedom
## Multiple R-squared: 0.3842, Adjusted R-squared: 0.3669
## F-statistic: 22.25 on 3 and 107 DF, p-value: 2.836e-11Si nos fijamos en los residuos, vemos que la mediana no es cero, y el primer cuartil no guarda simetría (con el cero como eje de simetría) con el tercer cuartil, ni el mínimo valor guarda simetría con el máximo. Esto es una primera indicación de que, seguramente, la hipótesis de residuos gaussianos no se cumple. Vamos a hacer un gráfico de los residuos, junto con un dibujo de la estimación de la densidad:
op <- par(mfrow = c(1, 2))
plot(modelo1$residuals, main="Residuals Plot")
abline(h=0)
plot(density(modelo1$residuals))
par(op)En el primer gráfico comprobamos que no existe simetría alrededor del cero, y, además, la varianza de los residuos o errores no es constante: vemos, por ejemplo que, en el intervalo \((20,60)\) los residuos llegan a tomar valores entre \(10\) y \(15\), mientras que en el intervalo de \(80\) en adelante no pasan de \(5\). La estimación de la densidad también refleja asimetría a la derecha.
En valor \(R^2\) ajustado es de \(0.36\). En un modelo de regresión lineal múltiple, se considera, de manera general, que tenemos un buen ajuste si \(R^2\) es mayor de \(0.7\) o \(0.8\).
4 Regresión logística
En el caso de que la variable \(Y\) (dependiente) únicamente tome valores \(0\) y \(1\) no tiene sentido plantear una relación lineal (simple o múltiple). Con el siguiente gráfico vemos que ocurre si tenemos datos bidimensionales donde la \(Y\) sea binaria:
sim_logistic_data = function(sample_size = 100, beta_0 = -2, beta_1 = 3) {
x = rnorm(n = sample_size)
eta = beta_0 + beta_1 * x
p = 1 / (1 + exp(-eta))
y = rbinom(n = sample_size, size = 1, prob = p)
data.frame(y, x)
}
set.seed(1)
example_data = sim_logistic_data()
# ordinary linear regression
fit_lm = lm(y ~ x, data = example_data)
# logistic regression
fit_glm = glm(y ~ x, data = example_data, family = binomial)
plot(y ~ x, data = example_data,
pch = 20, ylab = "Probabilidad",
main = "Regresion lineal vs Logistica")
grid()
abline(fit_lm, col = "darkorange")
curve(predict(fit_glm, data.frame(x), type = "response"),
add = TRUE, col = "dodgerblue", lty = 2)
legend("topleft", c("Ordinary", "Logistic", "Data"), lty = c(1, 2, 0),
pch = c(NA, NA, 20), lwd = 2, col = c("darkorange", "dodgerblue", "black"))
Podemos observar que cualquier regresión lineal va a hacer un ajuste extremadamente malo. Mientras tanto, observamos que puede trazarse una función no lineal (la azul, llamada logística) que va a estimar la probabilidad de que \(Y\) tome el valor \(0\) o el valor \(1\).
Llamemos \(p\) a la probabilidad de que la variable \(Y\) tome el valor \(1\). \(Y\) toma entonces el valor \(0\) con probabilidad \(1-p\).
La regresión logística consiste en ajustar un modelo de regresión lineal de la siguiente forma: \[ log\left(\dfrac{p}{1-p}\right) = b_0 + b_1X_1 +b_2X_2 +...b_pX_p \]
Fijémonos que el término de la izquierda es el logaritmo delodds o razón entre la probabilidad de que \(Y\) valga \(1\) contra que valga \(0\).
Si quitamos el logaritmo en la última ecuación, quedaría: \[ \left(\dfrac{p}{1-p}\right) = e^{b_0 + b_1X_1 +b_2X_2 +...b_pX_p} \] y, despejando \(p\) \[ p=\dfrac{e^{b_0 + b_1X_1 +b_2X_2 +...b_pX_p}}{1+e^{b_0 + b_1X_1 +b_2X_2 +...b_pX_p}} \]
La función \(f(z)=\dfrac{e^z}{1+e^z}\) se llama función logística, y de ahí el nombre de la regresión.
4.1 Interpretación de los coeficientes
Ahora los coeficientes \(b_1,b_2,...,b_p\) tienen otra interpretación, porque lo que se expresa linealmente en función de las variables \(X_j\) no es directamente la variable \(Y\), sino el logaritmo del odds.
Igual que hicimos en la regresión lineal múltiple, vamos a ver como cambia la variable dependiente cuando una variable \(X_j\) cualquiera aumenta una unidad, es decir pasa de valer \(x_0\) a valer \(x_0+1\). Estamos pasando de un estado donde la probabilidad de \(Y=1\) la llamamos \(p_0\), a un estado donde llamamos \(p_1\) a la misma probabilidad.
\[ log\left(\dfrac{p_{0}}{1-p_{0}}\right) = b_0+b_1X_1+...+b_jx_0+...b_pX_p \]
Si pasamos del valor \(x_0\) a \(x_0+1\) (y el resto de variables permanece constante) \[ log\left(\dfrac{p_{1}}{1-p_{1}}\right) = b_0+b_1X_1+...+b_j(x_0+1)+...b_pX_p \] Restamos la primera ecuación a la segunda: \[ log\left(\dfrac{p_{1}}{1-p_{1}}\right) - log\left(\dfrac{p_{0}}{1-p_{0}}\right) = b_j, \] esto es, \[ log\left(\dfrac{p_{1}/(1-p_1)}{p_0/(1-p_0)}\right) = b_j. \]
Como vemos, el coeficiente \(b_j\) es el logaritmo del odds-ratio entre los momentos \(x_0+1\) y \(x_0\), es decir \(e^{b_j}\) es el Odds-Ratio.
Es precisamente este análisis el que permite, por un lado, incluir variables binarias o categóricas en este modelo. Una variable podría ser, por ejemplo, el sexo, o estar enfermo o no. El coeficiente \(e^{b_j}\) significa entonces el odds-ratio entre un estado frente al otro.
En ocasiones, se indica esta proporción en tanto por ciento, es decir, si pasamos de un valor \(z\) a un valor \(z+\frac{\alpha}{100}\):
\[\dfrac{z+\frac{\alpha }{100} \cdot z}{z} =e^{b_j} \ \rightarrow \dfrac{z\cdot \left(1+\frac{\alpha }{100}\right)}{z}=e^{b_j} \ \rightarrow 1+\frac{\alpha }{100} =e^{b_j} \]
y entonces \[ \alpha = \left(e^{b_j}-1\right)\cdot 100 \]
4.2 Ejemplo con variables categóricas
Con nuestra base de datos hemos considerado dos grupos: los que caen más de lo normal cuando caen más de 3 veces por año (al ser esa la media) y los que caen dentro de lo normal. Esta variable la llamamos CAE.
Vamos a realizar una regresión logística entre esta variable y las seleccionadas como influyentes. En R se utiliza el comando glm (generalized linear model).
Bdatos$CAE=as.logical(Bdatos$trozos=="1")
cae_glm<-glm(CAE ~EDAD + TINETTI.TOTAL+ SUMA.ALT.MARCHA + NUM.MEDIC+
GRADO.DI+ DEF.ATENCI, family="binomial", data=Bdatos)
summary(cae_glm)##
## Call:
## glm(formula = CAE ~ EDAD + TINETTI.TOTAL + SUMA.ALT.MARCHA +
## NUM.MEDIC + GRADO.DI + DEF.ATENCI, family = "binomial", data = Bdatos)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.21800 -0.45806 -0.21549 -0.07331 3.08998
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.46587 3.79522 0.123 0.9023
## EDAD 0.08296 0.03882 2.137 0.0326 *
## TINETTI.TOTAL -0.26952 0.12097 -2.228 0.0259 *
## SUMA.ALT.MARCHA -0.10585 0.29731 -0.356 0.7218
## NUM.MEDIC 0.21066 0.12699 1.659 0.0971 .
## GRADO.DIModerado -2.10117 1.01709 -2.066 0.0388 *
## GRADO.DIGrave 0.34572 0.85729 0.403 0.6867
## GRADO.DIProfundo -0.53564 1.60084 -0.335 0.7379
## DEF.ATENCISi 1.50330 0.74545 2.017 0.0437 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 120.844 on 110 degrees of freedom
## Residual deviance: 63.028 on 102 degrees of freedom
## AIC: 81.028
##
## Number of Fisher Scoring iterations: 6Las variables que aparecen como significativas son la EDAD, TINETTI.TOTAL, NUM.MEDIC (no al nivel \(0.05\)), el Deficit de atención (el valor “si”) y la categoria “Moderado” del grado de discapacidad.
El que otras variables no aparezcan como significativas supone que no son de relevancia a la hora de predecir el cambio de un grupo (“cae poco” a “cae mucho”), cuando se consideran dentro del conjunto total de variables explicativas.
Como dijimos antes, ahora los coeficientes de la regresión logística presentan una explicación, pero no por sí mismos, sino al considerar el número \(e\) elevado a estos. Se puede hacer directamente (y obteniendo también los intervalos de confianza)
cae_glm$coefficients## (Intercept) EDAD TINETTI.TOTAL SUMA.ALT.MARCHA
## 0.46587287 0.08295639 -0.26951847 -0.10584926
## NUM.MEDIC GRADO.DIModerado GRADO.DIGrave GRADO.DIProfundo
## 0.21066275 -2.10116863 0.34572479 -0.53563922
## DEF.ATENCISi
## 1.50329932coef_log<- cae_glm$coefficients
exp(cbind(OR = coef_log, confint(cae_glm)))## Waiting for profiling to be done...## OR 2.5 % 97.5 %
## (Intercept) 1.5934044 0.0008926378 3388.6349951
## EDAD 1.0864944 1.0118943226 1.1811646
## TINETTI.TOTAL 0.7637472 0.5893627678 0.9562584
## SUMA.ALT.MARCHA 0.8995602 0.4851461556 1.6070419
## NUM.MEDIC 1.2344959 0.9698475637 1.6059839
## GRADO.DIModerado 0.1223134 0.0125808395 0.7592801
## GRADO.DIGrave 1.4130137 0.2584402947 7.9117937
## GRADO.DIProfundo 0.5852950 0.0205812933 12.1612813
## DEF.ATENCISi 4.4965000 1.0607278470 20.6213662La edad tiene como OR un valor de \(1.08\), lo que significa que, por cada año que aumenta, el riesgo de entrar en el grupo de “Caedores” aumenta \(1.08\).
Por cada unidad que aumenta el valor de TINETTI.TOTAL, el riesgo crece \(0.76\) (disminuye).
En las variables categóricas se razona igual pero pasando de una categoría a la otra. Vemos que DEF.ATENCIsi es \(5.69\), entonces, al pasar de la categoria DEF.ATENCI=no a sí el Odds-Ratio crece \(5.69\).
Podemos observar estas relaciones mediante lo que se llama gráfico de efectos:
eff <- allEffects( cae_glm)
plot(eff , type = "response")
Por último, dentro de la opción summary nos aparecen unos valores llamados devianza (deviance) y AIC. Estos valores sirven, en cierto modo, para calcular el ajuste del modelo a los datos (parecido al \(R^2\)).
El término Null deviance es una medida de la variabilidad cuando se considera la variable \(Y\) sin considerar las variables \(X\), de manera que el término residual deviance mide lo mismo pero considerando todas las variables \(X\). Fijémonos que el modelo tiene \(8\) variables. El término Null deviance vale \(120.84\) con \(110\) grados de libertad (corresponde al número de datos menos \(1\)). El término residual deviance vale \(83.02\) y tiene \(102\) grados de libertad (correspondiente a \(110-8\) variables consideradas). La comparación entre el primer término y el segundo nos dice que la devianza disminuye sensiblemente al considerar las variables en el modelo (baja un \(68\) por ciento). El problema es que no tenemos, como en el caso de la correlación lineal, ningún valor de referencia (recordemos que \(R^2\) varía entre \(0\) y \(1\)). Los números que aparecen aquí dependen de las unidades de medida, y sólo nos van a servir para comparación con otros modelos.
Igual ocurre con el valor AIC (Akaike’s information criterion), que sirve para comparar entre dos o más modelos, siendo preferible el que tenga un valor más bajo.
Vamos ahora a comparar este modelo con otro donde solo introducimos las variables que nos aparecieron antes como significativas (probamos a introducir también NUM.MEDIC aunque su significatividad era baja (el \(p-valor\) era \(0.09\)):
cae_glm<-glm(CAE ~EDAD + TINETTI.TOTAL + NUM.MEDIC+
GRADO.DI+ DEF.ATENCI, family="binomial", data=Bdatos)
summary(cae_glm)##
## Call:
## glm(formula = CAE ~ EDAD + TINETTI.TOTAL + NUM.MEDIC + GRADO.DI +
## DEF.ATENCI, family = "binomial", data = Bdatos)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.17332 -0.43761 -0.22221 -0.08025 3.04496
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.37350 2.96812 -0.126 0.89986
## EDAD 0.08081 0.03820 2.116 0.03439 *
## TINETTI.TOTAL -0.23890 0.08392 -2.847 0.00442 **
## NUM.MEDIC 0.20629 0.12547 1.644 0.10015
## GRADO.DIModerado -2.04918 0.99613 -2.057 0.03967 *
## GRADO.DIGrave 0.31401 0.84938 0.370 0.71161
## GRADO.DIProfundo -0.60886 1.60010 -0.381 0.70356
## DEF.ATENCISi 1.49076 0.74371 2.004 0.04502 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 120.844 on 110 degrees of freedom
## Residual deviance: 63.156 on 103 degrees of freedom
## AIC: 79.156
##
## Number of Fisher Scoring iterations: 6coef_log<- cae_glm$coefficients
exp(cbind(OR = coef_log, confint(cae_glm)))## Waiting for profiling to be done...## OR 2.5 % 97.5 %
## (Intercept) 0.6883233 0.001792805 243.9201819
## EDAD 1.0841660 1.010787425 1.1768066
## TINETTI.TOTAL 0.7874938 0.656136694 0.9163467
## NUM.MEDIC 1.2291052 0.967799946 1.5927047
## GRADO.DIModerado 0.1288408 0.014046111 0.7782740
## GRADO.DIGrave 1.3689077 0.252983612 7.4879961
## GRADO.DIProfundo 0.5439694 0.019497443 11.3083292
## DEF.ATENCISi 4.4404579 1.049991040 20.2712927En comparación con el modelo anterior, el valor de AIC es prácticamente igual al anterior, aunque baje 2 unidades.
El gráfico de los efectos:
eff <- allEffects( cae_glm)
plot(eff , type = "response")
Por último, planteemos un modelo sin la variable NUM.MEDIC:
cae_glm<-glm(CAE ~EDAD + TINETTI.TOTAL +
GRADO.DI+ DEF.ATENCI,family="binomial", data=Bdatos)
summary(cae_glm)##
## Call:
## glm(formula = CAE ~ EDAD + TINETTI.TOTAL + GRADO.DI + DEF.ATENCI,
## family = "binomial", data = Bdatos)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.32216 -0.49093 -0.21398 -0.06459 2.86685
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.06743 2.84872 -0.024 0.98112
## EDAD 0.10185 0.03626 2.809 0.00498 **
## TINETTI.TOTAL -0.26302 0.08083 -3.254 0.00114 **
## GRADO.DIModerado -2.01632 0.95516 -2.111 0.03477 *
## GRADO.DIGrave 0.35026 0.82622 0.424 0.67162
## GRADO.DIProfundo -0.88878 1.62007 -0.549 0.58327
## DEF.ATENCISi 1.73940 0.72021 2.415 0.01573 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 120.844 on 110 degrees of freedom
## Residual deviance: 66.003 on 104 degrees of freedom
## AIC: 80.003
##
## Number of Fisher Scoring iterations: 6# coef_log<- cae_glm$coefficients
# exp(cbind(OR = coef_log, confint(cae_glm)))
# eff <- allEffects( cae_glm)
# plot(eff , type = "response")Como vemos, los valores de la devianza residual y el AIC son prácticamente iguales, aunque suben algo con respecto al modelo anterior. Puede ocurrir que la variable NUM.MEDIC sea realmente significativa, pero el \(p-valor\) no valga menos de \(0.05\) por no tener suficiente tamaño muestral.
5 Regresión de Poisson
Dado que el número de caídas es una variable de contar (discreta con valores no negativos), el modelo que debe utilizarse para tratar de explicar esta variable es el llamado modelo de regresión de contar o de Poisson.
El origen del nombre procede de la variable aleatoria de Poisson.
Un proceso de Poisson es un experimento aleatorio donde se observa la aparición de un suceso concreto (éxito) sobre un soporte continuo (generalmente el tiempo). Además, debe cumplirse que los sucesos ocurren de forma independiente y con media estable (el número medio de sucesos por unidad de medida es constante).
Ejemplos de procesos de Poisson son: clientes que acuden a una tienda por unidad de tiempo, llamadas por unidad de tiempo a una centralita, defectos por metro de cable, baches por kilometro de autopista, caídas a lo largo de un periodo de tiempo…
En un proceso de Poisson, la variable \(X\)=número de éxitos en un intervalo se dice que sigue una distribución de Poisson de parámetro \(\lambda .\) Se escribe \(X\in Pois(\lambda ).\)
Su distribución de probabilidad es \[ P(X=k)=e^{-\lambda }\frac{\lambda ^{k}}{k!},\ \ \ k=0,1,2,... \]
Se verifica que la media (o esperanza matemática) de la variable es el parámetro \(\lambda\), y también es la varianza. \[ E(X)=Var(X)=\lambda. \]
La variable aleatoria de Poisson nos va a servir para comprobar si los sucesos considerados son aleatorios o no lo son. En una situación normal, el número de caídas (o accidentes en general) a lo largo de un año deberían de ser muy pocos (salvo personas con actividades de riesgo u otras peculiaridades). En nuestro caso,
fr <- table(Bdatos$CAIDAS.ANO.ANT) %>% data.frame
names(fr) <- c('caidas', 'freq')
fr$caidas <- as.numeric(as.character(fr$caidas)) #convert factor to numeric
ggplot(fr, aes(x = caidas, y = freq)) +
geom_col() +
theme_bw() +
lims(y = c(0, 50)) +
geom_line() +
labs(x = "Numero de caidas", y = "Frecuencia") +
geom_text(aes(x = caidas, y = freq, label = freq, vjust = -1)) +
theme(axis.title.y = element_text(angle = 0)) 
El número medio de caídas en el año \(2018\) es 2.9009009 que es aproximadamente \(3\). Si comparamos gráficamente la distribución de caidas con la que correspondería si siguiera una distribución de Poisson:
library(vcd)
gf <- goodfit(Bdatos$CAIDAS.ANO.ANT, "poisson" ,par = list(lambda = 3))
plot(gf, type = "standing", scale = "raw")
Otro dibujo (me parece que queda mejor):
x = 0:max(Bdatos$Caidas)
Poisson = dpois(x,lambda=3)
dt<-data.frame(x,Poisson)
obs <- table(Bdatos$Caidas) %>% prop.table() %>% data.frame #Observed
names(obs) <- c("x", 'Observed')
obs$x<-as.numeric(as.character(obs$x))
N_total=sum(table(Bdatos$Caidas))
mm <- dplyr::inner_join(obs, dt, by = 'x' )
mm <- mutate(mm, Observados=Observed*N_total)
mm <- mutate(mm, Teoricos=round(Poisson*N_total))
ggplot(data=mm)+
geom_bar(aes(x=x,y=Observados), stat="identity",
fill = "lightblue" , width=0.8 )+
geom_bar(aes(x=x,y=Teoricos), stat="identity",
fill = "grey" , width=0.4 )+
labs(x = "Caidas en 2018", y = 'Individuos ',
title = "Modelo para caidas por individuo",
subtitle="Azul: real. Gris: Poisson. ") 
Como vemos en este gráfico, si la variable número de caídas por año siguiera una distribución de Poisson (en gris), el número máximo de caídas por año, en una situación de aleatoriedad, debería ser de \(2\) o \(3\) como máximo, y luego las probabilidades de tener más caídas de \(3\) van descendiendo rápidamente.
Por otro lado, la varianza de la variable es \(25\). Como se observa también en el gráfico, la mayoría de las personas estudiadas tienen \(0\) caídas, pero también ocurre que hay casos extremadamente altos (personas con \(10, 11\), incluso con \(25\) caídas). Esta variabilidad tan alta se conoce como sobredispersión (overdispersion). Recordemos que una variable de Poisson tiene su media igual a su varianza.
Todo esto nos índica que la distribución de Poisson, en principio, no va a ser la más adecuada para modelar el número de caídas, y habrá que probar otro tipo de modelos dentro de la regresión de contar.
Vamos a crear un data.frame con las variables que usaremos, por comodidad.
attach(Bdatos)
rc=data.frame( Caidas, EDAD, NUM.MEDIC ,
SUMA.ALT.MARCHA,
GRADO.DI,
TINETTI.TOTAL,
DEF.ATENCI)El modelo de regresión de contar intenta predecir el logaritmo del parámetro \(\lambda\) (número medio) en función de las variables explicativas que se consideren (continuas o categóricas, o ambas):
\[ log(\lambda) = b_0 + b_1X_1+b_2X_2+...+b_pX_p \]
Para calcular el modelo de regresión de poisson se escribe:
fm_pois <- glm(Caidas ~ ., data = rc, family = poisson)Dos observaciones:
El comando
glmfue el mismo que utilizamos en el modelo de regresión logística, pero antes usamosfamily=binomial.Al escribir
Caidas ~ .estamos diciendo que queremos usar como variables predictoras o independientes todas las del conjuntorc(por eso hicimos un data.frame con pocas variables).
Veamos los resultados:
summary(fm_pois)##
## Call:
## glm(formula = Caidas ~ ., family = poisson, data = rc)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -3.756 -1.529 -1.048 0.524 5.392
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.968543 0.634930 3.100 0.001933 **
## EDAD 0.034567 0.006853 5.044 4.56e-07 ***
## NUM.MEDIC 0.037205 0.019392 1.919 0.055034 .
## SUMA.ALT.MARCHA -0.070213 0.049192 -1.427 0.153484
## GRADO.DIModerado -0.452132 0.172500 -2.621 0.008766 **
## GRADO.DIGrave 0.289646 0.176579 1.640 0.100938
## GRADO.DIProfundo 0.322401 0.224081 1.439 0.150216
## TINETTI.TOTAL -0.122306 0.018413 -6.642 3.09e-11 ***
## DEF.ATENCISi 0.497892 0.145674 3.418 0.000631 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 689.10 on 110 degrees of freedom
## Residual deviance: 366.21 on 102 degrees of freedom
## AIC: 574.93
##
## Number of Fisher Scoring iterations: 6De la tabla de Coeficientes observamos que, en presencia de las otras variables, SUMA.ALT.MARCHA no aparece como significativa y podría eliminarse del modelo.
5.1 Interpretación de los coeficientes
El significado de los resultados de la regresión de contar es parecido al de la regresión logística (tanto para las variables numéricas como para las variables categóricas). Recordemos que el modelo es: \[ log(\lambda) = b_0+b_1X_1+...+b_jX_j+...b_pX_p \]
Vamos a ver como cambia la variable dependiente cuando una variable \(X_j\) cualquiera aumenta una unidad, es decir pasa de valer \(X_j\) a valer \(X_j+1\). Llamemos \(\lambda_0\) al valor \(\lambda\) en el estado inicial, y \(\lambda_1\) al valor \(\lambda\) en \(X_j+1\) (suponemos que el resto de variables permanece constante).
\[ log(\lambda_1)- log(\lambda_0) = b_j \ \ \rightarrow \ \ log \left( \dfrac{\lambda_1}{\lambda_0}\right) = b_j \ \ \rightarrow \ \ \left( \dfrac{\lambda_1}{\lambda_0}\right) = e^{b_j}. \]
Entonces, si elevamos \(e\) a cada coeficiente de la regresión \(b_j\), tenemos el cambio en el número medio de la variable \(Y\) cuando \(X_j\) crece una unidad.
Entonces, en este caso, el término edad tiene coeficiente \(0.03\), lo que indica que, por unidad aumentada en edad (año) aumenta el número medio de caídas \(e^{0.03}=1.03\) (digamos una caida por año, en media).
Con una unidad aumentada en TINETTI.TOTAL el número medio de caídas aumenta \(e^{-0.12}=0.88\) (es decir, desciende).
Con las variables binarias o categóricas se razona igual que en la regresión logística. Por ejemplo, al pasar del estado Deficit.atencion=“no” al estado “sí” el número medio de caidas aumenta \(e^{0.49}=1.63\).
5.2 Test de sobre-dispersion
Con el siguiente test se chequea si existe sobre-dispersión (algo que aquí está bastante claro, y no sería necesario)
dispersiontest(fm_pois)##
## Overdispersion test
##
## data: fm_pois
## z = 3.5077, p-value = 0.000226
## alternative hypothesis: true dispersion is greater than 1
## sample estimates:
## dispersion
## 3.827756Si no existiera esta sobredispersión, y el modelo de Poisson pudiera ser adecuado, los resultados de la estimación de los coeficientes pueden mejorarse mediante
coeftest(fm_pois , vcov = sandwich)##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.968543 1.621384 1.2141 0.224705
## EDAD 0.034567 0.011700 2.9545 0.003132 **
## NUM.MEDIC 0.037205 0.045806 0.8122 0.416659
## SUMA.ALT.MARCHA -0.070213 0.120324 -0.5835 0.559535
## GRADO.DIModerado -0.452132 0.442620 -1.0215 0.307022
## GRADO.DIGrave 0.289646 0.382423 0.7574 0.448812
## GRADO.DIProfundo 0.322401 0.427881 0.7535 0.451160
## TINETTI.TOTAL -0.122306 0.047919 -2.5523 0.010700 *
## DEF.ATENCISi 0.497892 0.335021 1.4862 0.137239
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 15.3 Modelo quasi-poisson
Cuando un modelo de Poisson no se ajusta bien a los datos, el siguiente que debe probarse es el llamado Quasi-Poisson, que es una generalización donde la varianza puede ser más grande que la media. Se estima como:
fm_qpois <- glm(Caidas~ ., data = rc, family = quasipoisson)
summary(fm_qpois)##
## Call:
## glm(formula = Caidas ~ ., family = quasipoisson, data = rc)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -3.756 -1.529 -1.048 0.524 5.392
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.96854 1.28657 1.530 0.12909
## EDAD 0.03457 0.01389 2.489 0.01442 *
## NUM.MEDIC 0.03721 0.03929 0.947 0.34595
## SUMA.ALT.MARCHA -0.07021 0.09968 -0.704 0.48279
## GRADO.DIModerado -0.45213 0.34954 -1.294 0.19876
## GRADO.DIGrave 0.28965 0.35780 0.810 0.42011
## GRADO.DIProfundo 0.32240 0.45406 0.710 0.47930
## TINETTI.TOTAL -0.12231 0.03731 -3.278 0.00143 **
## DEF.ATENCISi 0.49789 0.29518 1.687 0.09471 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for quasipoisson family taken to be 4.105945)
##
## Null deviance: 689.10 on 110 degrees of freedom
## Residual deviance: 366.21 on 102 degrees of freedom
## AIC: NA
##
## Number of Fisher Scoring iterations: 6Si nos fijamos en los resultados de la devianza y del AIC para saber cuál de los dos modelos es mejor, vemos que las devianzas son iguales. El AIC para este modelo no puede calcularse, con lo cual no se puede comparar. En principio, ambos modelos serían iguales. Además, puede realizarse una prueba de tipo ANOVA para comparar los modelos:
anova(fm_pois, fm_qpois, test="Chisq")## Analysis of Deviance Table
##
## Model 1: Caidas ~ EDAD + NUM.MEDIC + SUMA.ALT.MARCHA + GRADO.DI + TINETTI.TOTAL +
## DEF.ATENCI
## Model 2: Caidas ~ EDAD + NUM.MEDIC + SUMA.ALT.MARCHA + GRADO.DI + TINETTI.TOTAL +
## DEF.ATENCI
## Resid. Df Resid. Dev Df Deviance Pr(>Chi)
## 1 102 366.21
## 2 102 366.21 0 0que es un test de diferencia entre devianzas, para contrastar si las diferencias podrían ser significativas. En este caso, al ser iguales, ambos modelos son iguales.
5.4 Modelo de variable Binomial Negativa
El paso siguiente es intentar ver como se adecua el modelo, en vez de a una variable de Poisson, a una variable binomial negativa. La definición de variable binomial negativa aparece en cualquier manual de estadística básica. A nosotros nos interesa que es una variable que, como la de Poisson, en principio toma valores desde \(0\) en adelante, pero la varianza puede ser bastante más grande que la media (que es lo que nos ocurre). Para estimar el modelo, debemos cambiar el comando glm por glm.nb.
fm_nbin <- glm.nb(Caidas ~ ., data = rc)
summary(fm_nbin)##
## Call:
## glm.nb(formula = Caidas ~ ., data = rc, init.theta = 0.9078070053,
## link = log)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.8575 -1.1042 -0.5851 0.3121 2.9844
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.402208 1.496275 0.937 0.3487
## EDAD 0.032884 0.013797 2.383 0.0171 *
## NUM.MEDIC 0.113835 0.046668 2.439 0.0147 *
## SUMA.ALT.MARCHA -0.001528 0.125714 -0.012 0.9903
## GRADO.DIModerado -0.426750 0.322946 -1.321 0.1864
## GRADO.DIGrave 0.290116 0.366806 0.791 0.4290
## GRADO.DIProfundo -0.341921 0.634352 -0.539 0.5899
## TINETTI.TOTAL -0.116244 0.045498 -2.555 0.0106 *
## DEF.ATENCISi 0.562301 0.320853 1.753 0.0797 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for Negative Binomial(0.9078) family taken to be 1)
##
## Null deviance: 191.3 on 110 degrees of freedom
## Residual deviance: 111.0 on 102 degrees of freedom
## AIC: 426.77
##
## Number of Fisher Scoring iterations: 1
##
##
## Theta: 0.908
## Std. Err.: 0.212
##
## 2 x log-likelihood: -406.770De entrada, comparando las devianzas y el AIC de este modelo con los términos correspondientes en el modelo de Poisson, se reducen bastante los números. El AIC pasa de \(574.93\) a \(426.77\), con lo que este último modelo es mejor que el anterior.
Observamos que no aparecen como variables significativas ni SUM.ALT.MARCHA ni el grado de discapacidad. Vamos entonces a ver si mejora el modelo quitando estas variables (ahora debemos escribirlas, puesto que no son todas las del conjunto de datos que estamos usando):
fm_nbin2<- glm.nb(Caidas ~ EDAD+NUM.MEDIC+TINETTI.TOTAL+
DEF.ATENCI, data = rc)
summary(fm_nbin2) ##
## Call:
## glm.nb(formula = Caidas ~ EDAD + NUM.MEDIC + TINETTI.TOTAL +
## DEF.ATENCI, data = rc, init.theta = 0.8358529869, link = log)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.6749 -1.0997 -0.7848 0.3664 2.4576
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.84171 0.95617 1.926 0.0541 .
## EDAD 0.03017 0.01321 2.283 0.0224 *
## NUM.MEDIC 0.11074 0.04680 2.367 0.0180 *
## TINETTI.TOTAL -0.13122 0.02822 -4.651 3.31e-06 ***
## DEF.ATENCISi 0.49220 0.31834 1.546 0.1221
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for Negative Binomial(0.8359) family taken to be 1)
##
## Null deviance: 182.12 on 110 degrees of freedom
## Residual deviance: 110.99 on 106 degrees of freedom
## AIC: 423.39
##
## Number of Fisher Scoring iterations: 1
##
##
## Theta: 0.836
## Std. Err.: 0.190
##
## 2 x log-likelihood: -411.393Los resultados de AIC y devianzas son muy similares al modelo anterior. Ahora, en cambio, no aparece como significativo el deficit de atención.
Vamos a hacer un test ANOVA para comparar entre los dos modelos
anova(fm_nbin2, fm_nbin, test="Chisq")## Likelihood ratio tests of Negative Binomial Models
##
## Response: Caidas
## Model
## 1 EDAD + NUM.MEDIC + TINETTI.TOTAL + DEF.ATENCI
## 2 EDAD + NUM.MEDIC + SUMA.ALT.MARCHA + GRADO.DI + TINETTI.TOTAL + DEF.ATENCI
## theta Resid. df 2 x log-lik. Test df LR stat. Pr(Chi)
## 1 0.835853 106 -411.3932
## 2 0.907807 102 -406.7699 1 vs 2 4 4.623277 0.3281793El \(p-valor\) nos dice que los modelos no pueden considerarse diferentes.
Por último, vamos a probar un modelo sin la variable Grado de discapacidad.
fm_nbin3<- glm.nb(Caidas ~ EDAD+NUM.MEDIC+TINETTI.TOTAL, data = rc)
summary(fm_nbin3) ##
## Call:
## glm.nb(formula = Caidas ~ EDAD + NUM.MEDIC + TINETTI.TOTAL, data = rc,
## init.theta = 0.7960430615, link = log)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.7477 -1.1076 -0.7665 0.2744 2.3429
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 2.28570 0.95016 2.406 0.01615 *
## EDAD 0.02952 0.01335 2.211 0.02706 *
## NUM.MEDIC 0.12253 0.04592 2.668 0.00762 **
## TINETTI.TOTAL -0.14484 0.02803 -5.167 2.38e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for Negative Binomial(0.796) family taken to be 1)
##
## Null deviance: 176.84 on 110 degrees of freedom
## Residual deviance: 110.74 on 107 degrees of freedom
## AIC: 423.95
##
## Number of Fisher Scoring iterations: 1
##
##
## Theta: 0.796
## Std. Err.: 0.178
##
## 2 x log-likelihood: -413.950Volvemos a obtener unos resultados de AIC muy similares a los dos de antes.
5.5 Modelos con exceso de ceros
El siguiente paso que se debe dar es considerar los llamados modelos que tienen “excesos de ceros”. Se utilizan cuando, como en nuestro caso, la frecuencia de aparición de ceros es bastante más grande que la frecuencia de aparición del resto de valores. Sin embargo, hay ocasiones en que tienen mucho mayor interés. Por ejemplo, el número de partes a un seguro al año por sus asegurados (la mayoría son cero y las probabilidades de valores mayores van disminuyendo). En seguros de salud, las visitas al centro hospitalario. El número de muertos por una enfermedad al año entre quienes la padecieron….
Este tipo de situaciones se modelizan con los llamados “Zero inflated” models (modelos inflados de ceros) y los “Hurdle models” (modelos con salto). Un modelo con salto es un modelo con exceso de ceros donde además puede existir una diferencia sustancial entre el valor cero y el resto de valores. Por ejemplo, si se realiza una encuesta en una población sobre el número de cigarrillos consumidos en una semana, hay dos subpoblaciones: los que fuman y los que no. Los que no fuman siempre van a dar el valor cero, y los que fuman pueden dar el valor cero o mayor. Si consideramos un grupo de personas con discapacidad y otro sin ella, el valor cero de caidas sería posiblemente diferente para cada subpoblación y habría que estudiarlo.
Los modelos con exceso de ceros se pueden analizar de manera equivalente a la que hemos hecho. Al estudiarlo con estos datos, no se han observado mejoras con respecto al modelo que utiliza la variable binomial negativa.
Dado que, por un lado, estos modelos son matemáticamente más complejos, suele ser necesario un número mayor de datos del que tenemos en esta ocasión para poder realizar estimaciones fiables. Por otro lado, al no haber obtenido mejoras con el modelo de variable binomial negativa, el principio de parsimonia (si se tienen varios modelos equivalentes, nos quedamos con el más sencillo) nos dice que nos quedemos con este último.
No sería necesario, puesto que hemos visto ya que el modelo de binomial negativa es mejor que el de Poisson, pero vamos a realizar una comparación entre ellos (Poisson, Quasi-Poisson, fm_nbin y fm_nbin2):
## ML-Pois Quasi-Pois NB NB_2
## (Intercept) 1.96854344 1.96854344 1.402207747 1.84171006
## EDAD 0.03456708 0.03456708 0.032883934 0.03017313
## NUM.MEDIC 0.03720516 0.03720516 0.113834902 0.11074251
## SUMA.ALT.MARCHA -0.07021307 -0.07021307 -0.001527688 -0.13122020
## GRADO.DIModerado -0.45213222 -0.45213222 -0.426749869 0.49220387
## GRADO.DIGrave 0.28964594 0.28964594 0.290116498 NA
## GRADO.DIProfundo 0.32240053 0.32240053 -0.341921453 NA
## TINETTI.TOTAL -0.12230619 -0.12230619 -0.116244234 NAConteo de ceros:
## zero-countsb
round(c("Obs" = sum(Bdatos$Caidas < 1),
"ML-Pois" = sum(dpois(0, fitted(fm_pois))),
"NB" = sum(dnbinom(0, mu = fitted(fm_nbin), size = fm_nbin$theta)),
"NB_2" = sum(dnbinom(0, mu = fitted(fm_nbin2), size = fm_nbin$theta))
))## Obs ML-Pois NB NB_2
## 47 26 46 45Este último resultado nos dice que, en los datos, se observaron \(47\) ceros. El modelo de Poisson solo predeciría \(26\), mientras que los dos de binomial negativa predecirían \(46\) y \(45\).
## esto simplemente para comparar los AIC de los modelos
## y elegir el menor
tmp <- c( AIC(fm_pois), AIC(fm_qpois), AIC(fm_nbin), AIC(fm_nbin2) )
( BestAIC<- which.min(tmp) )## [1] 4Predicción datos según modelos (comparamos Poisson con los dos binomial negativa)
po.p <- predprob(fm_pois) %>% colMeans
# po.qp <- predprob(fm_qpois) %>% colMeans
po.nb <- predprob(fm_nbin) %>% colMeans
po.nb.2 <- predprob(fm_nbin2) %>% colMeans
df <- data.frame(x = 0:max(Bdatos$Caidas), Poisson = po.p,
# QPoisson=po.qp,
NegBin = po.nb,
NegBin2 = po.nb.2)
obs <- table(Bdatos$Caidas) %>% prop.table() %>% data.frame #Observed
names(obs) <- c("x", 'Observed')
N_total=sum(table(Bdatos$Caidas))
linear <- glm(Caidas ~ ., data = rc) #identity link, OLS
p1 <- predict(linear) %>% round() %>% table %>% prop.table %>% data.frame #for OLS
names(p1) <- c('x', 'OLS')
tmp <- merge(p1, obs, by = 'x', all = T)
tmp$x <- as.numeric(as.character(tmp$x))
comb <- merge(tmp, df, by = 'x', all = T)
comb[is.na(comb)] <- 0
comb2 <- comb[1:45, ] #just for the first 11 results, including zero
mm <- melt(comb2, id.vars = 'x', value.name = 'prob', var.name = 'Model')
names(mm)<-c("x","Model","prob")
mm <- filter(mm, Model != "OLS") #can include the linear model too if you want
mm<-dplyr::filter(mm, mm$x >= 0)
mm <- mutate(mm, Cites_predicted=prob*N_total)
mm.obs<-dplyr::filter(mm, mm$Model=="Observed")
mm.modelo<-dplyr::filter(mm, mm$Model%in% c("Poisson","NegBin","NegBin2"))
ggplot()+
geom_bar(stat="identity", colour = "lemonchiffon1", fill = "lemonchiffon1" , width=0.4, aes(x=x,y=Cites_predicted),data=mm.obs)+
geom_line(aes(x=x,y=Cites_predicted,color=Model),size=0.8, data=mm.modelo)+
geom_point(aes(x=x,y=Cites_predicted, color=Model),data=mm.modelo )+
labs(x = "Number of falls", y = 'Individuals ',
title = "Models for number of falls per individual. Bars are the real number. ") 
Si se quiere hacer esto mismo, pero con una tabla en vez de gráfico, quitar los comentarios al siguiente código y ejecutarlo:
# mm_tabla<- comb2
# mm_tabla<-dplyr::filter(mm_tabla, x >= 0)
#
# N1=nrow(mm_tabla)-1
# NN=N_total
#
# ss<-data.frame(Caidas)
# names(ss)<-"R"
# s2<-(ss %>%
# mutate(R = factor(R, levels=0:N1)) %>%
# group_by(R) %>%
# summarise(freq=n()) %>%
# complete(R, fill=list(freq=0)) )
#
# tabla_comp<- cbind(s2,mm_tabla)
# Tabla_Final <- (tabla_comp %>%
# mutate( Caidas=R,
# Real_data = freq,
# Observed=Observed*NN,
# Poisson=round(Poisson*NN),
# NegBin=round(NegBin*NN),
# NegBin2=round(NegBin2*NN) ) %>%
# select(Caidas,Real_data,Poisson,NegBin,NegBin2))
# (Tabla_Final)
#
# (Error_Poisson<- sum( abs(Tabla_Final$Real_data-Tabla_Final$Poisson)))
# (Error_Negbin <-sum( abs(Tabla_Final$Real_data-Tabla_Final$NegBin)))
# (Error_Negbin2<-sum( abs(Tabla_Final$Real_data-Tabla_Final$NegBin2)))5.6 Interpretacion de coeficientes
Si nos quedamos con el modelo fm_nbin podemos centrarnos ya en los coeficientes y sus intervalos de confianza. Primero obtenemos los coeficientes y sus intervalos de confianza. Luego, exponenciamos ambos porque es lo que nos interesa.
cof=fm_nbin$coefficients
cbind(Coeficientes=cof , confint(fm_nbin)) # coeficientes y sus IC## Waiting for profiling to be done...## Coeficientes 2.5 % 97.5 %
## (Intercept) 1.402207747 -1.573895454 4.44075434
## EDAD 0.032883934 0.004966131 0.06180279
## NUM.MEDIC 0.113834902 0.019643639 0.21118514
## SUMA.ALT.MARCHA -0.001527688 -0.221300257 0.22875786
## GRADO.DIModerado -0.426749869 -1.062083438 0.20630417
## GRADO.DIGrave 0.290116498 -0.459494263 1.06021775
## GRADO.DIProfundo -0.341921453 -1.588176953 0.99986885
## TINETTI.TOTAL -0.116244234 -0.208407290 -0.02748428
## DEF.ATENCISi 0.562300592 -0.059772139 1.21396780exp(cbind(Coeficientes =cof, confint(fm_nbin))) # exponencial de los coef. e IC de estos## Waiting for profiling to be done...## Coeficientes 2.5 % 97.5 %
## (Intercept) 4.0641627 0.2072363 84.838915
## EDAD 1.0334306 1.0049785 1.063753
## NUM.MEDIC 1.1205671 1.0198378 1.235141
## SUMA.ALT.MARCHA 0.9984735 0.8014760 1.257038
## GRADO.DIModerado 0.6526268 0.3457347 1.229127
## GRADO.DIGrave 1.3365832 0.6316030 2.887000
## GRADO.DIProfundo 0.7104040 0.2042977 2.717925
## TINETTI.TOTAL 0.8902578 0.8118763 0.972890
## DEF.ATENCISi 1.7547047 0.9419791 3.366817La ecuación de regresión podría expresarse:
\(log(\) número medio de caidas) = 1.4022077 + (0.0328839) \(\cdot EDAD\) + (0.1138349) \(\cdot NUM.MEDIC\) + (-0.0015277) \(\cdot SUM.ALT.MARCHA\) + (-0.4267499) \(\cdot GRADO.DImoderado\) + (0.2901165) \(\cdot GRADO.DIGrave\) + (-0.3419215) \(\cdot GRADO.DIProfundo\) + (-0.1162442) \(\cdot TINETTI.TOTAL\) + (0.5623006) \(\cdot DEF.ATENCIsi\)
## Gráfico de efectos
plot(allEffects(fm_nbin2))
5.7 Colinealidad
Se dice que existe colinealidad en el modelo (o entre las variables explicativas) si existe dependencia entre las variables explicativas. Cualquiera de los modelos (regresión lineal, logística o de contar) puede presentar problemas de alta colinealidad, lo que se traduce en malas estimaciones de los coeficientes o en modelos inadecuados.
Una forma de medir la colinealidad es el VIF (Variance Inflation Factor). Un VIF menor a \(5\) indica baja colinealidad de una variable explicativa con el resto (en la mayoría de los casos, se considera un valor inferior a \(2\) como prácticamente nulo. Un valor entre \(5\) y \(10\) indica una colinealidad moderada, y valores mayores que \(10\) indicarían un exceso que haría inadecuada la presencia de la variable con dicho valor.
Se puede calcular (y presentar gráficamente) mediante
(x<-check_collinearity(fm_nbin))## # Check for Multicollinearity
##
## Low Correlation
##
## Parameter VIF Increased SE
## EDAD 1.29 1.14
## NUM.MEDIC 1.44 1.20
## SUMA.ALT.MARCHA 3.11 1.76
## GRADO.DI 1.56 1.25
## TINETTI.TOTAL 3.24 1.80
## DEF.ATENCI 1.26 1.12plot(x)
Vemos que todas las variables toman valores bajos. Los más altos son los VIF de SUM.ALT.MARCHA y TINETTI.TOTAL. Este valor quizá pudiera servir para considerar no incluir SUM.ALT.MARCHA, puesto que tampoco salía significativa.