Carga y limpieza de datos

library(readxl)
setwd("~/Desktop")
datos_brutos <- read_excel("EncuestaSitiosComida.xlsx")

Numero de registros en la base de datos original

cat("Registros originales en el archivo:", nrow(datos_brutos), "\n")
## Registros originales en el archivo: 1200

Generacion de muestra aleatoria (290) y semilla 4793.Asi al correr el codigo la muestra sera siempre la misma 

set.seed(4793)
tam_muestra <- 290
idx <- sample.int(n = nrow(datos_brutos), size = tam_muestra, replace = FALSE) #generar 290 valores aleatorios de 1 a 1200 sin repetir 
muestra_290 <- datos_brutos[idx, ] #generacion de la muestra 
cat("Tamaño de la muestra usada en el análisis:", nrow(muestra_290))
## Tamaño de la muestra usada en el análisis: 290

Seleccion de variables numericas

vars_num <- select(muestra_290, where(is.numeric))

A.Medidas descriptivas y tablas

Creacion de funcion para encontrar la moda

moda_simple <- function(v) {
  v_ok <- v[!is.na(v)]
  if (!length(v_ok)) return(NA_real_)
  t <- table(v_ok)
  as.numeric(names(t)[which.max(t)])
}
# Variable que almacena los nombres de las variables numericas
nom_var <- colnames(vars_num)

#Creacion de  Dataframe con resumen estadistico

res_desc <- data.frame(
  Variable = nom_var,
  Media    = sapply(vars_num, function(x) mean(x, na.rm = TRUE)),
  Mediana  = sapply(vars_num, function(x) median(x, na.rm = TRUE)),
  Moda     = sapply(vars_num, moda_simple),
  DesvEst  = sapply(vars_num, function(x) sd(x, na.rm = TRUE)),
  CV_pct   = NA_real_,
  row.names = NULL
)

res_desc$CV_pct <- round(100 * res_desc$DesvEst / res_desc$Media, 2)

1)Resumen estadistico

cat(">>> Resumen básico de variables cuantitativas:\n")
## >>> Resumen básico de variables cuantitativas:
flextable(res_desc)

Variable

Media

Mediana

Moda

DesvEst

CV_pct

ID

570.131034

564.5

4.0

351.317693

61.62

V2_TiempoConoce_meses

32.420690

24.0

21.0

25.813992

79.62

V5_Atencion_1a4

2.444828

3.0

3.0

1.064718

43.55

V6_VisitasSemana

2.900000

2.5

2.0

1.927851

66.48

V7_PrecioAlmuerzo

22,782.700000

22,933.5

15,256.0

4,705.617244

20.65

V8_CompanerosAlmuerzo

1.910345

2.0

2.0

1.411359

73.88

V9_Edad

37.165517

35.0

32.0

9.973330

26.83

V11_Ingresos

4,001,613.458621

3,195,199.5

1,500,000.0

2,550,634.349144

63.74

V12_Gastos

2,399,030.955172

1,907,369.5

500,000.0

1,648,285.544614

68.71

V13_PesoKg

69.432414

69.7

66.6

9.287331

13.38

V14_EstaturaCm

166.466552

166.1

150.0

8.634289

5.19

Cuartiles, RIQ, P10 y P90

extra_cuartiles <- t(sapply(
  vars_num,
  function(x) {
    q <- quantile(x, probs = c(0.25, 0.5, 0.75, 0.10, 0.90), na.rm = TRUE)
    c(Q1 = q[1], Q2 = q[2], Q3 = q[3], P10 = q[4], P90 = q[5], RIQ = q[3] - q[1])
  }
))

cuartiles_df <- data.frame(
  Variable = rownames(extra_cuartiles),
  extra_cuartiles,
  row.names = NULL
)

2)Cuartiles, RIQ y percentiles 10 y 90:

flextable(cuartiles_df)

Variable

Q1.25.

Q2.50.

Q3.75.

P10.10.

P90.90.

RIQ.75.

ID

261.50

564.5

855.750

83.90

1,090.10

594.250

V2_TiempoConoce_meses

14.00

24.0

41.750

9.00

66.10

27.750

V5_Atencion_1a4

1.25

3.0

3.000

1.00

4.00

1.750

V6_VisitasSemana

1.00

2.5

4.000

1.00

6.00

3.000

V7_PrecioAlmuerzo

18,622.25

22,933.5

26,761.750

16,500.20

29,118.80

8,139.500

V8_CompanerosAlmuerzo

1.00

2.0

3.000

0.00

4.00

2.000

V9_Edad

30.00

35.0

42.000

26.00

52.00

12.000

V11_Ingresos

2,182,824.75

3,195,199.5

5,121,989.000

1,500,000.00

7,663,144.10

2,939,164.250

V12_Gastos

1,252,848.00

1,907,369.5

3,036,766.500

914,772.30

4,644,757.20

1,783,918.500

V13_PesoKg

63.30

69.7

75.575

56.27

82.33

12.275

V14_EstaturaCm

160.00

166.1

173.275

155.28

177.90

13.275

3)Tabla de frecuencias para la V7 Precio del almuerzo

precio <- as.numeric(as.character(muestra_290$V7_PrecioAlmuerzo))
precio <- precio[!is.na(precio)]
#creacion de las clases o intervalos,en total son 8 
clases_precio <- cut(
  x = precio,
  breaks = 8,
  include.lowest = TRUE,
  right = FALSE,
  dig.lab=8
)
# Creacion de la tabla 
tabla_precio <- as.data.frame(table(clases_precio))
names(tabla_precio) <- c("Clase", "Frecuencia")

# Creacion de las columnas de la tabla
tabla_precio$Prop_rel   <- tabla_precio$Frecuencia / sum(tabla_precio$Frecuencia)
tabla_precio$Prop_acum  <- cumsum(tabla_precio$Prop_rel)
tabla_precio$Porc_rel   <- round(100 * tabla_precio$Prop_rel, 2)
tabla_precio$Porc_acum  <- round(100 * tabla_precio$Prop_acum, 2)

Tabla de distribución de V7_PrecioAlmuerzo (8 clases):

flextable(tabla_precio)

Clase

Frecuencia

Prop_rel

Prop_acum

Porc_rel

Porc_acum

[12988.539,15314.625)

11

0.03793103

0.03793103

3.79

3.79

[15314.625,17622.25)

38

0.13103448

0.16896552

13.10

16.90

[17622.25,19929.875)

46

0.15862069

0.32758621

15.86

32.76

[19929.875,22237.5)

40

0.13793103

0.46551724

13.79

46.55

[22237.5,24545.125)

33

0.11379310

0.57931034

11.38

57.93

[24545.125,26852.75)

53

0.18275862

0.76206897

18.28

76.21

[26852.75,29160.375)

40

0.13793103

0.90000000

13.79

90.00

[29160.375,31486.461]

29

0.10000000

1.00000000

10.00

100.00

kable(tabla_precio)
Clase Frecuencia Prop_rel Prop_acum Porc_rel Porc_acum
[12988.539,15314.625) 11 0.0379310 0.0379310 3.79 3.79
[15314.625,17622.25) 38 0.1310345 0.1689655 13.10 16.90
[17622.25,19929.875) 46 0.1586207 0.3275862 15.86 32.76
[19929.875,22237.5) 40 0.1379310 0.4655172 13.79 46.55
[22237.5,24545.125) 33 0.1137931 0.5793103 11.38 57.93
[24545.125,26852.75) 53 0.1827586 0.7620690 18.28 76.21
[26852.75,29160.375) 40 0.1379310 0.9000000 13.79 90.00
[29160.375,31486.461] 29 0.1000000 1.0000000 10.00 100.00 
datatable(tabla_precio)

4)Tabla de de contingencia V4 atencion categoria y V1 Sitio Zona 

tab_at_zona <- with(muestra_290,
                    table(V4_Atencion_cat, V1_SitioZona))

tab_at_zona
##                V1_SitioZona
## V4_Atencion_cat Zona Centro Zona Norte Zona Sur
##       Buena              41         28       23
##       Excelente          23         23        9
##       Mala               32         18       23
##       Regular            32         25       13

Proporciones por fila (Atención -> distribución en zonas)

prop_filas <- prop.table(tab_at_zona, margin = 1)
round(prop_filas, 3)
##                V1_SitioZona
## V4_Atencion_cat Zona Centro Zona Norte Zona Sur
##       Buena           0.446      0.304    0.250
##       Excelente       0.418      0.418    0.164
##       Mala            0.438      0.247    0.315
##       Regular         0.457      0.357    0.186

5)Ingresos por sexo

ing_x_sexo <- muestra_290 |>
  group_by(V15_Sexo) |>
  summarise(
    Media  = mean(V11_Ingresos, na.rm = TRUE),
    SD     = sd(V11_Ingresos, na.rm = TRUE),
    CV_pct = round(100 * SD / Media, 2),
    numero = n(),
    .groups = "drop"
  )

flextable(ing_x_sexo)

V15_Sexo

Media

SD

CV_pct

numero

Femenino

4,269,366

2,616,195

61.28

148

Masculino

3,722,547

2,458,666

66.05

142

B.Visualizacion

6)Histogramas y poligonos de frecuencia

for (v in colnames(vars_num)[-1]) {
  x <- vars_num[[v]]
  hist_obj <- hist(
    x,
    breaks = 12,
    freq=FALSE,
    main   = paste("Histograma y polígono -", v),
    xlab   = v,
    col    = "lightblue",
    border = "gray40",
    ylab="Frecuencia relativa"
  )
  # polígono de frecuencias normalizado
  lines(hist_obj$mids, hist_obj$density, type = "b", col = "darkblue", lwd = 2)
}

7)Cajas y bigotes

valores_largos <- pivot_longer(
  data = vars_num,
  cols = dplyr::everything(),
  names_to = "Variable",
  values_to = "Valor"
)

ggplot(valores_largos, aes(x = Variable, y = Valor)) +
  geom_boxplot(fill = "#fdd0a2", colour = "#d94801", outlier.alpha = 0.6) +
  theme_bw(base_size = 11) +
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
  labs(
    title = "Comparación de cajas para variables numéricas",
    x = "Variable",
    y = "Valor observado"
  )

8)QQ-plots y normalidad (Ingresos y Edad)

par(mfrow = c(1, 2))
qqnorm(ing <- na.omit(muestra_290$V11_Ingresos),
       main = "QQ-plot: V11_Ingresos",
       pch = 19, col = "#1f78b4")
qqline(ing, col = "#1f78b4", lwd = 2)

qqnorm(edad <- na.omit(muestra_290$V9_Edad),
       main = "QQ-plot: V9_Edad",
       pch = 19, col = "#33a02c")
qqline(edad, col = "#33a02c", lwd = 2)

par(mfrow = c(1, 1))

sh_ing  <- shapiro.test(ing)
sh_edad <- shapiro.test(edad)

cat("\n>>> Shapiro-Wilk para V11_Ingresos:\n"); print(sh_ing)
## 
## >>> Shapiro-Wilk para V11_Ingresos:
## 
##  Shapiro-Wilk normality test
## 
## data:  ing
## W = 0.8482, p-value = 0.0000000000000003335
cat("\n>>> Shapiro-Wilk para V9_Edad:\n"); print(sh_edad)
## 
## >>> Shapiro-Wilk para V9_Edad:
## 
##  Shapiro-Wilk normality test
## 
## data:  edad
## W = 0.95579, p-value = 0.0000001079

C.Conteo y probabilidad

9)Combinatoria (5 con exactamente 2 de Zona Centro)

n_centro <- sum(muestra_290$v1_sitio_zona == "Zona Centro")
## Warning: Unknown or uninitialised column: `v1_sitio_zona`.
n_no_centro <- nrow(muestra_290) - n_centro

formas_5 <- choose(n_centro, 2) * choose(n_no_centro, 3)

cat("\nNúmero de subgrupos posibles (5 personas, 2 de Zona Centro):", formas_5, "\n")
## 
## Número de subgrupos posibles (5 personas, 2 de Zona Centro): 0

10)\(P(Recomienda = Sí | \text{ V5_Atencion_1a4} \geq 3)\)

#subset se usa seleccionar un subconjunto de una base de datos si se pide que se satisfaga determinada condicion ,en otras palbras se usa para filtrar bases de datos o dataframes en R
subset_buena_at <- subset(muestra_290, V5_Atencion_1a4 >= 3)

p_cond <- mean(subset_buena_at$V3_Recomienda == "Sí", na.rm = TRUE)

cat("P(Recomienda = Sí | V5_Atencion_1a4 >= 3):",
    round(p_cond, 4), "\n")
## P(Recomienda = Sí | V5_Atencion_1a4 >= 3): 0.5374

11)Bayes: P(Zona Norte | Recomienda = Sí)

p_rec   <- mean(muestra_290$V3_Recomienda == "Sí",        na.rm = TRUE)
p_norte <- mean(muestra_290$V1_SitioZona   == "Zona Norte", na.rm = TRUE)

p_rec_y_norte <- mean(
  muestra_290$V3_Recomienda == "Sí" &
    muestra_290$V1_SitioZona == "Zona Norte",
  na.rm = TRUE
)

p_norte_dado_rec <- p_rec_y_norte / p_rec

cat("P(Zona Norte | Recomienda = Sí):", round(p_norte_dado_rec, 4), "\n")
## P(Zona Norte | Recomienda = Sí): 0.2761

12)Binomial: \(P(X\geq 8)\) con \(p\) = proporción muestral

p_hat <- mean(muestra_290$v3_recomienda == "Sí", na.rm = TRUE)
## Warning: Unknown or uninitialised column: `v3_recomienda`.
p_bin_ge8 <- 1 - pbinom(q = 7, size = 10, prob = p_hat)

cat("P(X >= 8) para X ~ Bin(10, p̂):", round(p_bin_ge8, 4), "\n")
## P(X >= 8) para X ~ Bin(10, p̂): NaN

13)Poisson: \(P(X \geq 4)\) visitas en Zona Sur

lambda_sur <- with(muestra_290,
                   mean(V6_VisitasSemana[V1_SitioZona == "Zona Sur"], na.rm = TRUE))

#ppois hace referencia a la probabilidad acumulada para una variable aleatoria distribuida Poisson (discreta) con promedio de ocurrencia del evento de interes conocido
p_pois_ge4 <- 1 - ppois(q = 3, lambda = lambda_sur)

cat("P(X >= 4) para visitas en Zona Sur (Poisson):", round(p_pois_ge4, 4), "\n")
## P(X >= 4) para visitas en Zona Sur (Poisson): 0.5283

14)Hipergeométrica: 6 ‘Excelente’ en 20 (N=1200)

N_poblacion <- 1200
K_exc <- sum(datos_brutos$V4_Atencion_cat == "Excelente", na.rm = TRUE)
#calculo de la poblacion que NO califica la atencion de manera escelente.Como el total debe ser la suma de :
#califican excelente+califican de otra manera(no excelente)=total--> despejamos no excelente
N_no_exc <- N_poblacion - K_exc

p_hip <- dhyper(
  x = 6,
  m = K_exc,
  n = N_no_exc,
  k = 20
)

cat("P(6 'Excelente' en una muestra de 20 sin reemplazo):", round(p_hip, 6), "\n")
## P(6 'Excelente' en una muestra de 20 sin reemplazo): 0.110361

15)Normal: \(P(20000 < Precio < 24000)\)

mu_prec <- mean(muestra_290$v7_precio_almuerzo, na.rm = TRUE)
## Warning: Unknown or uninitialised column: `v7_precio_almuerzo`.
## Warning in mean.default(muestra_290$v7_precio_almuerzo, na.rm = TRUE): argument
## is not numeric or logical: returning NA
sd_prec <- sd(muestra_290$v7_precio_almuerzo, na.rm = TRUE)
## Warning: Unknown or uninitialised column: `v7_precio_almuerzo`.
p_norm_20_24 <- pnorm(24000, mean = mu_prec, sd = sd_prec) -pnorm(20000, mean = mu_prec, sd = sd_prec)

cat("P(20000 < Precio < 24000) asumiendo Normal:", round(p_norm_20_24, 4), "\n")
## P(20000 < Precio < 24000) asumiendo Normal: NA

16)Exponencial: \(P(T > 10 ~días)\)

lambda_diaria <- mean(muestra_290$v6_visitas_semana, na.rm = TRUE) / 7
## Warning: Unknown or uninitialised column: `v6_visitas_semana`.
## Warning in mean.default(muestra_290$v6_visitas_semana, na.rm = TRUE): argument
## is not numeric or logical: returning NA
p_exp_10 <- exp(-lambda_diaria * 10)

cat("P(T > 10 días) bajo modelo Exponencial:", round(p_exp_10, 4), "\n")
## P(T > 10 días) bajo modelo Exponencial: NA

D: Intervalos de confianza y pruebas de hipotesis

17)IC 95% para media de V7_PrecioAlmuerzo

x_precio <- as.numeric(as.character(muestra_290$V7_PrecioAlmuerzo))
x_precio <- x_precio[!is.na(x_precio)]

ic_mean_precio <- t.test(x_precio, conf.level = 0.95)
cat("\nIC 95% para la media de V7_PrecioAlmuerzo:\n")
## 
## IC 95% para la media de V7_PrecioAlmuerzo:
print(ic_mean_precio$conf.int)
## [1] 22238.84 23326.56
## attr(,"conf.level")
## [1] 0.95

18)IC 95% para proporción de recomendación

indicador_reco <- muestra_290$V3_Recomienda == "Sí"

ic_prop_reco <- prop.test(
  x = sum(indicador_reco, na.rm = TRUE),
  n = sum(!is.na(indicador_reco)),
  conf.level = 0.95
)

cat("\nIC 95% para proporción que recomienda (V3_Recomienda = 'Sí'):\n")
## 
## IC 95% para proporción que recomienda (V3_Recomienda = 'Sí'):
print(ic_prop_reco$conf.int)
## [1] 0.4038718 0.5212874
## attr(,"conf.level")
## [1] 0.95

19)IC 95% para varianza de V11_Ingresos

n_ing  <- sum(!is.na(muestra_290$V11_Ingresos))
s2_ing <- var(muestra_290$V11_Ingresos, na.rm = TRUE)
alfa   <- 0.05

chi_inf <- qchisq(alfa / 2,     df = n_ing - 1)
chi_sup <- qchisq(1 - alfa / 2, df = n_ing - 1)

ic_var_ing <- c(
  (n_ing - 1) * s2_ing / chi_sup,
  (n_ing - 1) * s2_ing / chi_inf
)

cat("\nIC 95% para la varianza de V11_Ingresos:\n")
## 
## IC 95% para la varianza de V11_Ingresos:
print(ic_var_ing)
## [1] 5562834712230 7711847066837

20)Prueba t: \(H0: \mu= 22000\) (Precio)

x_precio <- as.numeric(as.character(muestra_290$V7_PrecioAlmuerzo))
x_precio <- x_precio[!is.na(x_precio)]

prueba_t_22000 <- t.test(
  x_precio,
  mu = 22000,
  alternative = "two.sided",
  conf.level = 0.95
)

cat("\nPrueba t bilateral para H0: mu = 22000 (V7_PrecioAlmuerzo):\n")
## 
## Prueba t bilateral para H0: mu = 22000 (V7_PrecioAlmuerzo):
print(prueba_t_22000)
## 
##  One Sample t-test
## 
## data:  x_precio
## t = 2.8326, df = 289, p-value = 0.004943
## alternative hypothesis: true mean is not equal to 22000
## 95 percent confidence interval:
##  22238.84 23326.56
## sample estimates:
## mean of x 
##   22782.7

21)Chi-cuadrado independencia Atención vs Zona

tabla_indep <- with(muestra_290,
                    table(V4_Atencion_cat, V1_SitioZona))

chi_res <- chisq.test(tabla_indep)

cat("\nChi-cuadrado de independencia entre V4_Atencion_cat y V1_SitioZona:\n")
## 
## Chi-cuadrado de independencia entre V4_Atencion_cat y V1_SitioZona:
print(chi_res)
## 
##  Pearson's Chi-squared test
## 
## data:  tabla_indep
## X-squared = 7.3173, df = 6, p-value = 0.2925

22)Normalidad de Precio (KS + Shapiro)

# Normalidad de Precio (KS + Shapiro)

# 1) Vector numérico de precios sin NA
precio_sin_na <- as.numeric(as.character(muestra_290$V7_PrecioAlmuerzo))
precio_sin_na <- precio_sin_na[!is.na(precio_sin_na)]

# 2) Parámetros muestrales
mu_prec <- mean(precio_sin_na)
sd_prec <- sd(precio_sin_na)

# 3) Precio estandarizado
precio_standar <- (precio_sin_na - mu_prec) / sd_prec

# 4) Kolmogorov–Smirnov contra N(0,1)
ks_norm <- ks.test(precio_standar, "pnorm")
## Warning in ks.test.default(precio_standar, "pnorm"): ties should not be present
## for the one-sample Kolmogorov-Smirnov test
cat("\nKolmogorov–Smirnov para Precio estandarizado ~ Normal(0,1):\n")
## 
## Kolmogorov–Smirnov para Precio estandarizado ~ Normal(0,1):
print(ks_norm)
## 
##  Asymptotic one-sample Kolmogorov-Smirnov test
## 
## data:  precio_standar
## D = 0.080692, p-value = 0.0458
## alternative hypothesis: two-sided
# 5) Shapiro-Wilk en una submuestra (máx. 500 datos)
set.seed(4793)
sub_muestra_prec <- sample(
  x     = precio_sin_na,
  size  = min(500, length(precio_sin_na)),
  replace = FALSE
)

sw_prec <- shapiro.test(sub_muestra_prec)
cat("\nShapiro-Wilk para V7_PrecioAlmuerzo (submuestra):\n")
## 
## Shapiro-Wilk para V7_PrecioAlmuerzo (submuestra):
print(sw_prec)
## 
##  Shapiro-Wilk normality test
## 
## data:  sub_muestra_prec
## W = 0.95883, p-value = 0.0000002602

E.Regresion y asociacion

23)Regresión lineal simple: Precio ~ Edad

#lm en ingles linear model se usa para generar el modelo de regresion lineal simple y/o multiple
mod_simple <- lm(V7_PrecioAlmuerzo ~ V9_Edad, data = muestra_290)
cat("\nRegresión lineal simple: V7_PrecioAlmuerzo ~ V9_Edad\n")
## 
## Regresión lineal simple: V7_PrecioAlmuerzo ~ V9_Edad
print(summary(mod_simple))
## 
## Call:
## lm(formula = V7_PrecioAlmuerzo ~ V9_Edad, data = muestra_290)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -9666.6 -3992.4    64.8  3808.7  8875.4 
## 
## Coefficients:
##             Estimate Std. Error t value            Pr(>|t|)    
## (Intercept) 23840.61    1067.77  22.327 <0.0000000000000002 ***
## V9_Edad       -28.46      27.75  -1.026               0.306    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 4705 on 288 degrees of freedom
## Multiple R-squared:  0.00364,    Adjusted R-squared:  0.0001801 
## F-statistic: 1.052 on 1 and 288 DF,  p-value: 0.3059
r_pe <- cor(muestra_290$V7_PrecioAlmuerzo, muestra_290$V9_Edad, use = "complete.obs")
cat("Correlación (Precio, Edad):", round(r_pe, 4), "\n")
## Correlación (Precio, Edad): -0.0603
cat("R^2 simple:", round(summary(mod_simple)$r.squared, 4), "\n")
## R^2 simple: 0.0036

24)Regresión múltiple 

mod_multiple <- lm(
  V7_PrecioAlmuerzo ~ V9_Edad + V11_Ingresos + V6_VisitasSemana,
  data = muestra_290
)

cat("\nRegresión múltiple: Precio ~ Edad + Ingresos + VisitasSemana\n")
## 
## Regresión múltiple: Precio ~ Edad + Ingresos + VisitasSemana
print(summary(mod_multiple))
## 
## Call:
## lm(formula = V7_PrecioAlmuerzo ~ V9_Edad + V11_Ingresos + V6_VisitasSemana, 
##     data = muestra_290)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -9464.5 -3561.7    79.5  3329.0  9545.4 
## 
## Coefficients:
##                        Estimate     Std. Error t value             Pr(>|t|)    
## (Intercept)      20831.58939915  1105.92254315  18.836 < 0.0000000000000002 ***
## V9_Edad              7.24663380    24.71146305   0.293             0.769544    
## V11_Ingresos         0.00077573     0.00009665   8.026   0.0000000000000262 ***
## V6_VisitasSemana  -490.47664695   127.91114176  -3.835             0.000155 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 4134 on 286 degrees of freedom
## Multiple R-squared:  0.2363, Adjusted R-squared:  0.2283 
## F-statistic:  29.5 on 3 and 286 DF,  p-value: < 0.00000000000000022

25)Correlacion de Pearson entre ingresos y gastos 

#cor calcula la correlacion entre variables
r_ing_gas <- cor(muestra_290$V11_Ingresos,  muestra_290$V12_Gastos, use = "complete.obs")
cat("\nCorrelación de Pearson entre V11_Ingresos y V12_Gastos:", round(r_ing_gas, 4), "\n")
## 
## Correlación de Pearson entre V11_Ingresos y V12_Gastos: 0.9393

F.ANOVA

26)ANOVA de una vía: Precio ~ Zona 

anova_uno <- aov(V7_PrecioAlmuerzo ~ V1_SitioZona, data = muestra_290)
cat("\nANOVA de una vía: V7_PrecioAlmuerzo ~ V1_SitioZona\n")
## 
## ANOVA de una vía: V7_PrecioAlmuerzo ~ V1_SitioZona
print(summary(anova_uno))
##               Df     Sum Sq    Mean Sq F value              Pr(>F)    
## V1_SitioZona   2 4813692513 2406846257   435.7 <0.0000000000000002 ***
## Residuals    287 1585586410    5524691                                
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
cat("\nVerificación rápida de supuestos (ANOVA 1 vía):\n")
## 
## Verificación rápida de supuestos (ANOVA 1 vía):
print(shapiro.test(residuals(anova_uno)))
## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(anova_uno)
## W = 0.99748, p-value = 0.9365
print(bartlett.test(V7_PrecioAlmuerzo ~ V1_SitioZona, data = muestra_290))
## 
##  Bartlett test of homogeneity of variances
## 
## data:  V7_PrecioAlmuerzo by V1_SitioZona
## Bartlett's K-squared = 34.96, df = 2, p-value = 0.00000002561
ggplot(muestra_290, aes(x = V1_SitioZona, y = V7_PrecioAlmuerzo, fill = V1_SitioZona)) +
  geom_boxplot(alpha = 0.85, colour = "gray20") +
  scale_fill_manual(values = c("#1b9e77", "#d95f02", "#7570b3")) +
  theme_minimal(base_size = 11) +
  theme(legend.position = "none") +
  labs(
    title = "Precio del almuerzo por zona (ANOVA 1 vía)",
    x = "Zona",
    y = "Precio del almuerzo"
  )

ANOVA de dos vías: Precio ~ Zona * Atención 

anova_dos <- aov(V7_PrecioAlmuerzo ~ V1_SitioZona * V4_Atencion_cat, data = muestra_290)
cat("\nANOVA de dos vías: V7_PrecioAlmuerzo ~ Zona * Atención\n")
## 
## ANOVA de dos vías: V7_PrecioAlmuerzo ~ Zona * Atención
print(summary(anova_dos))
##                               Df     Sum Sq    Mean Sq F value
## V1_SitioZona                   2 4813692513 2406846257 458.020
## V4_Atencion_cat                3   61745817   20581939   3.917
## V1_SitioZona:V4_Atencion_cat   6   62979390   10496565   1.997
## Residuals                    278 1460861203    5254896        
##                                            Pr(>F)    
## V1_SitioZona                 < 0.0000000000000002 ***
## V4_Atencion_cat                           0.00918 ** 
## V1_SitioZona:V4_Atencion_cat              0.06616 .  
## Residuals                                            
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
cat("\nPruebas de supuestos (ANOVA 2 vías):\n")
## 
## Pruebas de supuestos (ANOVA 2 vías):
print(shapiro.test(residuals(anova_dos)))
## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(anova_dos)
## W = 0.99524, p-value = 0.5145
print(car::leveneTest(V7_PrecioAlmuerzo ~ V1_SitioZona * V4_Atencion_cat,
                      data = muestra_290))
## Levene's Test for Homogeneity of Variance (center = median)
##        Df F value    Pr(>F)    
## group  11  3.2164 0.0003849 ***
##       278                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
interaction.plot(
  x.factor     = muestra_290$V1_SitioZona,
  trace.factor = muestra_290$V4_Atencion_cat,
  response     = muestra_290$V7_PrecioAlmuerzo,
  col          = c("#1b9e77", "#d95f02", "#7570b3", "#e7298a"),
  lwd          = 2,
  xlab         = "Zona",
  ylab         = "Precio promedio",
  trace.label  = "Categoría de atención",
  main         = "Interacción Zona x Atención",
  bty="l"
)