#----Install Packages-----
#install.packages("car")
#install.packages("mlogit")
#------And then load these packages, along with the boot package.-----

library(car)
## Loading required package: carData
library(mlogit)
## Loading required package: dfidx
## 
## Attaching package: 'dfidx'
## The following object is masked from 'package:stats':
## 
##     filter
library(Rcmdr)
## Loading required package: splines
## Loading required package: RcmdrMisc
## Loading required package: sandwich
## Loading required package: effects
## Registered S3 methods overwritten by 'lme4':
##   method                          from
##   cooks.distance.influence.merMod car 
##   influence.merMod                car 
##   dfbeta.influence.merMod         car 
##   dfbetas.influence.merMod        car
## lattice theme set by effectsTheme()
## See ?effectsTheme for details.
## La interfaz R-Commander sólo funciona en sesiones interactivas
## 
## Attaching package: 'Rcmdr'
## The following object is masked from 'package:base':
## 
##     errorCondition

Los datos están en Eel.dat.

Resultado(variable dependiente): Cured(curado o no curado);

Predictor(variable independiente): Intervention(intervención o ningún tratamiento);

Predictor(variable independiente): Duration(el número de días antes del tratamiento que el paciente tenía el problema).

Los datos se pueden descargar de

eel.dat

#********************* Eel Example ********************

#Cargando los datos
eelData<-read.delim("E:/eel.dat", header = TRUE)

#look at first 6 cases of data
head(eelData)

Tenga en cuenta que tenemos tres variables en diferentes columnas. Los datos categóricos han sido ingresado como texto. Por ejemplo, la variable Curado se compone de las frases Cured y Not Curado. Cuando leemos datos que son cadenas de texto como esta, R amablemente las convierte a factores. No nos dice que ha hecho esto, simplemente lo hace. Cuando hacemos regresión logística, queremos hacerlo con números, no con palabras. En creación los factores R también asignaron provechosamente algunos números a las variables. No hay fin a como útil R intentará serlo. El problema es que los números que R ha asignado pueden no ser los números que queremos. De hecho, R crea niveles del factor tomando las cadenas de texto en orden alfabético y asignándoles valores numéricos ascendentes. En otras palabras, para Curado tenemos dos categorías y R habrá ordenado estas categorías alfabéticamente (es decir, “curado” y “no curado”). Entonces, Cured será la categoría de referencia porque es la primera.

Asimismo, para Intervención las categorías fueron Intervención y No Tratamiento, por lo que dada la La intervención en orden alfabético será la categoría de referencia. Sin embargo, tiene más sentido codificar ambas variables al revés. Para Cured sería bueno si Not Cured fuera la línea de base, o la primera categoría, porque entonces sabría que los coeficientes del modelo reflejan la probabilidad de curarse (que es lo que queremos saber) en lugar de la probabilidad de no curarnos. Del mismo modo, para Intervención Sería útil si Sin tratamiento fuera la primera categoría (es decir, la línea de base). Afortunadamente, el la función relevel () nos permite especificar la categoría de línea base para un factor. Toma la forma general:

#Alternatively Re-orders the levels of the factyor so that Not Cured and No Treatment are the baseline categories
eelData$Cured<-factor(eelData$Cured, levels = c("Not Cured", "Cured"))
eelData$Intervention<-factor(eelData$Intervention, levels = c("No Treatment", "Intervention"))
# specify the baseline category
eelData$Cured<-relevel(eelData$Cured, "Not Cured")
eelData$Intervention<-relevel(eelData$Intervention, "No Treatment")
#Create the two hierarchical models:

eelModel.1 <- glm(Cured ~ Intervention, data = eelData, family = binomial())
eelModel.2 <- glm(Cured ~ Intervention + Duration, data = eelData, family = binomial())

summary(eelModel.1)
## 
## Call:
## glm(formula = Cured ~ Intervention, family = binomial(), data = eelData)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -1.5940  -1.0579   0.8118   0.8118   1.3018  
## 
## Coefficients:
##                          Estimate Std. Error z value Pr(>|z|)   
## (Intercept)               -0.2877     0.2700  -1.065  0.28671   
## InterventionIntervention   1.2287     0.3998   3.074  0.00212 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 154.08  on 112  degrees of freedom
## Residual deviance: 144.16  on 111  degrees of freedom
## AIC: 148.16
## 
## Number of Fisher Scoring iterations: 4
summary(eelModel.2)
## 
## Call:
## glm(formula = Cured ~ Intervention + Duration, family = binomial(), 
##     data = eelData)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -1.6025  -1.0572   0.8107   0.8161   1.3095  
## 
## Coefficients:
##                           Estimate Std. Error z value Pr(>|z|)   
## (Intercept)              -0.234660   1.220563  -0.192  0.84754   
## InterventionIntervention  1.233532   0.414565   2.975  0.00293 **
## Duration                 -0.007835   0.175913  -0.045  0.96447   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 154.08  on 112  degrees of freedom
## Residual deviance: 144.16  on 110  degrees of freedom
## AIC: 150.16
## 
## Number of Fisher Scoring iterations: 4

Recuerde que los valores más grandes de la estadística de desviación indican un ajuste deficiente modelos estadísticos. R proporciona dos estadísticas de desviación: la desviación nula y la desviación residual.

desviación. La desviación nula es la desviación del modelo que no contiene otros predictores que la constante, en otras palabras, -2LL (línea de base). La desviación residual es la desviación para el modelo, en otras palabras, –2LL (nuevo).

En esta etapa del análisis, el valor de la desviación del modelo debe ser menor que el valor para el modelo nulo (cuando solo se incluyó la constante en el modelo) porque valores más bajos de -2LL indican que el modelo predice la variable de resultado más precisamente. Para el modelo nulo, −2LL = 154.08, pero cuando se ha incluido Intervención este valor se ha reducido a 144,16. Esta reducción nos dice que el modelo es mejor en predecir si alguien estaba curado de lo que estaba antes de que se agregara la Intervención.

La cuestión de cuánto mejor predice el modelo la variable de resultado puede ser evaluado utilizando el modelo estadístico de chi-cuadrado, que mide la diferencia entre los modelo tal como está actualmente y el modelo cuando solo se incluyó la constante. El valor del estadístico chi-cuadrado del modelo funciona según este principio y, por lo tanto, es igual a −2LL con Intervención incluida menos el valor de −2LL cuando solo el La constante estaba en el modelo (154,08 - 144,16 = 9,92). Este valor tiene una distribución de chi-cuadrado por lo que su significación estadística se puede calcular fácilmente. En este ejemplo, el valor es significativo a un nivel de .05, por lo que podemos decir que, en general, el modelo predice si un paciente se cura o no mejora significativamente de lo que estaba con solo la constante incluida. El modelo chi-cuadrado es un análogo de la prueba F para la regresión lineal. En un mundo ideal, nos gustaría ver un −2LL general no significativo (lo que indica que el cantidad de datos inexplicables es mínima) y un modelo estadístico de chi-cuadrado altamente significativo (lo que indica que el modelo que incluye los predictores es significativamente mejor que sin los predictores). Sin embargo, en realidad es posible que ambas estadísticas sean muy significativas. Podemos usar R para calcular automáticamente el chi-cuadrado del modelo y su significado. Nosotros podemos hacer esto tratando el modelo de salida como datos. El objeto eelModel.1 tiene una serie de

#Just to prove what the null deviance is
eelModel.0 <- glm(Cured ~ 1, data = eelData, family = binomial())
summary(eelModel.0)
## 
## Call:
## glm(formula = Cured ~ 1, family = binomial(), data = eelData)
## 
## Deviance Residuals: 
##    Min      1Q  Median      3Q     Max  
## -1.309  -1.309   1.052   1.052   1.052  
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)
## (Intercept)   0.3032     0.1903   1.593    0.111
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 154.08  on 112  degrees of freedom
## Residual deviance: 154.08  on 112  degrees of freedom
## AIC: 156.08
## 
## Number of Fisher Scoring iterations: 4

variables asociadas a él. Dos de estos son la desviación y la desviación nula. Restando la desviación de la desviación nula podemos encontrar la mejora, lo que nos da un chi cuadrado. Podemos hacer referencia a estas variables de la misma manera que a cualquier otra: agregando el nombre de variable al modelo usando el signo de dólar. Para calcular este valor ejecute:

modelChi <- eelModel.1$null.deviance - eelModel.1$deviance
chidf <- eelModel.1$df.null - eelModel.1$df.residual
chisq.prob <- 1 - pchisq(modelChi, chidf)
modelChi; chidf; chisq.prob
## [1] 9.926201
## [1] 1
## [1] 0.001629425

En otras palabras, el valor p es .002 (redondeado a tres lugares decimales); porque esta probabilidad es menor que .05, podemos rechazar la hipótesis nula de que el modelo no es mejor que oportunidad de predecir el resultado. Este valor es el valor p de la razón de verosimilitud del modelo porque solo teníamos un predictor en el modelo. Podemos informar que incluida la Intervención produjo una mejora significativa en el ajuste del modelo, \(\chi^2(1) = 9.93,\quad p = .002\).

modelChi2 <- eelModel.2$null.deviance - eelModel.2$deviance
chidf2 <- eelModel.2$df.null - eelModel.2$df.residual
chisq.prob2 <- 1 - pchisq(modelChi2, chidf2)
modelChi2; chidf2; chisq.prob2
## [1] 9.928185
## [1] 2
## [1] 0.006984287

El estadístico crucial es el estadístico \(z\), que tiene una distribución normal y nos dice si el coeficiente \(b\) para ese predictor es significativamente diferente de cero. Si el coeficiente es significativamente diferente de cero, entonces podemos suponer que el predictor está haciendo una contribución significativa a la predicción del resultado (Y). La estadística \(z\) debe usarse con precaución porque cuando el coeficiente de regresión \((b)\) es grande, el error estándar tiende a inflarse, lo que resulta en \(z\) - subestimación de la estadística(véase Menard, 1995). Sin embargo, para estos datos, parece indicar que tener la intervención (o no) es un predictor significativo de si el paciente está curado(tenga en cuenta que la importancia de la estadística \(z\) es menor que \(0.05\)). Podríamos decir que la Intervención fue un predictor significativo de curación, \(b = 1.23,\quad z = 3.07, \quad p <.002\).

Podemos calcular un análogo de \(R\) usando la ecuación (8.7). Para estos datos, la estadística \(z\) y su \(df\) se pueden leer en la salida R (\(3.074\) y 1, respectivamente), y la desviación del modelo nulo fue \(154.08\). Por lo tanto, \(R\) se puede calcular como: \[ R = \sqrt{\frac{3.074^{2} -2\times 1}{154.08}} = 0.22 \] Sabemos que la medida de Hosmer y Lemeshow \(\left(R_{\mathrm{L}}^{2} \right)\) se calcula dividiendo el chi-cuadrado del modelo por el \(-2LL original.\) In En este ejemplo, el chi-cuadrado del modelo después de que se ingresó la Intervención en el modelo es \(9.93\) (calculado como modelChi arriba), y el \(-2 LL\) original (antes de ingresar cualquier variable) era \(154.08\) (la desviación.null). Entonces, \(R_{\mathrm{L}}^{2} = 9.93 /154.08 = .06\), que es diferente del valor que obtendríamos al elevar al cuadrado el valor de \(R\) dado arriba de \(\left(R^{ 2} =. 22^{2} = 0.05 \right)\). Podemos hacer que \(\mathrm{R}\) haga este cálculo por nosotros ejecutando:

R2.hl<-modelChi/eelModel.1$null.deviance
R2.hl
## [1] 0.06442071
R.cs <- 1 - exp ((eelModel.1$deviance - eelModel.1$null.deviance)/113)
R.cs
## [1] 0.08409487
R.n <- R.cs /( 1- ( exp (-(eelModel.1$null.deviance/ 113))))
R.n
## [1] 0.112992

El primer comando simplemente toma el valor del chi-cuadrado del modelo (que ya hemos calculado como modelChi, y lo divide por \(-2 L L\) para el modelo original (eelModel. 1\(\$\)null. Deviance)). Este es un análogo directo de la ecuación dada anteriormente en este capítulo. El segundo comando muestra el valor, que es \(.064\).

También vimos otras dos medidas de \(R^{2}\) , Cox y Snell y Nagelkerke. Hay funciones disponibles en R para calcularlas, pero son un poco difíciles de encontrar y usar. Sin embargo, es bastante fácil escribir comandos en R para calcularlos. Podemos escribir la ecuación para el estadístico de Cox y Snell como:

R.cs <- 1 - exp ((eelModel.1$deviance - eelModel.1$null.deviance) /113)
R.cs
## [1] 0.08409487

El primer comando usa −2LL para el modelo (eelModel.1\(\$\)deviance) y el modelo nulo (eelModel.1$null.deviance) y divide la diferencia por el tamaño de la muestra (en este caso 113, pero necesitará cambiar este valor para otros conjuntos de datos). El segundo comando se mostrará el resultado: un valor de .084. Podemos usar este valor para calcular la estimación de Nagelkerke además. Nuevamente, simplemente escribimos la ecuación en lenguaje R:

R.n<-R.cs /(1-(exp(-(eelModel.1$null.deviance/113))))
R.n
## [1] 0.112992
###############################################
# This section creates a function called      #
# logisticPseudoR2s().  To use it             #
# type logisticPseudoR2s(myLogisticModel)     #
###############################################
 logisticPseudoR2s <- function(LogModel) {
    dev <- LogModel$deviance 
    nullDev <- LogModel$null.deviance 
    modelN <-  length(LogModel$fitted.values)
    R.l <-  1 -  dev / nullDev
    R.cs <- 1- exp ( -(nullDev - dev) / modelN)
    R.n <- R.cs / ( 1 - ( exp (-(nullDev / modelN))))
    cat("Pseudo R^2 for logistic regression\n")
    cat("Hosmer and Lemeshow R^2  ", round(R.l, 3), "\n")
    cat("Cox and Snell R^2        ", round(R.cs, 3), "\n")
    cat("Nagelkerke R^2           ", round(R.n, 3),    "\n")
      }
############End of function ######################

logisticPseudoR2s(eelModel.1)
## Pseudo R^2 for logistic regression
## Hosmer and Lemeshow R^2   0.064 
## Cox and Snell R^2         0.084 
## Nagelkerke R^2            0.113

Lo último en lo que debemos pensar son las razones de probabilidades, que se describieron en . Para calcular el cambio en las probabilidades que resulta de un cambio de unidad en el predictor de este ejemplo, primero debemos calcular las probabilidades de que un paciente se cure dado que no recibió la intervención. Luego calculamos las probabilidades de que un paciente se cure dado que sí se sometió a la intervención. Finalmente, calculamos el cambio proporcional en estas dos probabilidades. Para calcular el primer conjunto de probabilidades, necesitamos calcular la probabilidad de que un paciente se cure dado que no se sometió a la intervención. La codificación del parámetro al comienzo de la salida nos dijo que los pacientes que no se sometieron a la intervención se codificaron con un 0, por lo que podemos usar este valor en lugar de \(X\). El valor de \(b_{1}\) se ha estimado para nosotros como \(1.229\) (ver Coeficientes: en la Salida), y el coeficiente para la constante puede tomarse de la misma tabla y es \(-0.288.\) Podemos calcular el probabilidades como:

\[\begin{array}{l} P(\text { Cured })=\frac{1}{1+e^{-\left(b_{0}+b_{1} X_{1}\right)}}=\frac{1}{1+e^{-[-0.288+(1.229 \times 0)]}}=.428 \\ P(\text { Not Cured })=1-P(\text { Cured })=1-.428=.527 \\ \text { odds }=\frac{.428}{.572}=0.748 \end{array}\]

Ahora calculamos lo mismo después de que la variable predictora haya cambiado en una unidad. En En este caso, debido a que la variable predictora es dicotómica, necesitamos calcular las probabilidades de un paciente que está siendo curado, dado que ha tenido la intervención. Entonces, el valor de la la variable de intervención, X, es ahora 1 (en lugar de 0). Los cálculos resultantes son los siguientes:

\[\begin{array}{l} P(\text { Cured })=\frac{1}{1+e^{-\left(b_{0}+b_{1} X_{1}\right)}}=\frac{1}{1+e^{-[-0.288+(1.229 \times 1)]}}=.719 \\ P(\text { Not Cured })=1-P(\text { Cured })=1-.719=.281 \\ \text { odds }=\frac{.719}{.281}=2.559 \end{array}\]

Ahora conocemos las probabilidades antes y después de un cambio de unidad en la variable predictora. Esto es ahora una simple cuestión para calcular el cambio proporcional en las probabilidades dividiendo las probabilidades después de un cambio de unidad en el predictor por las probabilidades antes de ese cambio:

\[\begin{aligned} \Delta \text{Odds} & = \frac{\text {Odds después de un cambio de unidad en el predictor}}{\text{Odds originales}} \\ & = \frac {2.56} {0.75} \\ & = 3,41 \end{aligned}\] 
#Compute odds ratio
exp(eelModel.1$coefficients)
##              (Intercept) InterventionIntervention 
##                 0.750000                 3.416667
exp(confint(eelModel.1))
## Waiting for profiling to be done...
##                              2.5 %   97.5 %
## (Intercept)              0.4374531 1.268674
## InterventionIntervention 1.5820127 7.625545

La razón de probabilidades para la Intervención (3.417) es la misma que calculamos anteriormente (teniendo en cuenta las diferencias en el redondeo). Podemos interpretar la razón de probabilidades en términos del cambio en las probabilidades. Si el valor es mayor que 1, indica que a medida que aumenta el predictor, aumentan las probabilidades de que se produzca el resultado. Por el contrario, un valor menor que 1 indica que a medida que aumenta el predictor, las probabilidades de que se produzca el resultado disminuyen. En este ejemplo, podemos decir que las probabilidades de que un paciente que recibe tratamiento se cure son \(3.42\) veces más altas que las de un paciente que no recibe tratamiento.

También podemos calcular intervalos de confianza para las razones de probabilidades. Para obtener los intervalos de confianza de los parámetros, usamos la función confint (), tal como lo hicimos para la regresión ordinaria. También podemos exponenciarlos con la función \(exp ()\). Para que se ejecuten los intervalos de confianza:

exp(confint(eelModel.1))
## Waiting for profiling to be done...
##                              2.5 %   97.5 %
## (Intercept)              0.4374531 1.268674
## InterventionIntervention 1.5820127 7.625545

Lo importante de este intervalo de confianza es que no cruza 1 (los valores en cada extremo del intervalo son mayores que 1). Esto es importante porque los valores mayores que 1 quieren decir que a medida que aumenta la variable predictora, también aumentan las probabilidades de (en este caso) curarse.

Los valores inferiores a 1 significan lo contrario: a medida que aumenta la variable predictora, las probabilidades de ser curado disminuyen . El hecho de que tanto el límite inferior como el superior de nuestro intervalo de confianza son arriba 1 nos da la confianza de que la dirección de la relación que hemos observado es cierto en la población (es decir, es probable que tener una intervención en comparación con no aumente las probabilidades de curarse). Si el límite inferior hubiera estado por debajo de 1, entonces nos diría que no es una posibilidad de que en la población la dirección de la relación sea la opuesta a lo que hemos observado. Esto significaría que no podríamos confiar en que nuestra intervención aumente las probabilidades de curarse.

#compare model1 and model 2
modelChi <- eelModel.1$deviance - eelModel.2$deviance
chidf <- eelModel.1$df.residual - eelModel.2$df.residual
chisq.prob <- 1 - pchisq(modelChi, chidf)
modelChi; chidf; chisq.prob
## [1] 0.001983528
## [1] 1
## [1] 0.9644765
anova(eelModel.1, eelModel.2)
#Diagnostics for model 1

eelData$predicted.probabilities<-fitted(eelModel.1)
eelData$standardized.residuals<-rstandard(eelModel.1)
eelData$studentized.residuals<-rstudent(eelModel.1)
eelData$dfbeta<-dfbeta(eelModel.1)
eelData$dffit<-dffits(eelModel.1)
eelData $leverage<-hatvalues(eelModel.1)
head(eelData[, c("Cured", "Intervention", "Duration", "predicted.probabilities")])
eelData[, c("leverage", "studentized.residuals", "dfbeta")]

El ajuste general del modelo final se muestra mediante el estadístico de desviación y su estadístico chi-cuadrado asociado. Si el significado de la estadística de chi-cuadrado es menor que .05, entonces el modelo es un ajuste significativo a los datos.