I. Regresión lineal simple y polinomial:

1. Datos

library(readxl)
tg <-mod_empiricos <- read_excel("mod_empiricos.xlsx", 
    sheet = "trigliceridos")
kable(tg)

cintura	trigliceridos
32.95	28
33.00	40
33.50	50
33.50	61
34.50	66
37.00	70
35.00	79
36.85	85
36.45	87
40.65	89
40.50	90
43.50	107
42.00	109
39.40	111
40.20	112
44.90	119
43.50	121
44.80	130
41.50	149
45.00	159

2. Selección de la varible dependiente y varible independiente

En un estudio de los factores asociados a las enfermedades cardiovasculares Christian, A. et al. encontraron que los pacientes con una circunferencia de cintura mayor a la de los niveles recomendados son más propensos a tener hipertensión, bajos niveles de colesterol HDL, altos niveles de triglicéridos, entre otros en comparación con aquellos con una circunferencia de la cintura normal. Estos resultados concuerdan con la hipótesis alterna: el número de triglicéridos en la sangre aumenta con la circunferencia de la cintura. Por lo tanto, los triglicéridos en la sangre (mg/dl) son la variable dependiente.

3. Gráfica de puntos

plot(tg$cintura, tg$trigliceridos, type = "p",
     xlim = c(30,46),
     ylim = c(27, 160),
     xlab = "Circunferencia de la cintura", 
     ylab = "Triglicéridos en la sangre", 
     asp = NA)

Figura 1. Esta gráfica muestra la relación entre los niveles de triglicéridos en la sangre (mg/dl) y la circunferencia de la cintura (en pulgadas).

4. Cálculo de estadísticos

Estimadores de los coeficientes de regresión para un modelo lineal simple y para un modelo polinomial

Modelo lineal simple

rls <- lm(tg$trigliceridos ~ tg$cintura)
summary(rls)

## 
## Call:
## lm(formula = tg$trigliceridos ~ tg$cintura)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -21.814 -11.236  -4.454  10.281  37.349 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -188.4947    33.3154  -5.658 2.29e-05 ***
## tg$cintura     7.2324     0.8508   8.501 1.02e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 15.94 on 18 degrees of freedom
## Multiple R-squared:  0.8006, Adjusted R-squared:  0.7895 
## F-statistic: 72.27 on 1 and 18 DF,  p-value: 1.021e-07

Modelo polinomial

rpol <- lm(trigliceridos ~ cintura + I(cintura^2), data=tg)
summary(rpol)

## 
## Call:
## lm(formula = trigliceridos ~ cintura + I(cintura^2), data = tg)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -20.592 -12.224  -0.563   8.832  34.320 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)
## (Intercept)  -612.0084   395.7948  -1.546    0.140
## cintura        29.2756    20.5454   1.425    0.172
## I(cintura^2)   -0.2835     0.2640  -1.074    0.298
## 
## Residual standard error: 15.88 on 17 degrees of freedom
## Multiple R-squared:  0.8133, Adjusted R-squared:  0.7913 
## F-statistic: 37.02 on 2 and 17 DF,  p-value: 6.389e-07

library(tab)
library(sjPlot)

## Learn more about sjPlot with 'browseVignettes("sjPlot")'.

p.tg<- lm(trigliceridos ~ cintura + I(cintura^2), data=tg)
m.tg <- lm(trigliceridos ~ cintura, data = tg)
tab_model(m.tg, p.tg)

	trigliceridos			trigliceridos
Predictors	Estimates	CI	p	Estimates	CI	p
(Intercept)	-188.49	-258.49 – -118.50	<0.001	-612.01	-1447.06 – 223.05	0.140
cintura	7.23	5.45 – 9.02	<0.001	29.28	-14.07 – 72.62	0.172
I(cintura^2)				-0.28	-0.84 – 0.27	0.298
Observations	20			20
R² / R² adjusted	0.801 / 0.790			0.813 / 0.791

Tabla 1. Tabla comparativa de los valores de los parámetros para ambos modelos (modelo lineal simple y modelo polinomial).

5. Gráfica de residuales

Modelo lineal simple

plot(tg$cintura,residuals(rls),
     xlab="Cintura (pulgadas)",
     ylab="Residuales")

Modelo polinomial

plot(tg$cintura,residuals(rpol),
     xlab="Cintura (pulgadas)",
     ylab="Residuales")

Las gráficas de residuales demuestran que los errores de las predicciones se distribuyen de forma aleatoria y no se observan tendencias en los datos.

6. Gráfica de líneas de regresión con intervalo de confianza

Modelo lineal simple (Figura 2)

library(ggplot2)
ggplot(data=tg, aes(x=cintura, y=trigliceridos)) +
  geom_point(pch=19, color="green", size=2) +
  geom_smooth(method="lm", color="purple", linetype=2) +
  labs(x="Cintura (pulgadas)", y="Trigliceridos en la sangre (mg/dl)")

Modelo polinomial (Figura 3)

library(ggplot2)
ggplot(tg, aes(x = cintura, y = trigliceridos)) +
  geom_point(aes(x = cintura, y = trigliceridos)) +
  stat_smooth(aes(), method = "lm", formula = y ~ poly(x,2), se =TRUE, size = 2,color="cyan") +
  labs(x="Cintura (pulgadas)", y="Trigliceridos en la sangre (mg/dl)")

7. Explicación y discusión de resultados

En la primera gráfica (Figura 1) se observó que las variables (Circunferencia de la cintura y triglicéridos en la sangre) están correlacionadas ya que se observa una tendencia en ella. Esta correlación no es completamente lineal ya que se observó dispersión en los datos de la gráfica. Al realizar los cálculos estadísticos se encontró que ambos modelos mostraron valores similares de R cuadrado y valores de p menores a cero. Al realizar las gráficas de los residuales se observó que los errores de las predicciones se distribuyen de forma aleatoria y no hay tendencias en los datos.

En la Figura 2 se graficaron los datos junto al intervalo de confianza y se logró observar una mejor relación entre los datos (niveles de triglicéridos y circunferencia de la cintura).Se puede establecer que gran parte de la variación en la circunferencia de la cintura se puede explicar por variaciones en la masa de grasa y las áreas de grasa abdominal (Maki, 1997), las cuales estan directamente relacionadas a los niveles de triglicéridos. Al utilizarse el modelo polinomial (Figura 3) se pudo ver una ligera mejora en la correlación entre ambas variables. También, al utilizarse el modelo polinomial se pudo observar mayor cantidad de valores cercanos al intervalo de confianza, lo cual permite determinar que es un mejor modelo. En conclusión, el modelo polinomial (Figura 3) es mejor ya que explica de mejor forma la relación entre la circunferencia de la cintura y los triglicéridos en la sangre, los cuales aumentan al ser mayor la circunferencia.

Referencias (regresión lineal)

Christian, A. H., Mochari, H., & Mosca, L. J. (2009). Waist circumference, body mass index, and their association with cardiometabolic and global risk. Journal of the cardiometabolic syndrome, 4(1), 12–19. doi:10.1111/j.1559-4572.2008.00029.x

K C Maki, K Kritsch, S Foley, I Soneru & M H Davidson (1997) Age-dependence of the relationship between adiposity and serum low density lipoprotein cholesterol in men., Journal of the American College of Nutrition, 16:6, 578-583.

II. Regresión múltiple

1.Datos

library(readxl)
death <- read_excel("death_small_cities.xlsx", 
    sheet = "data")
kable(death)

death1K	doctor100K	hosp100K	income1K	density1
8.0	78	284	9.1	109
9.3	68	433	8.7	144
7.5	70	739	7.2	113
8.9	96	1792	8.9	97
10.2	74	477	8.3	206
8.3	111	362	10.9	124
8.8	77	671	10.0	152
8.8	168	636	9.1	162
10.7	82	329	8.7	150
11.7	89	634	7.6	134
8.5	149	631	10.8	292
8.3	60	257	9.5	108
8.2	96	284	8.8	111
7.9	83	603	9.5	182
10.3	130	686	8.7	129
7.4	145	345	11.2	158
9.6	112	1357	9.7	186
9.3	131	544	9.6	177
10.6	80	205	9.1	127
9.7	130	1264	9.2	179
11.6	140	688	8.3	80
8.1	154	354	8.4	103
9.8	118	1632	9.4	101
7.4	94	348	9.8	117
9.4	119	370	10.4	88
11.2	153	648	9.9	78
9.1	116	366	9.2	102
10.5	97	540	10.3	95
11.9	176	680	8.9	80
8.4	75	345	9.6	92
5.0	134	525	10.3	126
9.8	161	870	10.4	108
9.8	111	669	9.7	77
10.8	114	452	9.6	60
10.1	142	430	10.7	71
10.9	238	822	10.3	86
9.2	78	190	10.7	93
8.3	196	867	9.6	106
7.3	125	969	10.5	162
9.4	82	499	7.7	95
9.4	125	925	10.2	91
9.8	129	353	9.9	52
3.6	84	288	8.4	110
8.4	183	718	10.4	69
10.8	119	540	9.2	57
10.1	180	668	13.0	106
9.0	82	347	8.8	40
10.0	71	345	9.2	50
11.3	118	463	7.8	35
11.3	121	728	8.2	86
12.8	68	383	7.4	57
10.0	112	316	10.4	57
6.7	109	388	8.9	94

2. Correlaciones

#Matriz de correlación:
cor(death)

##               death1K  doctor100K   hosp100K    income1K    density1
## death1K     1.0000000  0.11576504 0.11059019 -0.17199239 -0.27760696
## doctor100K  0.1157650  1.00000000 0.29562836  0.43328796 -0.01993791
## hosp100K    0.1105902  0.29562836 1.00000000  0.02750354  0.18661628
## income1K   -0.1719924  0.43328796 0.02750354  1.00000000  0.12874370
## density1   -0.2776070 -0.01993791 0.18661628  0.12874370  1.00000000

3. Modelos de regresión múltiple

A. Modelos de regresión múltiple y B. sus cálculos estadísticos:

a. Modelo 1 (Variable dependiente: death1K, variables independientes: doctor100K, hosp100K, income1K)

modRM1 <- lm(death1K~ doctor100K + hosp100K + income1K, 
              data = death)
summary(modRM1)

## 
## Call:
## lm(formula = death1K ~ doctor100K + hosp100K + income1K, data = death)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -5.7478 -1.0013  0.1878  1.2049  3.1652 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 11.9402252  2.0688314   5.771 5.25e-07 ***
## doctor100K   0.0094482  0.0070414   1.342   0.1858    
## hosp100K     0.0002713  0.0007231   0.375   0.7092    
## income1K    -0.4124113  0.2370644  -1.740   0.0882 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.646 on 49 degrees of freedom
## Multiple R-squared:  0.07681,    Adjusted R-squared:  0.02029 
## F-statistic: 1.359 on 3 and 49 DF,  p-value: 0.2663

b.Modelo 2 (Variable dependiente: death1K, variables independientes: interacción de doctor100K con hosp100K, income1K)

modRM2 <- lm(death1K~ doctor100K*hosp100K + income1K, 
              data = death)
summary(modRM2)

## 
## Call:
## lm(formula = death1K ~ doctor100K * hosp100K + income1K, data = death)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -5.8346 -1.0294  0.1467  1.2956  3.0475 
## 
## Coefficients:
##                       Estimate Std. Error t value Pr(>|t|)    
## (Intercept)          1.307e+01  2.851e+00   4.587 3.24e-05 ***
## doctor100K          -1.806e-03  2.057e-02  -0.088    0.930    
## hosp100K            -1.620e-03  3.326e-03  -0.487    0.628    
## income1K            -4.104e-01  2.387e-01  -1.719    0.092 .  
## doctor100K:hosp100K  1.759e-05  3.018e-05   0.583    0.563    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.657 on 48 degrees of freedom
## Multiple R-squared:  0.0833, Adjusted R-squared:  0.006907 
## F-statistic:  1.09 on 4 and 48 DF,  p-value: 0.3719

c.Modelo 3 (Variable dependiente: death1K, variables independientes: doctor100K, density e income1K)

modRM3 <- lm(death1K~ density1 + income1K +doctor100K, 
              data = death)
summary(modRM3)

## 
## Call:
## lm(formula = death1K ~ density1 + income1K + doctor100K, data = death)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -5.7852 -0.6675  0.4010  0.9348  2.7495 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 12.565900   1.978836   6.350 6.74e-08 ***
## density1    -0.008580   0.004746  -1.808   0.0768 .  
## income1K    -0.359140   0.230990  -1.555   0.1264    
## doctor100K   0.009284   0.006504   1.428   0.1598    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.596 on 49 degrees of freedom
## Multiple R-squared:  0.132,  Adjusted R-squared:  0.0789 
## F-statistic: 2.485 on 3 and 49 DF,  p-value: 0.07163

d. Modelo 4 (Variable dependiente: death1K, variables independientes: interacción de doctor100K con hosp100K con income1K con density)

modRM4 <- lm(death1K~ income1K * density1 * doctor100K *hosp100K,
              data = death)
summary(modRM4)

## 
## Call:
## lm(formula = death1K ~ income1K * density1 * doctor100K * hosp100K, 
##     data = death)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -5.4478 -0.4731  0.2674  1.0079  2.3411 
## 
## Coefficients:
##                                         Estimate Std. Error t value
## (Intercept)                            1.756e+02  1.153e+02   1.522
## income1K                              -1.767e+01  1.230e+01  -1.437
## density1                              -1.375e+00  1.150e+00  -1.196
## doctor100K                            -1.398e+00  9.621e-01  -1.454
## hosp100K                              -3.517e-01  2.334e-01  -1.507
## income1K:density1                      1.468e-01  1.229e-01   1.194
## income1K:doctor100K                    1.502e-01  1.018e-01   1.475
## density1:doctor100K                    1.104e-02  9.224e-03   1.197
## income1K:hosp100K                      3.657e-02  2.438e-02   1.500
## density1:hosp100K                      3.054e-03  2.251e-03   1.357
## doctor100K:hosp100K                    2.998e-03  1.837e-03   1.631
## income1K:density1:doctor100K          -1.200e-03  9.792e-04  -1.226
## income1K:density1:hosp100K            -3.204e-04  2.357e-04  -1.359
## income1K:doctor100K:hosp100K          -3.125e-04  1.910e-04  -1.636
## density1:doctor100K:hosp100K          -2.491e-05  1.742e-05  -1.430
## income1K:density1:doctor100K:hosp100K  2.628e-06  1.820e-06   1.444
##                                       Pr(>|t|)
## (Intercept)                              0.136
## income1K                                 0.159
## density1                                 0.239
## doctor100K                               0.154
## hosp100K                                 0.140
## income1K:density1                        0.240
## income1K:doctor100K                      0.149
## density1:doctor100K                      0.239
## income1K:hosp100K                        0.142
## density1:hosp100K                        0.183
## doctor100K:hosp100K                      0.111
## income1K:density1:doctor100K             0.228
## income1K:density1:hosp100K               0.182
## income1K:doctor100K:hosp100K             0.110
## density1:doctor100K:hosp100K             0.161
## income1K:density1:doctor100K:hosp100K    0.157
## 
## Residual standard error: 1.705 on 37 degrees of freedom
## Multiple R-squared:  0.2516, Adjusted R-squared:  -0.05179 
## F-statistic: 0.8293 on 15 and 37 DF,  p-value: 0.6408

C. Tabla comparativa de los estadísticos obtenidos para cada modelo:

library(tab)
library(sjPlot)
modRM1 <- lm(death1K~ doctor100K + hosp100K + income1K, 
              data = death)
tab_model(modRM1,modRM2,modRM3,modRM4)

	death1K			death1K			death1K			death1K
Predictors	Estimates	CI	p	Estimates	CI	p	Estimates	CI	p	Estimates	CI	p
(Intercept)	11.94	7.78 – 16.10	<0.001	13.07	7.34 – 18.81	<0.001	12.57	8.59 – 16.54	<0.001	175.60	-58.10 – 409.30	0.136
doctor100K	0.01	-0.00 – 0.02	0.186	-0.00	-0.04 – 0.04	0.930	0.01	-0.00 – 0.02	0.160	-1.40	-3.35 – 0.55	0.154
hosp100K	0.00	-0.00 – 0.00	0.709	-0.00	-0.01 – 0.01	0.628				-0.35	-0.82 – 0.12	0.140
income1K	-0.41	-0.89 – 0.06	0.088	-0.41	-0.89 – 0.07	0.092	-0.36	-0.82 – 0.11	0.126	-17.67	-42.59 – 7.25	0.159
doctor100K:hosp100K				0.00	-0.00 – 0.00	0.563				0.00	-0.00 – 0.01	0.111
density1							-0.01	-0.02 – 0.00	0.077	-1.38	-3.71 – 0.95	0.239
income1K:density1										0.15	-0.10 – 0.40	0.240
income1K:doctor100K										0.15	-0.06 – 0.36	0.149
density1:doctor100K										0.01	-0.01 – 0.03	0.239
income1K:hosp100K										0.04	-0.01 – 0.09	0.142
density1:hosp100K										0.00	-0.00 – 0.01	0.183
income1K:density1:doctor100K										-0.00	-0.00 – 0.00	0.228
income1K:density1:hosp100K										-0.00	-0.00 – 0.00	0.182
income1K:doctor100K:hosp100K										-0.00	-0.00 – 0.00	0.110
density1:doctor100K:hosp100K										-0.00	-0.00 – 0.00	0.161
income1K:density1:doctor100K:hosp100K										0.00	-0.00 – 0.00	0.157
Observations	53			53			53			53
R² / R² adjusted	0.077 / 0.020			0.083 / 0.007			0.132 / 0.079			0.252 / -0.052

D. Gráficas de los residuales:

spreadLevelPlot(modRM1)

## 
## Suggested power transformation:  -2.246865

spreadLevelPlot(modRM2)

## 
## Suggested power transformation:  -0.8045761

spreadLevelPlot(modRM3)

## 
## Suggested power transformation:  -0.4243738

spreadLevelPlot(modRM4)

## 
## Suggested power transformation:  1.742708

Discusión de los resultados:

El modelo 4 mostró un valor de coeficiente de determinación (R2) de 0.2516 (25.2%), este fue el valor obtenido de coeficiente de determinación mayor, por otro lado, este modelo tiene un valor de p de 0.6408, el cual es muy grande. Debido a esto el modelo más ideal es el modelo 3. El modelo 3, a pesar de tener un coeficiente de determinación menor: 0.132 (13.2%) mostró un valor de p más bajo: 0.07163.

4. Comparación de coeficientes de los modelos

library(ggplot2)
library(coefplot)

modRM1 <- lm(death1K~ doctor100K + hosp100K + income1K, 
              data = death)
modRM2 <- lm(death1K~ doctor100K*hosp100K + income1K, 
              data = death)
modRM3 <- lm(death1K~ density1 + income1K +doctor100K, 
              data = death)
modRM4 <- lm(death1K~ income1K * density1 * doctor100K *hosp100K,
              data = death)

multiplot(modRM1, modRM2, modRM3, modRM4, pointSize = 2, intercept=FALSE)

Esta gráfica indica que el mejor modelo es el que emplea la variable del ingreso.

5. Selección de modelos por pasos

modNulodeath <- lm(death1K ~ 1, data = death)
modFulldeath <- lm(death1K~ doctor100K + hosp100K + income1K*doctor100K + hosp100K + income1K+density1 +density1* income1K +density1* hosp100K+density1* +doctor100K* hosp100K+doctor100K* income1K +income1K *hosp100K, 
              data = death)

deathstep <- step(modNulodeath,
                scope = list(lower=modNulodeath, upper=modFulldeath,
                             direction="both"))

## Start:  AIC=54.87
## death1K ~ 1
## 
##              Df Sum of Sq    RSS    AIC
## + density1    1   11.0765 132.65 52.624
## <none>                    143.73 54.875
## + income1K    1    4.2517 139.48 55.283
## + doctor100K  1    1.9262 141.80 56.159
## + hosp100K    1    1.7578 141.97 56.222
## 
## Step:  AIC=52.62
## death1K ~ density1
## 
##              Df Sum of Sq    RSS    AIC
## <none>                    132.65 52.624
## + hosp100K    1    3.9272 128.72 53.031
## + income1K    1    2.7132 129.94 53.529
## + doctor100K  1    1.7471 130.91 53.921
## - density1    1   11.0765 143.73 54.875

death1 <- lm(formula = death1K ~ density1 , data = death)
summary(death1)

## 
## Call:
## lm(formula = death1K ~ density1, data = death)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -5.7119 -1.0315  0.1011  1.0413  2.9696 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 10.387997   0.569350  18.245   <2e-16 ***
## density1    -0.009782   0.004740  -2.064   0.0442 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.613 on 51 degrees of freedom
## Multiple R-squared:  0.07707,    Adjusted R-squared:  0.05897 
## F-statistic: 4.259 on 1 and 51 DF,  p-value: 0.04416

Utilizando el procediento step el valor estadístico de AIC más bajo que se obtuvo fue de 52.62. Por lo tanto, el modelo que mantiene el mínimo valor de AIC es el que utiliza la densidad (density1). Al comparar este modelo con los previamente analizados, se encontró que este modelo es más ideal ya que presenta un valor de p (0.04416) menor al resto de los modelos y el valor del coeficiente de determinación se mantiene cercano al resto de los modelos.

6. Selección de modelos con otros procedimientos

library(leaps)
modCp <- regsubsets(death1K ~ density1 + income1K + hosp100K + density1*income1K, data = death, nbest = 2)
plot(modCp, scale="adjr2")

plot(modCp, scale="Cp")

Según la gráfica de R-cuadrado, los modelos que utilizan la variable hosp100K y la variable de la densidad interactuando con el ingreso son mejores, es decir, presentan un valor de R-cuadrado mayor. En el caso del valor de Cp, el más bajo (0.88) está en el intercepto y en la interacción de densidad con ingreso. Ambos modelos sugeridos por estas últimas dos gráficas son similares entre sí ya que proponen la interacción de densidad con ingreso como un buen modelo, pero difieren con el valor estadístico de AIC, el cual sugiere que el mejor modelo es utilizando la variable de densidad.

7. Conclusiones

El valor estadístico de AIC sugirió que el mejor modelo era utilizando la variable de la densidad mientras que los otros dos métodos (R-cuadrado y Mallow’s Cp) tenían en común como mejor modelo la interacción entre densidad y el ingreso. Estos tres métodos concuerdan con los valores obtenidos en la matriz realizada previamente. A la misma vez, son similares a los resultados obtenidos en la gráfica de comparación de coeficientes de los modelos.

La interacción entre densidad e ingreso, parece ser una variable más adecuada para el modelo debido a que las personas con menor ingreso tienen menor acceso a servicios de salud y debido a esto son más propensos a sufrir de consecuencias graves. En conclusión los factores que pueden incidir de forma más significativa sobre la tasa muertes en las ciudades pequeñas de E.E.U.U. los son la interacción entre la densidad y el ingreso de las personas. Por lo tanto, las variables que son más determinantes de la variable dependiente lo son la densidad y el ingreso, siendo el mejor modelo de regresión múltiple el que incluye a estas dos variables.

III. Regresión logística

1. Datos

library(readxl)
Titanic <- read_excel("Titanic.xlsx")

2. División de los datos en “training” y “test”

set.seed(123)
training.samples <- Titanic$Survived %>% 
  createDataPartition(p = 0.8, list = FALSE)
train.data  <- Titanic[training.samples, ]
test.data <- Titanic[-training.samples, ]

3. Modelo y coeficientes

train.data$Age<- as.numeric(as.character(train.data$Age))

## Warning: NAs introduced by coercion

library(MASS)

## 
## Attaching package: 'MASS'

## The following object is masked from 'package:EnvStats':
## 
##     boxcox

## The following object is masked from 'package:dplyr':
## 
##     select

modelo <- glm(Survived ~., data = train.data, family = binomial)

summary(modelo)

## 
## Call:
## glm(formula = Survived ~ ., family = binomial, data = train.data)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -2.6208  -0.5753  -0.3451   0.4926   2.5372  
## 
## Coefficients:
##              Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  5.287599   0.542251   9.751  < 2e-16 ***
## Pclass      -1.149689   0.140959  -8.156 3.46e-16 ***
## Sexmale     -3.575132   0.229566 -15.573  < 2e-16 ***
## Age         -0.031339   0.007951  -3.941 8.10e-05 ***
## Nfamily     -0.224549   0.075009  -2.994  0.00276 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 1115.56  on 829  degrees of freedom
## Residual deviance:  650.48  on 825  degrees of freedom
##   (218 observations deleted due to missingness)
## AIC: 660.48
## 
## Number of Fisher Scoring iterations: 5

Según el análisis estadístico las variables más determinantes de la variable respuesta lo son la clase, el sexo (masculino) y la edad. La variable del número de familiares demuestra ser significante, aunque no del mismo modo que las variables ya mencionadas.

4. Gráficas binomiales con las variables más significativas

library(ggplot2)
#gráfica con curva logística
ggplot(Titanic, aes(x=Nfamily, y=Survived, na.rm = TRUE)) +
  geom_point() +
  geom_smooth(method = "glm", 
    method.args = list(family = "binomial"), 
    se = TRUE)

5. Cálculo del porcentaje de predicción utilizando los datos de “test”

#Predicciones
test.data$Age<- as.numeric(as.character(test.data$Age))

## Warning: NAs introduced by coercion

probabilities <- modelo %>% predict(test.data, type = "response", na.action=na.pass)
predicted.classes <- ifelse(probabilities > 0.5, "1", "0")
# Precisión de la predicción
observed.classes <- test.data$Survived
mean(predicted.classes == observed.classes, na.rm = TRUE)

## [1] 0.875

6. Repetición de los pasos 3 y 5 (Modelo y sus coeficientes junto a los porcentajes de predicción)

a. Modelo utilizando la variable de edad solamente

Coeficientes y valores estadísticos del modelo a

modelo1 <- glm(Survived ~Age, data = train.data, family = binomial)

summary(modelo1)

## 
## Call:
## glm(formula = Survived ~ Age, family = binomial, data = train.data)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -1.0810  -1.0178  -0.9721   1.3430   1.5009  
## 
## Coefficients:
##              Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.229995   0.160786   -1.43    0.153
## Age         -0.006306   0.004927   -1.28    0.201
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 1115.6  on 829  degrees of freedom
## Residual deviance: 1113.9  on 828  degrees of freedom
##   (218 observations deleted due to missingness)
## AIC: 1117.9
## 
## Number of Fisher Scoring iterations: 4

Porcentaje de predicción del modelo a

probabilities <- modelo1 %>% predict(test.data, type = "response", na.action=na.pass)
predicted.classes <- ifelse(probabilities > 0.5, "1", "0")
# Prediction accuracy
observed.classes <- test.data$Survived
mean(predicted.classes == observed.classes, na.rm = TRUE)

## [1] 0.625

b. Modelo utilizando la variable de sexo solamente

Coeficientes y valores estadísticos del modelo b

modelo2 <- glm(Survived ~Sex, data = train.data, family = binomial)

summary(modelo2)

## 
## Call:
## glm(formula = Survived ~ Sex, family = binomial, data = train.data)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -1.7794  -0.5358  -0.5358   0.6780   2.0060  
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)   1.3533     0.1278   10.59   <2e-16 ***
## Sexmale      -3.2217     0.1708  -18.86   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 1381.43  on 1047  degrees of freedom
## Residual deviance:  910.32  on 1046  degrees of freedom
## AIC: 914.32
## 
## Number of Fisher Scoring iterations: 4

Porcentaje de predicción del modelo b

probabilities <- modelo2 %>% predict(test.data, type = "response", na.action=na.pass)
predicted.classes <- ifelse(probabilities > 0.5, "1", "0")
# Prediction accuracy
observed.classes <- test.data$Survived
mean(predicted.classes == observed.classes, na.rm = TRUE)

## [1] 0.8697318

c. Modelo utilizando la variable de clase solamente

Coeficientes y valores estadísticos del modelo c

modelo3 <- glm(Survived ~Pclass, data = train.data, family = binomial)

summary(modelo3)

## 
## Call:
## glm(formula = Survived ~ Pclass, family = binomial, data = train.data)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -1.3202  -0.7728  -0.7728   1.0410   1.6458  
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  1.02229    0.18616   5.491 3.99e-08 ***
## Pclass      -0.69267    0.07892  -8.777  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 1381.4  on 1047  degrees of freedom
## Residual deviance: 1301.1  on 1046  degrees of freedom
## AIC: 1305.1
## 
## Number of Fisher Scoring iterations: 4

Porcentaje de predicción del modelo c

probabilities <- modelo3 %>% predict(test.data, type = "response", na.action=na.pass)
predicted.classes <- ifelse(probabilities > 0.5, "1", "0")
# Prediction accuracy
observed.classes <- test.data$Survived
mean(predicted.classes == observed.classes, na.rm = TRUE)

## [1] 0.6666667

d. Modelo utilizando la variable de número de familiares solamente

Coeficientes y valores estadísticos del modelo d

modelo4 <- glm(Survived ~Nfamily, data = train.data, family = binomial)

summary(modelo4)

## 
## Call:
## glm(formula = Survived ~ Nfamily, family = binomial, data = train.data)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -1.0963  -0.9488  -0.9488   1.4081   1.4247  
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept) -0.56482    0.07327  -7.709 1.27e-14 ***
## Nfamily      0.03709    0.03895   0.952    0.341    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 1381.4  on 1047  degrees of freedom
## Residual deviance: 1380.5  on 1046  degrees of freedom
## AIC: 1384.5
## 
## Number of Fisher Scoring iterations: 4

Porcentaje de predicción del modelo d

probabilities <- modelo4 %>% predict(test.data, type = "response", na.action=na.pass)
predicted.classes <- ifelse(probabilities > 0.5, "1", "0")
# Prediction accuracy
observed.classes <- test.data$Survived
mean(predicted.classes == observed.classes, na.rm = TRUE)

## [1] 0.6360153

e. Modelo utilizando las variables de sexo y edad

Coeficientes y valores estadísticos del modelo e

modelo5 <- glm(Survived ~Age+Sex, data = train.data, family = binomial)

summary(modelo5)

## 
## Call:
## glm(formula = Survived ~ Age + Sex, family = binomial, data = train.data)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -1.8413  -0.5654  -0.5600   0.6459   1.9776  
## 
## Coefficients:
##              Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  1.418187   0.232126   6.110 9.99e-10 ***
## Age          0.001484   0.006503   0.228    0.819    
## Sexmale     -3.221818   0.191373 -16.835  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 1115.56  on 829  degrees of freedom
## Residual deviance:  734.62  on 827  degrees of freedom
##   (218 observations deleted due to missingness)
## AIC: 740.62
## 
## Number of Fisher Scoring iterations: 4

Porcentaje de predicción del modelo e

probabilities <- modelo5 %>% predict(test.data, type = "response", na.action=na.pass)
predicted.classes <- ifelse(probabilities > 0.5, "1", "0")
# Prediction accuracy
observed.classes <- test.data$Survived
mean(predicted.classes == observed.classes, na.rm = TRUE)

## [1] 0.8703704

Discusión

En el análisis de los datos (pasajeros del RMS Titanic) se emplearon seis modelos distintos. El primer modelo incluyó todas las variables (Sexo, Edad, Clase y Número de familiares). Al realizar el cálculo del porcentaje de predicción se obtuvo un 87.5%. Luego se realizaron el resto de los modelos (a-Edad, b-Sexo, c-Clase, d-Número de familiares & e-Edad y Sexo), para estos se obtuvieron los porcentajes de predicción 62.5%, 87.0% , 66.7%, 63.6% y 87.0%, respectivamente. Esto nos indica que utilizar la variable de Sexo crea un mejor modelo ya que el resto de las variables no mostraron buenos modelos, pero esto sólo si se utilizan las varibles de forma individual. Al combinar las varibles de Sexo y Edad combinadas se mejoró ligeramente el modelo. El modelo que resultó ser mejor fue el modelo original (emplea todas las variables) ya que su porcentaje de predicción fue mayor. Al comparar el modelo original con los que solo emplearon una variable se puede observar que las variables presentan cierto grado de interacción ya que al ser combinadas en un modelo su porcentaje de predicción mejora. A pesar de que el modelo original es el mejor, se debe destacar que la variable de Sexo es buena predictora, es decir, es buena determinante de la variable respuesta y que influyó significativamente en la supervivencia de los pasajeros. Según el análisis estadístico del modelo el sexo masculino resultó mostrar mayor supervivencia lo que lleva a un cuestionamiento de la veracidad del código de conducta “mujeres y niños priemro”. En conclusión el mejor modelo es el original (emplea todas las variables) debido a que tiene el menor valor de AIC (660.48) y el porcentaje de predicción mayor.

Informe - Modelos Empíricos de Regresión

etand

9/5/2019

I. Regresión lineal simple y polinomial:

1. Datos

2. Selección de la varible dependiente y varible independiente

3. Gráfica de puntos

4. Cálculo de estadísticos

Estimadores de los coeficientes de regresión para un modelo lineal simple y para un modelo polinomial

5. Gráfica de residuales

6. Gráfica de líneas de regresión con intervalo de confianza

7. Explicación y discusión de resultados

Referencias (regresión lineal)

II. Regresión múltiple

1.Datos

2. Correlaciones

3. Modelos de regresión múltiple

A. Modelos de regresión múltiple y B. sus cálculos estadísticos:

a. Modelo 1 (Variable dependiente: death1K, variables independientes: doctor100K, hosp100K, income1K)

b.Modelo 2 (Variable dependiente: death1K, variables independientes: interacción de doctor100K con hosp100K, income1K)

c.Modelo 3 (Variable dependiente: death1K, variables independientes: doctor100K, density e income1K)

d. Modelo 4 (Variable dependiente: death1K, variables independientes: interacción de doctor100K con hosp100K con income1K con density)

C. Tabla comparativa de los estadísticos obtenidos para cada modelo:

D. Gráficas de los residuales:

Discusión de los resultados:

4. Comparación de coeficientes de los modelos

5. Selección de modelos por pasos

6. Selección de modelos con otros procedimientos

7. Conclusiones

III. Regresión logística

1. Datos

2. División de los datos en “training” y “test”

3. Modelo y coeficientes

4. Gráficas binomiales con las variables más significativas

5. Cálculo del porcentaje de predicción utilizando los datos de “test”

6. Repetición de los pasos 3 y 5 (Modelo y sus coeficientes junto a los porcentajes de predicción)

a. Modelo utilizando la variable de edad solamente

Coeficientes y valores estadísticos del modelo a

Porcentaje de predicción del modelo a

b. Modelo utilizando la variable de sexo solamente

Coeficientes y valores estadísticos del modelo b

Porcentaje de predicción del modelo b

c. Modelo utilizando la variable de clase solamente

Coeficientes y valores estadísticos del modelo c

Porcentaje de predicción del modelo c

d. Modelo utilizando la variable de número de familiares solamente

Coeficientes y valores estadísticos del modelo d

Porcentaje de predicción del modelo d

e. Modelo utilizando las variables de sexo y edad

Coeficientes y valores estadísticos del modelo e

Porcentaje de predicción del modelo e

Discusión