R Tarea Estadística
#Defino mi directorio de trabajo
setwd("~/Documents/Estadistica")
library(pacman)
p_load(haven,dplyr,ggplot2,MASS)
#Importar la base de datos
Datos <- read_sav("Datos hombro.sav")
#Definimos como factor la variable sexoN
Datos$sexoN <- factor(Datos$sexoN,
levels = c(1, 2),
labels = c("Hombre", "Mujer"))
#Para ver frecuencia por sexo
table (Datos$sexoN)
##
## Hombre Mujer
## 50 30
#Lo que sigue es resumir por sexo la longitud máxima de la clavícula derecha:
res_lmcd <- Datos %>%
group_by(sexoN) %>%
summarise(
n = sum(!is.na(LMCD)),
media = mean(LMCD, na.rm = TRUE),
sd = sd(LMCD, na.rm = TRUE)
) %>%
mutate(across(c(media, sd), ~round(.x, 2)))
res_lmcd
## # A tibble: 2 × 4
## sexoN n media sd
## <fct> <int> <dbl> <dbl>
## 1 Hombre 35 151. 10.2
## 2 Mujer 25 131. 7.56
#Y resumir por sexo la circunferencia en punto medio de la clavícula derecha:
res_c12cd <- Datos %>%
group_by(sexoN) %>%
summarise(
n = sum(!is.na(C12CD)),
media = mean(C12CD, na.rm = TRUE),
sd = sd(C12CD, na.rm = TRUE)
) %>%
mutate(across(c(media, sd), ~round(.x, 2)))
res_c12cd
## # A tibble: 2 × 4
## sexoN n media sd
## <fct> <int> <dbl> <dbl>
## 1 Hombre 46 35.9 2.33
## 2 Mujer 29 30.3 2.83
#Ahora unas gráficas comparativas:
library(ggplot2)
ggplot(Datos, aes(x = LMCD, fill = sexoN)) +
geom_density(alpha = 0.5) +
labs(
title = "Distribución de longitud máxima de la clavícula derecha",
x = "Longitud máxima de la clavícula derecha",
y = "Densidad"
) +
theme_minimal()
## Warning: Removed 20 rows containing non-finite outside the scale range
## (`stat_density()`).
#Box plot:
ggplot(Datos, aes(x = sexoN, y = LMCD, fill = sexoN)) +
geom_boxplot(alpha = 0.7) +
labs(
title = "Longitud máxima de la clavícula derecha por sexo",
x = "Sexo",
y = "Longitud máxima de la clavícula derecha)"
) +
theme_minimal()
## Warning: Removed 20 rows containing non-finite outside the scale range
## (`stat_boxplot()`).
#Siguiendo el punto 1 del ejercicio. Comparaciones por sexo de ambas variables:
#Alternativa no paramétrica: Prueba u de Mann-Whitney
wilcox.test(LMCD ~ sexoN, data = Datos)
## Warning in wilcox.test.default(x = DATA[[1L]], y = DATA[[2L]], ...): cannot
## compute exact p-value with ties
##
## Wilcoxon rank sum test with continuity correction
##
## data: LMCD by sexoN
## W = 831.5, p-value = 3.584e-09
## alternative hypothesis: true location shift is not equal to 0
#Esto quiere decir que hay una diferencia estadísticamente significativa, y se debería rechazar la hipótesis nula.
wilcox.test(C12CD ~ sexoN, data = Datos)
## Warning in wilcox.test.default(x = DATA[[1L]], y = DATA[[2L]], ...): cannot
## compute exact p-value with ties
##
## Wilcoxon rank sum test with continuity correction
##
## data: C12CD by sexoN
## W = 1258, p-value = 1.154e-10
## alternative hypothesis: true location shift is not equal to 0
#Esto quiere decir que hay una diferencia estadísticamente significativa, y se debería rechazar la hipótesis nula.
#Homogeneidad de varianza
var.test(LMCD ~ sexoN, data = Datos)
##
## F test to compare two variances
##
## data: LMCD by sexoN
## F = 1.8064, num df = 34, denom df = 24, p-value = 0.1344
## alternative hypothesis: true ratio of variances is not equal to 1
## 95 percent confidence interval:
## 0.8287493 3.7477045
## sample estimates:
## ratio of variances
## 1.806424
var.test(C12CD ~ sexoN, data = Datos)
##
## F test to compare two variances
##
## data: C12CD by sexoN
## F = 0.67915, num df = 45, denom df = 28, p-value = 0.2424
## alternative hypothesis: true ratio of variances is not equal to 1
## 95 percent confidence interval:
## 0.3353139 1.3032174
## sample estimates:
## ratio of variances
## 0.679154
#Siguiendo el punto 2 del ejercicio.
#Construir una función discriminante exponiendo el modelo, punto de corte y % de clasificación correcta usando la variable LMCD.
#Ejemplo de aplicación:
#Perlita, una estudiante de licenciatura en Antropología Física está decidiendo acerca de su tema de tesis.
#Está planteado la posibilidad de estimar el sexo mediante la medición de la longitud máxima de la clavícula derecha.
#Para conocer la viabilidad de lo que plantea realizará una serie de pruebas estadisticas.
#Quitamos los NA de la base de datos
Hombro_sinNA <- na.omit(Datos[, c("sexoN", "LMCD")])
#Se ajusta el modelo con la base depurada
# El modelo discriminante que queremos ajustar es D = a*LMCD + b
lda <- lda(sexoN ~ LMCD, data = Hombro_sinNA) #Ajustar Disciminante Lineal
#Para obtener los coeficientes a y b:
a <- coef(lda)[1]
lda
## Call:
## lda(sexoN ~ LMCD, data = Hombro_sinNA)
##
## Prior probabilities of groups:
## Hombre Mujer
## 0.5833333 0.4166667
##
## Group means:
## LMCD
## Hombre 151.0263
## Mujer 131.3026
##
## Coefficients of linear discriminants:
## LD1
## LMCD 0.1089485
#Predicción de valores mediante modelo:
PREMO <- predict(lda)
table(Real = Hombro_sinNA$sexoN, Predicho = PREMO$class)
## Predicho
## Real Hombre Mujer
## Hombre 32 3
## Mujer 3 22
#Porcentaje de acierto:
mean(PREMO$class == Hombro_sinNA$sexoN)*100
## [1] 90
#El porcentaje de calsificación correcta es de 90%
# valor de D cuando LMCD=0 y C12CD=0 (intercepto):
b <- PREMO$x
b
## LD1
## 1 0.56564625
## 2 1.65513096
## 3 -0.85068387
## 4 -1.97441778
## 5 -2.86843900
## 6 -1.34095199
## 7 -2.32148823
## 8 -2.81175635
## 9 -1.69503615
## 10 -1.28647776
## 11 -0.63278693
## 12 -1.06858081
## 13 0.12985237
## 14 -1.77674587
## 15 -0.08804458
## 16 -1.41565452
## 17 -2.37596246
## 18 -0.95963234
## 19 -1.63915983
## 20 -1.63915983
## 21 -1.39542623
## 22 -0.68726117
## 23 -1.06858081
## 24 -0.19699305
## 25 -0.52383846
## 26 0.02090390
## 27 -0.63278693
## 28 -0.85068387
## 29 1.44016152
## 30 -0.19699305
## 31 0.34774931
## 32 0.67459472
## 33 -0.41488999
## 34 0.56564625
## 35 0.56564625
## 36 -0.07462270
## 37 1.21933708
## 38 -1.12305505
## 39 1.69500093
## 40 0.34774931
## 41 1.11038861
## 42 0.78354319
## 43 0.56564625
## 44 0.18432660
## 45 1.31274182
## 46 3.56172920
## 47 -1.06858081
## 48 0.20475893
## 49 1.21933708
## 50 1.11038861
## 51 0.78354319
## 52 1.38275978
## 53 -0.08887492
## 54 1.21933708
## 55 0.67459472
## 56 3.34383226
## 57 1.76407943
## 58 0.80306301
## 59 3.07146108
## 60 2.74461567
#Para el punto de corte:
centroide_H <- mean(PREMO$x[Hombro_sinNA$sexoN == "Hombre"])
centroide_M <- mean(PREMO$x[Hombro_sinNA$sexoN == "Mujer"])
cutoff <- mean(c(centroide_H, centroide_M))
cutoff
## [1] -0.1790721
#El punto de corte es: -0.1790721
#El modelo discriminante es:
cat("Función discriminante: D(x) = ", round(a, 4), " * LMCD + ", round(b, 4), "\n")
## Función discriminante: D(x) = 0.1089 * LMCD + 0.5656 1.6551 -0.8507 -1.9744 -2.8684 -1.341 -2.3215 -2.8118 -1.695 -1.2865 -0.6328 -1.0686 0.1299 -1.7767 -0.088 -1.4157 -2.376 -0.9596 -1.6392 -1.6392 -1.3954 -0.6873 -1.0686 -0.197 -0.5238 0.0209 -0.6328 -0.8507 1.4402 -0.197 0.3477 0.6746 -0.4149 0.5656 0.5656 -0.0746 1.2193 -1.1231 1.695 0.3477 1.1104 0.7835 0.5656 0.1843 1.3127 3.5617 -1.0686 0.2048 1.2193 1.1104 0.7835 1.3828 -0.0889 1.2193 0.6746 3.3438 1.7641 0.8031 3.0715 2.7446
#La función discriminante es:
#D(x) = 0.1089 * LMCD + 0.5656 1.6551 -0.8507 -1.9744 -2.8684 -1.341 -2.3215 -2.8118 -1.695 -1.2865 -0.6328 -1.0686 0.1299 -1.7767 -0.088 -1.4157 -2.376 -0.9596 -1.6392 -1.6392 -1.3954 -0.6873 -1.0686 -0.197 -0.5238 0.0209 -0.6328 -0.8507 1.4402 -0.197 0.3477 0.6746 -0.4149 0.5656 0.5656 -0.0746 1.2193 -1.1231 1.695 0.3477 1.1104 0.7835 0.5656 0.1843 1.3127 3.5617 -1.0686 0.2048 1.2193 1.1104 0.7835 1.3828 -0.0889 1.2193 0.6746 3.3438 1.7641 0.8031 3.0715 2.7446
cat("Punto de corte=", round(cutoff, 4), "Si D>",round(cutoff, 4), "es Hombre")
## Punto de corte= -0.1791 Si D> -0.1791 es Hombre
#El punto de corte es: -0.1791. Entonces si D> -0.1791 es hombre
tabla_clas <- table(Observado = Hombro_sinNA$sexoN,
Predicho = PREMO$class)
tabla_clas
## Predicho
## Observado Hombre Mujer
## Hombre 32 3
## Mujer 3 22
prop_clas <- sum(diag(tabla_clas)) / sum(tabla_clas)*100
#Siguiendo el punto 3 del ejercicio.
#Construir una función discriminante exponiendo el modelo, punto de corte y % de clasificación correcta usando las variables LMCD y C12CD.
#Ejemplo de aplicación:
#Perlita, la misma estudiante de Antropología Física del ejemplo anterior, no esta
#convencida con los resultados que obtuvo de sus pruebas estadisticas. Así que, decidió
#no solo estimar utilizar la longitud máxima de la clavícula derecha para estimar el sexo, sino también
#incluir la medición de la circunferencia en un punto medio de la clavícula derecha, y así realizar
#nuevamente la serie de pruebas estadisticas para conocer la viabilida de lo que plantea.
# Eliminar NAs de forma explícita
Hombro_sinNA2 <- na.omit(Datos[, c("sexoN", "LMCD","C12CD")])
#Ajustar el modelo con la base depurada
#El modelo discriminante que queremos ajustar es D=a*LMCD+b*C12CD+c
lda2 <- lda(sexoN ~ LMCD+C12CD, data = Hombro_sinNA2) #Ajustar Disciminante Lineal
# Obtener los coeficientes a y b
a1 <- coef(lda2)[1]
a2 <- coef(lda2)[2]
lda2
## Call:
## lda(sexoN ~ LMCD + C12CD, data = Hombro_sinNA2)
##
## Prior probabilities of groups:
## Hombre Mujer
## 0.5833333 0.4166667
##
## Group means:
## LMCD C12CD
## Hombre 151.0263 35.77151
## Mujer 131.3026 30.31206
##
## Coefficients of linear discriminants:
## LD1
## LMCD -0.0727066
## C12CD -0.2559868
# Predicción de valores mediante modelo
pred2 <- predict(lda2)
# Matriz de confusión
table(Real = Hombro_sinNA2$sexoN, Predicho = pred2$class)
## Predicho
## Real Hombre Mujer
## Hombre 34 1
## Mujer 1 24
# Porcentaje de acierto
mean(pred2$class == Hombro_sinNA2$sexoN)*100
## [1] 96.66667
#El porcentaje de calsificación correcta es de 96.66%
#valor de D cuando LMCD=0 y C12CD=0 (intercepto)
b2 <- predict(lda2,
newdata = data.frame(LMCD = 0, C12CD = 0))$x
b2
## LD1
## 1 18.95782
#Para el punto de corte:
centroide_H2 <- mean(pred2$x[Hombro_sinNA2$sexoN == "Hombre"])
centroide_M2 <- mean(pred2$x[Hombro_sinNA2$sexoN == "Mujer"])
cutoff2 <- mean(c(centroide_H2, centroide_M2))
cutoff2
## [1] 0.2359658
#El punto de corte es: 0.2359658
#El modelo discriminante es:
cat("Función discriminante: D(x) = ", round(a1, 4), " * LMCD + ", round(a2, 4),"*C12CD +",round(b2, 4), "\n")
## Función discriminante: D(x) = -0.0727 * LMCD + -0.256 *C12CD + 18.9578
#La función discriminante es: D(x) = -0.0727 * LMCD + -0.256 *C12CD + 18.9578
cat("Punto de corte=", round(cutoff, 4), "Si D>",round(cutoff, 4), "es Hombre")
## Punto de corte= -0.1791 Si D> -0.1791 es Hombre
#El punto de corte es: -0.1791. Es decir, si D> -0.1791 es hombre.
#Graficando:
ggplot(Hombro_sinNA2, aes(x = LMCD, y = C12CD, color = sexoN)) +
geom_point(size = 3) +
theme_minimal() +
labs(title = "Discriminante lineal: LMCD vs C12CD")
a1 <- lda2$scaling["LMCD", 1]
a2 <- lda2$scaling["C12CD", 1]
b2 <- predict(lda2, newdata = data.frame(LMCD = 0, C12CD = 0))$x - cutoff
LD1 <- pred2$x[, 1]
#Graficando
ggplot(data.frame(LD1, sexoN = Hombro_sinNA2$sexoN),
aes(x = LD1, fill = sexoN)) +
geom_density(alpha = 0.4) +
geom_vline(xintercept = cutoff, linetype = "dashed") +
scale_x_continuous(limits = c(-5, 5)) + # <-- escala de -5 a 5
theme_minimal() +
labs(title = "Distribución sobre la función discriminante",
x = "LD1 (puntuación discriminante)",
y = "Densidad")
#Siguiendo el punto 4 del ejercicio.
#Ajustar un modelo logístico, calcular % de clasificación correcta usando la variable LMCD .
#Ejemplo de aplicación para determinar la probabilidad de hombre a partir de un valor hipotético de LMCD.
#Frecuencia del sexo
table (Datos$sexoN)
##
## Hombre Mujer
## 50 30
#Cambiamos las categorías para predecir p(Hombre)
Datos$sexoN <- factor(Datos$sexoN,
levels = c(2, 1),
labels = c("Mujer", "Hombre"))
#Ajustamos el modelo logístico binario:
modelo <- glm(sexoN ~ LMCD,
data = Hombro_sinNA,
family = binomial(link = "logit"))
summary(modelo)
##
## Call:
## glm(formula = sexoN ~ LMCD, family = binomial(link = "logit"),
## data = Hombro_sinNA)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 41.72697 10.79320 3.866 0.000111 ***
## LMCD -0.29945 0.07717 -3.880 0.000104 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 81.503 on 59 degrees of freedom
## Residual deviance: 32.430 on 58 degrees of freedom
## AIC: 36.43
##
## Number of Fisher Scoring iterations: 6
#Probabilidades de ser "Hombre":
Hombro_sinNA$prob_Hombre <- predict(modelo, type = "response")
#Clasificación usando 0.5 como punto de corte
Hombro_sinNA$predicho <- ifelse(Hombro_sinNA$prob_Hombre >= 0.5, "Hombre", "Mujer")
#mayor 0.5 hombre y menor mujer.
table(Real = Hombro_sinNA$sexoN, Predicho = Hombro_sinNA$predicho)
## Predicho
## Real Hombre Mujer
## Hombre 3 32
## Mujer 22 3
#Porcentaje de acierto:
mean(Hombro_sinNA$sexoN == Hombro_sinNA$predicho) * 100
## [1] 10
#El porcentaje de clasificación correcta es del 10%.
#Gráfica de probabilidades predichas
ggplot(Hombro_sinNA, aes(x = prob_Hombre, fill = sexoN)) +
geom_density(alpha = 0.4) +
geom_vline(xintercept = 0.5, linetype = "dashed") +
theme_minimal() +
labs(title = "Probabilidades predichas de ser Hombre",
x = "P(Hombre)", y = "Densidad")
#Coeficientes del modelo:
coef(modelo)
## (Intercept) LMCD
## 41.7269684 -0.2994501
#Ecuación logistica:
cat("Ecuación logística:\nlogit(p) = ",
round(coef(modelo)[1], 4), " + ",
round(coef(modelo)[2], 4), "*LMCD\n")
## Ecuación logística:
## logit(p) = 41.727 + -0.2995 *LMCD
#La ecuación es logit(p) = 41.727 -0.2995 *LMCD
#Ahora, Perlita nos proporciona el siguiente valor para probarlo con la ecuación que obtuvimos,
#la longitud máxima de la clavícula derecha es: 116.4797.
#Entonces tenemos que logit(p) = 41.727 -0.2995 *116.4797 = 6.8413.
#Si calculamos el logaritmo tenemos que e elevado a −6.8413 es igual a 0.00107.
#Si lo sustituimos en la fórmula de probalidad, donde p=1/(1+0.00107), tenemos que
#p es igual a 0.999, es decir, según el modelo la probabilidad
#de que esa clavícula pertenezca a un hombre es de 99.9%. Sin embargo, debemos recordar que
#el procentaje de clasificación correcta del modelo es del 10%.
#Siguiendo el punto 5 del ejercicio.
#Ajustar un modelo logístico, calcular % de clasificación correcta usando % de clasificación correcta usando las variables LMCD y C12CD.
#Hacer un ejemplo de aplicación para valores hipotéticos de LMCD y C12CD.
# Ajustar modelo logístico binario
modelo2 <- glm(sexoN ~ LMCD + C12CD,
data = Hombro_sinNA2,
family = binomial(link = "logit"))
summary(modelo2)
##
## Call:
## glm(formula = sexoN ~ LMCD + C12CD, family = binomial(link = "logit"),
## data = Hombro_sinNA2)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 69.0885 22.9591 3.009 0.00262 **
## LMCD -0.3067 0.1117 -2.747 0.00602 **
## C12CD -0.8022 0.3405 -2.356 0.01847 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 81.503 on 59 degrees of freedom
## Residual deviance: 18.731 on 57 degrees of freedom
## AIC: 24.731
##
## Number of Fisher Scoring iterations: 8
#Probabilidades de ser "Hombre"
Hombro_sinNA2$prob_Hombre <- predict(modelo2, type = "response")
#Clasificación usando 0.5 como punto de corte
Hombro_sinNA2$predicho <- ifelse(Hombro_sinNA2$prob_Hombre >= 0.5, "Hombre", "Mujer")
#mayor 0.5 hombre y menor mujer
table(Real = Hombro_sinNA2$sexoN, Predicho = Hombro_sinNA2$predicho)
## Predicho
## Real Hombre Mujer
## Hombre 1 34
## Mujer 24 1
#Porcentaje de acierto
mean(Hombro_sinNA2$sexoN == Hombro_sinNA2$predicho) * 100
## [1] 3.333333
#El porcentaje de clasificación correcta es del 3.33%.
#Gráfica de probabilidades predichas:
ggplot(Hombro_sinNA2, aes(x = prob_Hombre, fill = sexoN)) +
geom_density(alpha = 0.4) +
geom_vline(xintercept = 0.5, linetype = "dashed") +
theme_minimal() +
labs(title = "Probabilidades predichas de ser Hombre",
x = "P(Hombre)", y = "Densidad")
#Para sacar los coeficientes del modelo:
coef(modelo2)
## (Intercept) LMCD C12CD
## 69.0884686 -0.3067252 -0.8022468
#Ecuación logistica
cat("Ecuación logística:\nlogit(p) = ",
round(coef(modelo2)[1], 4), " + ",
round(coef(modelo2)[2], 4), "*LMCD + ",
round(coef(modelo2)[3], 4), "*C12CD\n")
## Ecuación logística:
## logit(p) = 69.0885 + -0.3067 *LMCD + -0.8022 *C12CD
#La ecuación es logit(p) = 69.0885 -0.3067 *LMCD -0.8022 *C12CD
#Ahora, Perlita nos proporciona el siguiente valor para probarlo con la ecuación que obtuvimos,
#la longitud máxima de la clavícula derecha es: 116.4797, y,
#la circunferencia en el punto medio de la clavícula es: 28.42324
#Entonces tenemos que logit(p) = 69.0885 -0.3067 *116.4797 -0.8022 *28.42324 = 10.56305
#Si calculamos el logaritmo tenemos que e elevado a −10.56305 es igual a 0.0000259.
#Si lo sustituimos en la fórmula de probalidad, donde p=1/(1+0.0000259), tenemos que
#p es igual a 0.99997, es decir, según el modelo la probabilidad
#de que esa clavícula pertenezca a un hombre es de 99.997%. Sin embargo, debemos recordar que
#el procentaje de clasificación correcta del modelo es del 3.33%.