Pauta de solución realizada con los datos de PAULA ANDREA GUTIÉRREZ.

Si alguien desea que genere esta pauta con sus datos, avíseme, y encantado lo haré.

Pregunta 1

Descargue del aula virtual el archivo de datos con su nombre: Nombre – Pregunta 1.csv. El archivo consiste en un conjunto de terrenos agrícolas, con cuatro variables:

Los valores son ficticios; de modo que no se preocupe si no parecen guardar relación con la realidad. Realice una regresión con las cuatro variables, y con su modelo, prediga el precio de un terreno de 10.000 hectáreas, con lluvia 500 mm en el lugar de Curanilahue.

Respuesta

Como se desea predecir el precio, éste es la variable dependiente, y las otras variables son las predictivas. La variable “Lugar” es categórica, por lo cual se debe generar variables dummy de ella. Sus valores son:

\[\{Rengo, \ Curanilahue, \ Perquilauquen\}\] Al realizar la regresión lineal con R, incluyendo la variable categórica “Lugar”, el software automáticamente elige uno de los valores como valor de referencia y crea una variable dummy para cada uno de los valores restantes. Para cada conjunto de datos es distinto. Por ejemplo, podría elegir “Rengo” como valor de referencia y entonces crea las siguientes dos variables dummy:

\[Curanilahue = \left\{ \begin{matrix} 1 & cuando \ Lugar = "Curanilahue" \\ 0 & en \ caso \ contrario \end{matrix} \right.\]

\[Perquilauquen = \left\{ \begin{matrix} 1 & cuando \ Lugar = "Perquilauquen" \\ 0 & en \ caso \ contrario \end{matrix} \right.\]

Entonces el modelo de regresión lineal es:

\[Precio = \beta_0 + \beta_1 Hectáreas + \beta_2 Lluvia + \gamma_1 Curanilahue + \gamma_2 Perquilauquen + \varepsilon\] Primero, se hace un gráfico con cada una de las variables predictoras. Con las continuas se hace un gráfico de dispersión, y con la variable categórica se hace un gráfico de columnas.

library(ggplot2)
df_data = as.data.frame(read.table(paste0(ruta.datos.estudiantes, nombre, " - Pregunta 1.csv"), header = T, sep = "|"))
ggplot(data = df_data, mapping = aes(x=Hectareas, y=Precio)) + geom_point() + geom_smooth(method = "lm")
## `geom_smooth()` using formula = 'y ~ x'

Se observa una relación clara entre Hectáreas y Precio.

ggplot(data = df_data, mapping = aes(x=Lluvia, y=Precio)) + geom_point() + geom_smooth(method = "lm")
## `geom_smooth()` using formula = 'y ~ x'

Podría haber una relación tenue entre Lluvia y Precio.

library(plyr)
library(dplyr)
## 
## Attaching package: 'dplyr'
## The following objects are masked from 'package:plyr':
## 
##     arrange, count, desc, mutate, rename, summarise, summarize
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
df_lugares = summarise(group_by(df_data, Lugar), mean(Precio))
names(df_lugares) = c("Lugar", "Precio")
ggplot(data = df_lugares, mapping = aes(x=reorder(Lugar, Precio), y=Precio)) + geom_col()

Podría haber una relación tenue entre Lugar y Precio.

El siguiente código R realiza la regresión y predice el valor del precio para los valores entregados en el data frame df_nuevo.

mod = lm(data = df_data, formula = Precio ~ Hectareas + Lluvia + Lugar)
df_nuevo = data.frame(Hectareas = 10000, Lluvia = 500, Lugar = "Curanilahue")
predict.lm(mod, newdata = df_nuevo)
##        1 
## 28070.95

El precio predicho es 28.071 UF.

Pregunta 2

Descargue del aula virtual el archivo de datos con su nombre: Nombre – Pregunta 2.csv. El archivo consiste en una serie de tiempo para 50 periodos con dos variables:

Genere un modelo para predecir el valor de y en el periodo t = 60.

El puntaje de esta pregunta se asigna de acuerdo al valor absoluto del error porcentual de su predicción relativo a la predicción del mejor modelo posible. \[Puntaje=30(1-\varepsilon)\] donde \(\varepsilon\) es el valor absoluto del error porcentual.

Respuesta

Primero se hace el gráfico de la serie de tiempo.

df_data = as.data.frame(read.table(paste0(ruta.datos.estudiantes, nombre, " - Pregunta 2.csv"), header = T, sep = "|"))
ggplot(data = df_data, mapping = aes(x=t, y=y)) + geom_line()

Se debe verificar si la serie de tiempo es estacionaria o no, para saber que método usar para estimar el modelo. Se usa la prueba Dickey-Fuller. A cada estudiante le sale un resultado distinto, porque si bien todos los conjuntos de datos fueron generados como un AR(1) estacionario, sus parámetros \(\alpha_0\) y \(\alpha_1\) son distintos. A continuación se muestra el método en cada uno de los distintos escenarios que pueden surgir.

library(tseries)
## Registered S3 method overwritten by 'quantmod':
##   method            from
##   as.zoo.data.frame zoo
adf.test(df_data$y, k=0)
## Warning in adf.test(df_data$y, k = 0): p-value smaller than printed p-value
## 
##  Augmented Dickey-Fuller Test
## 
## data:  df_data$y
## Dickey-Fuller = -5.1279, Lag order = 0, p-value = 0.01
## alternative hypothesis: stationary
p = adf.test(df_data$y, k=0)$p.value
## Warning in adf.test(df_data$y, k = 0): p-value smaller than printed p-value

El valor p es 0.01, menor que 0,05. Esto indica que la serie es estacionaria y entonces se procede a realizar el método Box-Jenkins. Se prueban modelos ARMA(p,q) con distintos valores de p y q. `

if (estacionario){ # Box-Jenkins
  df_resultados = data.frame(p = c(0,0,1,1,1,2,2,2), q = c(1,2,0,1,2,0,1,2), AIC = rep(NA, 8))
  for (j in 1:nrow(df_resultados)){
    df_resultados$AIC[j] = arima(df_data$y, order = c(df_resultados$p[j], 0, df_resultados$q[j]))$aic
  }
  menor.aic = min(df_resultados$AIC)
  p = df_resultados[df_resultados$AIC==menor.aic, "p"]
  q = df_resultados[df_resultados$AIC==menor.aic, "q"]
  df_resultados  
}
##   p q      AIC
## 1 0 1 305.4945
## 2 0 2 307.4694
## 3 1 0 305.5274
## 4 1 1 307.5238
## 5 1 2 308.2974
## 6 2 0 307.5267
## 7 2 1 306.8378
## 8 2 2 308.0270

Los valores de p y q del modelo con menor AIC son p=0 y q=1. Por consiguiente, se realiza la predicción con el modelo ARMA con estos valores.

if (estacionario){ # Box-Jenkins
  mod = arima(df_data$y, order = c(p, 0, q))
  df_nuevo = data.frame(t = 60)
  predict(mod, newdata = df_nuevo)$pred
}
## Time Series:
## Start = 51 
## End = 51 
## Frequency = 1 
## [1] 6.038072

La predicción en el tiempo t=60 es 6.04.

Si la serie es no estacionaria, se debe determinar si es de tendencia determinista o estocástica. Incluso, es posible que sea estacionaria pero la prueba Dickey-Fuller estándar puede no detectarla cuando hay pocos datos. Entonces se procede con la prueba Dickey-Fuller completo.

Se realiza una regresión al modelo: \[\Delta y_t = \alpha + \beta t + \gamma y_{t-1} + \varepsilon_t\]

n = nrow(df_data)
df_data$y_1[2:n] = df_data$y[1:(n-1)]
df_data$Dy = df_data$y - df_data$y_1
mod = lm(data = df_data, formula = Dy ~ t + y_1)
tau = abs(summary(mod)$coefficients[2,3])
tau
## [1] 2.192045

Se tiene: \[\tau = 2.19\]

Si \(\tau > 1,95\), \(\beta \ne 0\) y entonces el modelo es determinista. En este caso se realiza la predicción con el modelo:

\[y = \beta_0 + \beta_1t + \varepsilon\]

mod = lm(data = df_data, formula = y ~ t)
predict.lm(mod, newdata = df_nuevo)
##      1 
## 9.3097

En este caso, la predicción es 9.31.

Si \(\tau \le 1,95\), \(\beta = 0\) y entonces el modelo es estocástico. En este caso se realiza una nueva regresión con el modelo:

\[\Delta y_t = \alpha + \gamma y_{t-1} + \varepsilon_t\]

mod = lm(data = df_data, formula = Dy ~ y_1)
tau = abs(summary(mod)$coefficients[2,3])

\[\tau = 4.46\]

Si ahora \(\tau > 1,95\), \(\gamma \ne 0\) y entonces el modelo es AR(1). En este caso se realiza la predicción con el modelo:

\[y = \alpha_0 + \alpha_1 y_{t-1} + \varepsilon_t\]

mod = arima(df_data$y, order = c(1,0,0))
predict(mod, newdata = df_nuevo)$pred
## Time Series:
## Start = 51 
## End = 51 
## Frequency = 1 
## [1] 6.093616

En este caso, la predicción es 6.09.

En cambio, si ahora \(\tau \le 1,95\), \(\gamma = 0 0\) y entonces el modelo es de raíz unitaria. En este caso se realiza la predicción con el modelo:

\[y_t = \alpha_0 + y_{t-1} + \varepsilon_t\] Este es una arima(0,1,0) porque al diferenciarlo una vez se vuelve estacionario y sin rezagos en \(y\) y \(\varepsilon\).

Ninguno de los conjuntos de datos de la evaluación es de raíz unitaria, de modo que no se muestra la predicción para ese caso.

Pregunta 3

Descargue del aula virtual el archivo de datos con su nombre: Nombre – Pregunta 3.csv.

El archivo consiste en un conjunto de clientes con las siguientes características:

Realice un análisis de conglomerados con las variables Edad, Educación e Ingreso, usando el algoritmo k-medios. Determine la cantidad óptima de conglomerados con el método del codo y de la silueta. Genere al menos 3 conglomerados, y máximo 6. Escriba sus conglomerados en un archivo Excel con la siguiente estructura:

Conglomerado Centro Edad Centro Educación Centro Ingreso Descripción del Conglomerado

Conglomerado 1
Conglomerado 2
…
Conglomerado n

Suba su archivo Excel con su tabla de conglomerados a la plataforma.

Respuesta

Se cargan los datos y se implementa el método del codo y el método de la silueta para encontrar el número óptimo de conglomerados.

# Paquete
library(flexclust)
library(cluster)

# Cargar datos
df_data = as.data.frame(read.table(paste0(ruta.datos.estudiantes, nombre, " - Pregunta 3.csv"), header = T, sep = "|"))

# Preparar matriz para k-medios
X = as.matrix(df_data[, 2:4])

# Crear y poblar data frame para guardar sumas de cuadrados de distancias dentro de centroides
df_totwithinss = data.frame(  k = 2:10
                              , tot.withinss = NA
)
for (k in 2:10){
  modelo = kcca(X, k = k, family = kccaFamily("kmeans"))
  df_totwithinss[k-1, "tot.withinss"] = info(modelo, "distsum")
}

# Gráfico
ggplot(df_totwithinss, aes(y = tot.withinss, x = k)) +
  ggtitle("Método del Codo") +
  ylab("Suma de Cuadrados de Distancias") +
  geom_line(colour = "red") +
  scale_x_continuous(breaks = 1:10)

# Crear y poblar data frame para guardar siluetas medias para cada k
d = dist(X)
df_silinfoavgwidth = data.frame(  k = 2:10
                                , silinfo.avg.width = NA
                                )
for (k in 2:10){
  modelo = kcca(X, k = k, family = kccaFamily("kmeans"))
  cl = clusters(modelo)
  sil = silhouette(cl, d)
  df_silinfoavgwidth[k-1, "silinfo.avg.width"] = mean(sil[, "sil_width"])
}

# Gráfico
ggplot(df_silinfoavgwidth, aes(y = silinfo.avg.width, x = k)) +
  ggtitle("Método de la Silueta") +
  ylab("Ancho Medio de Silueta") +
  geom_line(colour = "red") +
  scale_x_continuous(breaks = 1:10)

Al observar a ambos gráficos, se postula k=3 como un número óptimo de conglomerados. Se generan los conglomerados con esta cantidad elegida. Luego se agrega el número de conglomerado al data frame original de datos, de manera de asignar a cada Rut a un conglomerado. Finalmente, se exporta a Excel a fin de analizar los conglomerados con una tabla dinámica.

library(openxlsx)
k = 3
# Generar conglomerados
rownames(df_data) = df_data$Rut
X = as.matrix(df_data[, 2:4])
modelo = kcca(X, k = k, family = kccaFamily("kmeans"))
df_conglomerados = data.frame(Rut = names(modelo@cluster), Conglomerado = modelo@cluster)

# Agregar número de conglomerado a data frame de datos
df_data = merge(df_data, df_conglomerados, all.x = T)

# Exportar a Excel
write.xlsx(df_data, file = "Personas Asignadas a Conglomerados.xlsx")

La tabla Excel y tabla dinámica con el resumen de conglomerados se muestran a continuación.

Fig. 1 Asignación de Conglomerados a Personas

Fig. 2 Resumen de Conglomerados

——————————————————————————–