Modelos lineales para problemas no lineales

En esta sesión se usa un modelo lineal para modelar un fenómeno lineal: los precios de la vivienda en Taiwan. El conjunto de datos es Real estate valuation data set Data Set.

Citación: Yeh, I. C., & Hsu, T. K. (2018). Building real estate valuation models with comparative approach through case-based reasoning. Applied Soft Computing, 65, 260-271.

Se comienza con la lectura de los datos:

datos <- read.csv2("real_estate_valuation_dataset.csv")

A continuación se inspecciona visualmente el conjunto de datos con la función head():

head(datos)
##   No       X1   X2         X3 X4       X5       X6    Y
## 1  1 2012.917 32.0   84.87882 10 24.98298 121.5402 37.9
## 2  2 2012.917 19.5  306.59470  9 24.98034 121.5395 42.2
## 3  3 2013.583 13.3  561.98450  5 24.98746 121.5439 47.3
## 4  4 2013.500 13.3  561.98450  5 24.98746 121.5439 54.8
## 5  5 2012.833  5.0  390.56840  5 24.97937 121.5425 43.1
## 6  6 2012.667  7.1 2175.03000  3 24.96305 121.5125 32.1

La definición de las variables es la siguiente:

Para este problema se seleccionan solo las variables X2, X3, X4 y Y:

datos_an<-subset(datos,select = c("X2","X3","X4","Y"))
head(datos_an)
##     X2         X3 X4    Y
## 1 32.0   84.87882 10 37.9
## 2 19.5  306.59470  9 42.2
## 3 13.3  561.98450  5 47.3
## 4 13.3  561.98450  5 54.8
## 5  5.0  390.56840  5 43.1
## 6  7.1 2175.03000  3 32.1

Veamos un resumen de los datos con la función summary():

summary(datos_an)
##        X2               X3                X4               Y         
##  Min.   : 0.000   Min.   :  23.38   Min.   : 0.000   Min.   :  7.60  
##  1st Qu.: 9.025   1st Qu.: 289.32   1st Qu.: 1.000   1st Qu.: 27.70  
##  Median :16.100   Median : 492.23   Median : 4.000   Median : 38.45  
##  Mean   :17.713   Mean   :1083.89   Mean   : 4.094   Mean   : 37.98  
##  3rd Qu.:28.150   3rd Qu.:1454.28   3rd Qu.: 6.000   3rd Qu.: 46.60  
##  Max.   :43.800   Max.   :6488.02   Max.   :10.000   Max.   :117.50

Ahora se muestran

panel.cor <- function(x, y, digits = 2, prefix = "", cex.cor, ...)
{
    usr <- par("usr"); on.exit(par(usr))
    par(usr = c(0, 1, 0, 1))
    r <- abs(cor(x, y,use="na.or.complete"))
    txt <- format(c(r, 0.123456789), digits = digits)[1]
    txt <- paste0(prefix, txt)
    if(missing(cex.cor)) cex.cor <- 0.8/strwidth(txt)
    text(0.5, 0.5, txt, cex = cex.cor * r)
}

pairs(datos_an, lower.panel = panel.smooth, upper.panel = panel.cor)

Transformación de variables

Ahora se explora la pertinencia de una transformación logarítmica:

par(mfrow=c(1,3))
plot(datos_an$X2+1,datos_an$Y,log="xy", # la respuesta en escala logarítmica
     xlab="Edad de la casa [ln(1+años)]",
     ylab="Precio [ln ($/m2)]")
plot(datos_an$X3,datos_an$Y,log="xy",
     xlab="Distancia al MRT [ln m]",
     ylab="Precio [ln ($/m2)]")
plot(datos_an$X4+1,datos_an$Y,log="yx",
     xlab="Comercios de conveniencia cercanos [ln(1+n)]",
     ylab="Precio [ln ($/m2)]")

El siguiente código crea las variables transformadas:

datos_an$X2log<-log(datos_an$X2+1) #ln(1+x)
datos_an$X3log<-log(datos_an$X3) # ln(x)
datos_an$X4log<-log(datos_an$X4+1) #ln(1+x)
datos_an$Ylog<-log(datos_an$Y) # ln(x)

Ahora se muestran las dispersiones por pares incluyendo las variables transformadas:

pairs(datos_an, lower.panel = panel.smooth, upper.panel = panel.cor, 
      main="Dispersión por pares con variables transformadas")

¿Qué permite concluir el gráfico anterior sobre la transformación de variables?

Escalamiento de las variables:

El siguiente código centra las variables en su respectiva media y las divide por su desviación estándar:

datos_an_scale<-scale(datos_an,center = TRUE,scale = TRUE)
medias<-attr(datos_an_scale,"scaled:center") # se guarda la media para procesos posteriores
desv_est<-attr(datos_an_scale,"scaled:scale") # se guarda la sd para procesos posteriores
datos_an_scale<-as.data.frame(datos_an_scale) # los datos escalados se ponen en un dataframe

A continuación se muestra el gráfico de dispersión por pares de los datos escalados:

pairs(datos_an_scale, lower.panel = panel.smooth, upper.panel = panel.cor, 
      main="Dispersión por pares con variables transformadas y escaladas")

Calibración de un modelo lineal

Primero se separa el conjunto de datos en entrenamiento y validación:

set.seed(20200521)
p_vl<-0.2
N<-dim(datos_an)[1]
ix_vl<-sample(N,round(0.2*N),replace = FALSE)
datos_tr<-datos_an_scale[-ix_vl,]
datos_vl<-datos_an_scale[ix_vl,]

Ajuste de un modelo lineal con interacciones

El siguiente código muestra el ajuste de un modelo lineal con interacciones usando las variables en escala logarítmica y estandarizadas:

modelo1<-lm(Ylog~-1+X2log+X3log+X4log+X4log:X3log,data=datos_tr)
summary(modelo1)
## 
## Call:
## lm(formula = Ylog ~ -1 + X2log + X3log + X4log + X4log:X3log, 
##     data = datos_tr)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -4.2364 -0.2461  0.0753  0.3482  2.7206 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## X2log       -0.18169    0.03607  -5.037 7.84e-07 ***
## X3log       -0.55947    0.04571 -12.238  < 2e-16 ***
## X4log        0.21399    0.04528   4.726 3.41e-06 ***
## X3log:X4log  0.09110    0.02916   3.124  0.00194 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.6128 on 327 degrees of freedom
## Multiple R-squared:  0.6249, Adjusted R-squared:  0.6203 
## F-statistic: 136.2 on 4 and 327 DF,  p-value: < 2.2e-16
Gráficos diagnósticos

Residuales versus valores predichos:

plot(modelo1,which=1)

Gráfico cuantil-cuantil de los residuales estandarizados vs los de una distribución normal:

plot(modelo1,which=2)

Residuales estandarizados versus valores predichos:

plot(modelo1,which=3)

Distancia de Cook de cada observación:

plot(modelo1,which=4)

Residuales estandarizados versus leverage:

plot(modelo1,which=5)

Distancia de Cook vs Leverage:

plot(modelo1,which=6)

Ajuste de un modelo lineal sin interacciones

El siguiente código ajusta un modelo lineal sin interacciones:

modelo2<-lm(Ylog~-1+X2log+X3log+X4log,data=datos_tr)
summary(modelo2)
## 
## Call:
## lm(formula = Ylog ~ -1 + X2log + X3log + X4log, data = datos_tr)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -4.2759 -0.3005 -0.0131  0.3134  2.8060 
## 
## Coefficients:
##       Estimate Std. Error t value Pr(>|t|)    
## X2log -0.18504    0.03653  -5.065 6.84e-07 ***
## X3log -0.55692    0.04631 -12.025  < 2e-16 ***
## X4log  0.23645    0.04530   5.220 3.18e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.6209 on 328 degrees of freedom
## Multiple R-squared:  0.6137, Adjusted R-squared:  0.6102 
## F-statistic: 173.7 on 3 and 328 DF,  p-value: < 2.2e-16

Para comparar ambos modelos podemos usar la función anova():

anova(modelo2,modelo1)
## Analysis of Variance Table
## 
## Model 1: Ylog ~ -1 + X2log + X3log + X4log
## Model 2: Ylog ~ -1 + X2log + X3log + X4log + X4log:X3log
##   Res.Df    RSS Df Sum of Sq     F   Pr(>F)   
## 1    328 126.46                               
## 2    327 122.79  1    3.6646 9.759 0.001944 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

De acuerdo a la prueba anterior, la interacción sí le aporta al modelo.

Construcción de la función de predicción en la escala natural

El siguiente código calcula los parámetros del modelo en la escala natural de las variables. Es decir, calcula los parámetros de un modelo no lineal (las explicaciones se dan en la clase):

m1coef<-coef(modelo1) # extrae los coeficientes
b2s<-m1coef[1]
b3s<-m1coef[2]
b4s<-m1coef[3]
b34s<-m1coef[4]
b2p<-b2s*desv_est[8]/desv_est[5]
b3p<-b3s*desv_est[8]/desv_est[6]
b4p<-b4s*desv_est[8]/desv_est[7]
b34p<-b34s*desv_est[8]/(desv_est[6]*desv_est[7])
sigma<-summary(modelo1)$sigma
b0<-medias[8]-b2p*medias[5]-b3p*medias[6]-b4p*medias[7]+b34p*medias[6]*medias[7]#+((sigma)^2)/2
b3<-b3p-b34p*medias[7]
b4<-b4p-b34p*medias[6]
b<-c(b0,b2p,b3,b4,b34p) # ector con los coeficientes del modelo no lineal

En el siguiente código se define la función de predicción:

fn_pred=function(x,b){
  # Función de predicción derivada de un modelo lineal con variables con
  # transformación logarítimica y centradas
  y<-exp(b[1])*((x[1]+1)^b[2])*(x[2]^b[3])*((x[3]+1)^b[4])*((x[3]+1)^(b[5]*log(x[2])))
  return(y)
}

El siguiente código lleva a cabo una prueba de escritorio de la función de predicción:

x0<-datos_tr[1,1:3]*desv_est[1:3]+medias[1:3] # este paso multiplica por la sd y suma la media
x1<-(datos_an[-ix_vl,1:3])[1,]
(fn_pred(x0,b))
##         X2
## 1 49.86693
(fn_pred(x0,b))
##         X2
## 1 49.86693

Cálculo de los valores predichos en entrenamiento:

y_pred_tr<-apply(datos_an[-ix_vl,1:3],1,fn_pred,b=b)
plot(datos_an$Y[-ix_vl],y_pred_tr,ylab="Predichos",xlab="Observados",
     xlim=c(0,120),ylim=c(0,120),las=1,
     main="Predichos vs observados (Entrenamiento)")
abline(a=0,b=1,lwd=2,col="red")
R_tr<-cor(datos_an$Y[-ix_vl],y_pred_tr)
R_tr<-format(R_tr,digits = 2)
rmse_tr<-sqrt(mean(datos_an$Y[-ix_vl]-y_pred_tr)^2)
rmse_tr<-format(rmse_tr,digits = 2,nsmall = 2)
grid()
legend("topleft",legend=paste0(c("Cor: ","RMSE: "),c(R_tr,rmse_tr)), bty="n")

Cálculo de los valores predichos para el conjunto de validación:

y_pred_vl<-apply(datos_an[ix_vl,1:3],1,fn_pred,b=b)
plot(datos_an$Y[ix_vl],y_pred_vl,ylab="Predichos",xlab="Observados",
     xlim=c(0,120),ylim=c(0,120),las=1,
     main="Predichos vs observados (Validación)")
abline(a=0,b=1,lwd=2,col="red")
R_vl<-cor(datos_an$Y[ix_vl],y_pred_vl)
R_vl<-format(R_vl,digits = 2)
rmse_vl<-sqrt(mean(datos_an$Y[ix_vl]-y_pred_vl)^2)
rmse_vl<-format(rmse_vl,digits = 2,nsmall = 2)
grid()
legend("topleft",legend=paste0(c("Cor: ","RMSE: "),c(R_vl,rmse_vl)), bty="n")

Construcción de la superficile de respuesta

Primero se generara una superficie de respuesta que depende de las variables X2 y X3 y para la cual se considerará un valor fijo para X4.

El siguiente código crea los puntos en los que se evaluará la superficie de respuesta:

X2<-seq(0,43.8,length.out=100)
X3<-seq(23.4,6488,length.out = 100)
# X4<-0:10
newdata<-expand.grid(X2,X3)
names(newdata)<-c("X2","X3")

El código anterior produjo el dataframe newdata con los valores sobre los que se quiere evaluar la superficie de respuesta. A continuación se agrega un valor fijo para la variable X4 (se fija en 5) y se calcula la respuesta para todos los puntos:

X4_0<-5 # referencia para X4
newdata$X4<-X4_0
newdata$Y<-apply(newdata,1,fn_pred,b=b)

A continuación se representa los valores de la superficie de respuesta en una matriz z:

z<-matrix(newdata$Y,ncol=length(X2),nrow = length(X3))

Con la matriz anterior se puede visualizar la superficie de respuesta con la función image():

par(mar=c(6,8,4,2),mgp=c(5,1,0))
z_img<-t(z)
image(z_img,xaxt='n',yaxt='n',xlab="Distancia al MRT [m]",ylab="Edad (años)",col=gray.colors(8,rev=TRUE))
title(main=c("Representación de la superficie", "de respuesta en función de X2 y X3","(X4=5)"),
      sub="X4=5")
labelsX2<-format(X2, digits = 1,decimal.mark = ".",big.mark = " ")
labelsX3<-format(X3, digits = 1,decimal.mark = ".",big.mark = " ")
axis(2,at=seq(0,1,by=1/(length(X2)-1)),labels =labelsX2,las=2)
axis(1,at=seq(0,1,by=1/(length(X3)-1)),labels = labelsX3,las=2)

contour(X2,X3,z,las=1,
        ylab="Distancia al MRT [m]",
        xlab="Edad (años)",
        main=c("Representación de las curvas de nivel de la superficie", "de respuesta en función de X2 y X3 (X4=5)"))

A continuación se presenta la superficie de respuesta en 3D estática:

persp(X2,X3,z,xlab="Edad de la casa",ylab="Distancia al MRT",zlab="Precio",main="Superficie de respuesta para un modelo con dos variables",theta=120,shade=0.25,sub="X4=5")

El siguiente código presenta la superficie de respuesta 3D con plotly:

library(plotly)
Z<-z
p<- plot_ly(x=X2,y=X3,z = z) 
p<- add_surface(p)
p<-layout(p,title='Precio vs dis MRT y número de comercios (X4=5)',
          xaxis=list(title="Edad"),
          yaxis=list(title="X3"),
          scene = list(camera=list(eye = list(x=1.87, y=0.88, z=-0.64))))
p

Intervalos de predicción

Univariados

Se creará un intervalo de predicción en función de la variable distancia al MRT (X2) para una casa que tenga 3 años de antigüedad y tenga cinco tiendas alrededor. Para ello se dan los siguientes pasos:

Se crea un dataframe con los valores de las variables explicativas sobre los que se calcularán las predicciones y sus intervalos:

datos_inter<-data.frame(X2=3,X3=X3,X4=5)

Se aplican las mismas transformaciones que requiere el modelo. La primera es la tranformación logarítmica:

datos_inter$X2log<-log(datos_inter$X2+1) #ln(1+x)
datos_inter$X3log<-log(datos_inter$X3) # ln(x)
datos_inter$X4log<-log(datos_inter$X4+1) #ln(1+x)

La segunda es el escalamiento:

datos_inter_scaled<-scale(datos_inter,
                          center=medias[-c(4,8)], # se quitan los valores de Y y Ylog
                          scale = desv_est[-c(4,8)]# se quitan los valores de Y y Ylog
                          )
datos_inter_scaled<-as.data.frame(datos_inter_scaled)

Ahora se pueden utilizar los datos para hacer la predicción:

pred_int<-predict(modelo1, newdata = datos_inter_scaled, interval = 'prediction')
pred_int<-as.data.frame(pred_int)
ymin<-min(pred_int)
ymax<-max(pred_int)
plot(datos_inter_scaled$X3log,pred_int$fit,las=1,type="l",
     col="red",ylim=c(ymin,ymax),xlab="ln(X3) estandarizada",
     ylab="ln(Precio) estandarizado")
lines(datos_inter_scaled$X3log,pred_int$lwr,col="blue",lty=2)
lines(datos_inter_scaled$X3log,pred_int$upr,col="blue",lty=2)
legend("topright",col=c("red","blue"),lty=1:2,legend=c("Predicción","95%CI"))

pred_int_escalanatural<-exp(pred_int*desv_est[8]+medias[8])
ymin<-min(pred_int_escalanatural)
ymax<-max(pred_int_escalanatural)
plot(X3,pred_int_escalanatural$fit,las=1,type="l",
     col="red",ylim=c(ymin,ymax),xlab="Distancia al MRT [m]",
     ylab="Precio [$/m2]",
     main=c("Precios vs distancia al MRT","Edad: 3 años, Tiendas: 5"))
lines(X3,pred_int_escalanatural$lwr,col="blue",lty=2)
lines(X3,pred_int_escalanatural$upr,col="blue",lty=2)
grid()
legend("topright",col=c("red","blue"),lty=1:2,
       legend=c("Predicción","95%CI"))

Actividades:

  1. Genere un gif animado con las superficies de respuesta para los distintos valores de X4
  2. Genere un video que muestre rotando las diferentes superficies para los distintos valores de X4
  3. Cree un intervalo de predicción para los precios en función del número de tiendas asumiendo que las casas tienen 3 años y están a 1000 m del MRT
  4. Visualice los datos atípicos en un mapa

Solución al punto 4

library(rgdal)
library(leaflet)
library(leaflet.extras)
CD<-cooks.distance(modelo1)
h<-influence(modelo1)$hat
datos_aug<-datos
datos_aug$CD<-NA
datos_aug$h<-NA
datos_aug[-ix_vl,"CD"]<-CD
datos_aug[-ix_vl,"h"]<-h
popup<-paste0("Edad: ",format(datos_aug$X2[-ix_vl],digits = 2),
              "<br>",
              "Dist MRT: ",format(datos_aug$X2[-ix_vl],digits = 2),
              "<br>",
              "# Tiendas: ",format(datos_aug$X3[-ix_vl],digits = 2),
              "<br>",
              "Precio: ", format(datos_aug$Y,digits=2),
              "<br>",
              "h: ",format(datos_aug$h,digits=2))
colorh<-ifelse(datos_aug$h>0.04 & !is.na(datos_aug$h),"#FF6666","#FFFF66")
mapa<-leaflet()
mapa<-addTiles(mapa,urlTemplate="https://tiles.stadiamaps.com/tiles/osm_bright/{z}/{x}/{y}{r}.png",
               attribution="<a href='https://www.openstreetmap.org/copyright'>OpenStreetMap</a> contributors"
               )
mapa<-addCircleMarkers(mapa,lat=datos_aug$X5[-ix_vl],lng = datos_aug$X6[-ix_vl],popup = popup,
                       color=colorh,
                        options=markerOptions(opacity = 0.3))
mapa

¿Dónde es más caro el precio de la propiedad?

mapa_p<-leaflet(data=datos)
mapa_p<-addTiles(mapa_p,urlTemplate="https://tiles.stadiamaps.com/tiles/osm_bright/{z}/{x}/{y}{r}.png",
               attribution="<a href='https://www.openstreetmap.org/copyright'>OpenStreetMap</a> contributors"
               )
mapa_p<-addHeatmap(mapa_p,lat=~X5,lng=~X6,intensity = ~Y^(-0.8),max=0.5,blur=20,radius=12)

mapa_p