library(readxl)
datos <- read_excel("datos_vivienda.xlsx")
attach(datos)

Analisís exploratorio

Análisis Exploratorio de Precio de Vivienda

summary(datos$precio_millon)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   240.0   251.2   305.0   332.1   395.0   480.0
Desviacion_Estandar_Precio <-sd(datos$precio_millon)
Desviacion_Estandar_Precio
## [1] 82.14423

Según estos datos podemos afirmar que la variable Precio de Vivienda (millones de pesos), cuenta con valores entre 240.0 siendo su mínimo y 480.0 siendo su máximo, además cuenta con una media de 332.1.0 y una desviación estándar de 82.14423 millones de pesos.

plot(datos$precio_millon, 
     main = "Distribución de Precio de Vivienda",
     xlab = "Índice del dato",
     ylab = "Millones de pesos"
     )

hist(datos$precio_millon,
     main = "Frecuencia de Precio de Vivienda",
     xlab = "Millones de pesos",
     ylab = "Frecuencia"
     )

Análisis Exploratorio de Área Construida

summary(datos$Area_contruida)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    80.0    86.0    97.0   115.7   130.0   195.0
Desviacion_Estandar_Area <- sd(datos$Area_contruida)
Desviacion_Estandar_Area
## [1] 35.54332

Según estos datos podemos afirmar que la variable Área Construida (metros cuadrados), cuenta con valores entre 80.0 siendo su mínimo y 195.0 siendo su máximo, además cuenta con una media de 115.7 y una desviación estándar de 35.54332 metros cuadrados.

plot(datos$Area_contruida, 
     main = "Distribución de Área Construida",
     xlab = "Índice del dato",
     ylab = "Metros cuadrados"
     )

hist(datos$Area_contruida,
     main = "Frecuencia de Área Construida",
     xlab = "Metros cuadrados",
     ylab = "Frecuencia"
     )

Modelo \(\small y= \beta _{0} + \beta _{1}{x} +e\)

Análisis Exploratorio Bivariado (Y = Precio, X = area)

summary(datos)
##  Area_contruida  precio_millon  
##  Min.   : 80.0   Min.   :240.0  
##  1st Qu.: 86.0   1st Qu.:251.2  
##  Median : 97.0   Median :305.0  
##  Mean   :115.7   Mean   :332.1  
##  3rd Qu.:130.0   3rd Qu.:395.0  
##  Max.   :195.0   Max.   :480.0
cor(datos$Area_contruida,datos$precio_millon)
## [1] 0.9190295

Según estos datos, podemos afirmar que nuestro modelo puede usar los datos en (80.0, 195.0), y obtener una respuesta en (240.0, 480.0), intentar observar o predecir valores fuera de estos rangos sería una extrapolación para la que el modelo no aporta datos y por ende no es correcto. El Coeficiente de Correlación de Pearson nos arroja un valor de 0.9190295, este coeficiente cuenta con valores en el rango (-1, 1), mientras más positivo, las variables están más relacionadas de manera directamente proporcional, por ende, podemos concluir que al tener un Coeficiente de Correlación de 0.9190295, estas dos variables tienen una relación directamente proporcional, a más aumente una, más aumenta la otra.

plot(datos$Area_contruida,datos$precio_millon
     , main = "Precio de Vivienda vs. Área Construida",
     xlab = "Metros cuadrados",
     ylab = "Millones de pesos"
     )

En el diagrama se denota una relación creciente, pero con una aparente curvatura leve que podría indicar que el modelo no estaría correctamente adaptado para hacer uso de la Regresión Lineal Simple y podría necesitar una transformación, esto se corroborará con la validación de supuestos más adelante.

#x y
x <- Area_contruida
y <- precio_millon

Calculo del modelo de regresión lineal simple y coeficientes \(\small \beta _{0}\) y \(\small \beta _{1}\)

Nuestro modelo de regresión lineal simple del precio en función de una función de área y la suma del error es \(precio= 86.234 + 2.124\cdot area +e\).

v <- summary(datos)
mod=lm(y~x)
summary(mod)
## 
## Call:
## lm(formula = y ~ x)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -51.673 -25.612  -6.085  24.875  67.650 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   86.234     22.479   3.836 0.000796 ***
## x              2.124      0.186  11.422 3.45e-11 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 33.05 on 24 degrees of freedom
## Multiple R-squared:  0.8446, Adjusted R-squared:  0.8381 
## F-statistic: 130.5 on 1 and 24 DF,  p-value: 3.45e-11

Si bien \(\small \beta _{0}\) es significativa no tiene sentido interpretar un apartamento de 0 metros cuadrados, \(\small \beta _{1}\) por su parte también es significativa y podemos interpretarle como que por cada aumento unitario en el \(area\) el \(precio\) aumenta en promedio 2.124 millones, lo que quiere decir que su relación es positiva.

b <- summary(mod)
n <- b$df[2]+2
bs <- data.frame(b$coefficients)
b0 <- data.frame(b$coefficients)$Estimate[1]
b1 <- data.frame(b$coefficients)$Estimate[2]
b0d <- c(c(bs$Std..Error[1]),c(bs$t.value[1]),c(bs$Pr...t..[1]))
b1d <- c(c(bs$Std..Error[2]),c(bs$t.value[2]),c(bs$Pr...t..[2]))

sigb0 <- function(a){
  ta <- qt(a/2,n-2,lower.tail = F)
  if(abs(b0d[2])>ta)return('significativa')
  if(abs(b0d[2])<=ta)return('no significativa')
}
sigb1 <- function(a){
  ta <- qt(a/2,n-2,lower.tail = F)
  if(abs(b1d[2])>ta)return('significativa')
  if(abs(b1d[2])<=ta)return('no significativa')
}
#nos regresa si los valores $\beta _{0}$ y $\beta _{1}$  son significativos
sigb0(0.05)
## [1] "significativa"
sigb1(0.05)
## [1] "significativa"

Calculo de intervalo de confianza e interpretación para \(\small\beta _{1}\)

un intervalo de confianza al 95% para el coeficiente \(\small\beta _{1}\) es \(\small(1.740, 2.508)\) esto quiere decir que el verdadero \(\small\beta _{1}\) esta entre esos valores con un 95% de confianza, 0 no está dentro de este intervalo por eso podemos decir con confianza que el valor es diferente de 0 y por lo tanto significativa.

icb1 <- function(a){
  ta <- qt(a/2,n-2,lower.tail = F)
  ib1 <- c(b1-(ta*b1d[1]),b1+(ta*b1d[1]))
  return(ib1)
}
icb1(0.05)# devuelve un intervalo con confianza a para  b1
## [1] 1.740170 2.507771

Haciendo una prueba de significancia con \(\small t\) para nuestro \(\small\beta _{1}\) nuestra hipótesis inicial es \(\small\beta _{1}\) es igual a 0 y la alternativa es que \(\small\beta _{1}\) es diferente de 0 y tenemos que nuestro \(\small t _{0}\) es 11.42 y nuestro t es 2.06, ahora bien nuestra hipótesis inicial es rechazada, así que nuestra \(\small\beta _{1}\) es diferente de 0 y por lo tanto significativa.

b1d[2]
## [1] 11.4217
qt(0.05/2,n-2,lower.tail = F)
## [1] 2.063899
sigb1 <- function(a){
  ta <- qt(a/2,n-2,lower.tail = F)
  if(abs(b1d[2])>ta)return('significativa')
  if(abs(b1d[2])<=ta)return('no significativa')
}
sigb1(0.05)
## [1] "significativa"

Calculo e interpretación del indicador de bondad y ajuste \(\small R^{2}\)

el \(\small R^{2}\) es igual a 0.845 lo que significa que el modelo explica el 84.5% de la variabilidad del precio de las viviendas , y que modelo tiene un ajuste relativamente bueno a la recta de regresión, sin embargo no significa que cumpla con los supuestos de linealidad, aunque podría mejorar.

b$r.squared
## [1] 0.8446152

Precio promedio estimado para un apartamento de 110 metros cuadrados y consideraciones adicionales

Para poder hacer el calculo del precio promedio estimado para este apartamento primero tenemos que saber sí la \(x\) es interpolada, y ya que esta dentro del rango de las \(x\), si podemos proceder a hacer cálculos.

xo <- 110
c(xo>min(x),xo<max(x))
## [1] TRUE TRUE

El precio promedio estimado para un apartamento de 110 metros cuadrados esta en un intervalo de 306.3133 millones a 333.4279 millones, por lo tanto, un apartamento con un precio de 200 millones es una gran oferta en esta zona, una consideración adicional a tener en en cuenta es que la razón por la que tenga tan buen precio puede estar relacionada con el estado del apartamento, su poca deseabilidad en el mercado, o que hayan sucedido cosas terribles en el lugar.

ta <- qt(0.05/2,n-2,lower.tail = F)
yi <- predict(mod, newdata = data.frame(x=c(xo)), se.fit = T)
yo <- yi$fit
ys <- yi$se.fit
c(yo-(ta*ys),yo+(ta*ys))
##        1        1 
## 306.3133 333.4279

Validación de supuestos del modelo \(\small y= \beta _{0} + \beta _{1}\cdot x +e\)

Supuesto de Media Cero

mod = lm(datos$precio_millon~datos$Area_contruida)
summary(mod)
## 
## Call:
## lm(formula = datos$precio_millon ~ datos$Area_contruida)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -51.673 -25.612  -6.085  24.875  67.650 
## 
## Coefficients:
##                      Estimate Std. Error t value Pr(>|t|)    
## (Intercept)            86.234     22.479   3.836 0.000796 ***
## datos$Area_contruida    2.124      0.186  11.422 3.45e-11 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 33.05 on 24 degrees of freedom
## Multiple R-squared:  0.8446, Adjusted R-squared:  0.8381 
## F-statistic: 130.5 on 1 and 24 DF,  p-value: 3.45e-11
ei = mod$residuals
yi_mod = mod$fitted.values

round(mean(ei),3)
## [1] 0

Gracias a esto validamos que la media es efectivamente 0, se valida el supuesto de Media Cero.

Supuesto de Varianza Constante

plot(yi_mod,ei,main="Residuales vs. Valores Ajustados", xlab = "Valores Ajustados", ylab = "Residuales")
abline(h=0,lty=2,col=2)

require(ggplot2)
## Loading required package: ggplot2
## Warning: package 'ggplot2' was built under R version 4.0.4

data_val=data.frame(yi_mod,ei)

ggplot(data_val,aes(x=yi_mod,y=ei))+geom_point()+theme_bw()+geom_smooth()
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'

Aquí podemos observar que no se estaría cumpliendo gráficamente los requerimientos para validar una varianza constante, el gráfico sólo nos dice que sospechemos la no linealidad del modelo, no se valida el supuesto de Varianza constante.

Supuesto de Normalidad

ei
##          1          2          3          4          5          6          7 
## -18.895336  48.137608  32.649962 -51.672533 -28.895336  25.617018  32.691142 
##          8          9         10         11         12         13         14 
## -22.135041 -26.771366   2.691142 -31.019307  47.245540  21.104664 -26.771366 
##         15         16         17         18         19         20         21 
## -18.135041 -18.895336 -18.895336  32.649962  14.154080  -1.151513  67.649962 
##         22         23         24         25         26 
## -11.019307  22.649962  14.732752 -50.408120 -37.308858
hist(ei,col="gray", main = "Histograma de Residuales", xlab ="Residuales", ylab = "Frecuencia")

qqnorm(ei)
qqline(ei,col="red")

shapiro.test(ei)
## 
##  Shapiro-Wilk normality test
## 
## data:  ei
## W = 0.95489, p-value = 0.3009
alpha = 0.05

Se observa que el gráfico de normalidad presenta curvaturas y una tendencia a alejarse de la recta teórica, la prueba de Shapiro-Wilk nos arroja un valor P de 0.3009, mayor que el valor alfa de 0.05, entre la prueba gráfica y el test de Shapiro-Wilk prima la validación grafica, que en este caso no se puede otorgar por inconsistencia con los requerimientos de validación, no se valida el supuesto de Normalidad.

Supuesto de Independencia de Errores

Dado que todos estos registros no corresponden a datos en el tiempo no se tiene un orden para realizar la validación del supuesto de Independencia de variables. Se valida por definición del tipo de datos de corte transversal.

Supuesto de la linealidad

mod.lof<-lm(datos$precio_millon~as.factor(datos$Area_contruida)) 
anova(mod,mod.lof)
## Analysis of Variance Table
## 
## Model 1: datos$precio_millon ~ datos$Area_contruida
## Model 2: datos$precio_millon ~ as.factor(datos$Area_contruida)
##   Res.Df     RSS Df Sum of Sq      F  Pr(>F)  
## 1     24 26212.2                              
## 2     12  5363.4 12     20849 3.8872 0.01305 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Para un nivel de confianza del 95%, se rechaza la Hipótesis de que el modelo sea lineal, por ende, se le debe hacer alguna transformación para que valide supuestos y linealidad.

Vemos que nuestro modelo actual incumple varios supuestos asi que hemos decidido cambiar de modelo.

Validación de supuestos del modelo \(\small y= \beta _{0} + \beta _{1}\cdot\frac{1}{x} +e\)

area_transformada = 1 / datos$Area_contruida
modelo_transformado <- lm(datos$precio_millon~area_transformada)
summary(modelo_transformado)
## 
## Call:
## lm(formula = datos$precio_millon ~ area_transformada)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -36.987 -16.743  -5.023  18.547  44.379 
## 
## Coefficients:
##                    Estimate Std. Error t value Pr(>|t|)    
## (Intercept)          635.35      18.27   34.77  < 2e-16 ***
## area_transformada -32464.72    1895.32  -17.13 5.84e-15 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 23.05 on 24 degrees of freedom
## Multiple R-squared:  0.9244, Adjusted R-squared:  0.9212 
## F-statistic: 293.4 on 1 and 24 DF,  p-value: 5.839e-15

Supuesto de Media Cero

El supuesto de la media cero se cumple por defecto.

Supuesto de Varianza Constante

En el gráfico de residuales vs predichos podemos observar que se cumple el supuesto de varianza constante y que mejora con respecto al modelo orginal.

plot(fitted(modelo_transformado), residuals(modelo_transformado), main = "Residuales vs valores ajustados")
abline(h = 0, lty = 2,col= 2)

Supuesto de Normalidad

No hay evidencia fuerte en contra del supuesto de que los residuales distribuyen normal, ademas se observa que el gráfico de probabilidad normal mejora con respecto al módelo original.

plot(modelo_transformado, 2)

Dado que el valor p = 0.5871 del test shapiro es suficientemente grande confirmamos que se cumple la hipotesis de normalidad.

test_shapiro = shapiro.test(residuals(modelo_transformado))
test_shapiro
## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(modelo_transformado)
## W = 0.96859, p-value = 0.5871

Supuesto de Independencia de Errores

Dado que estos registros no corresponden a datos en el tiempo no se tiene un orden temporal para realizar la validación de este supuesto. Se valida por definición del tipo de datos de corte

Supuesto de la linealidad

De la prueba de falta de ajuste podemos observar que el valor F ajuste(1.38) es inferior a f k.2,n-k(0.37) al 95% por tanto no hay evidencia que rechace la linealidad del modelo.

mod.aux <- lm(datos$precio_millon~as.factor(area_transformada)) 
anova(modelo_transformado,mod.aux)
## Analysis of Variance Table
## 
## Model 1: datos$precio_millon ~ area_transformada
## Model 2: datos$precio_millon ~ as.factor(area_transformada)
##   Res.Df     RSS Df Sum of Sq      F Pr(>F)
## 1     24 12755.6                           
## 2     12  5363.4 12    7392.2 1.3783 0.2935
e <- anova(modelo_transformado,mod.aux)
k<- e$Df[2]+2
fb <- e$F[2]
rlt <- function(a){
  fa <- qf(a,k-2,n-k,lower.tail = F)
  if(fb>fa) return('no lineal')
  if(fb<fa) return('lineal')
}
rlt(0.05)
## [1] "lineal"

Comparación de ambos modelos

Decidimos entonces usar el modelo \(\small y= \beta _{0} + \beta _{1}\cdot\frac{1}{x} +e\) porqué dado que tiene un \(\small R^{2}\) más alto tiene un mejor ajuste en su recta de regresión y alto puede hacer predicciones más útiles, además de que cumple todos los supuestos de una regresión lineal.

Modelo \(\small y= \beta _{0} + \beta _{1}\cdot\frac{1}{x} +e\)

#x y
x <- 1/Area_contruida
y <- precio_millon

Veamos la correlación entre el precio en millones de pesos y el reciproco del area en metros cuadrados. Se observa que a mayor cantidad en el area transformada el precio de las viviendas disminuye y su relación es fuerte de acuerdo con el coeficiente de correlación de Pearson (-0,96) mejorando la relación con respecto modelo original con coeficiente de correlación de Pearson(0.92).

  area_transformada = 1 / datos$Area_contruida
  plot(datos$precio_millon, area_transformada, ylab = "Precio en Millones",
     xlab = "1/Area" , main = "Gráfico de dispersión Y = β0 + β1(-1/X) + ε")

  cor(datos$precio_millon, area_transformada)
## [1] -0.9614495

Calculo del modelo de regresión lineal simple y coeficientes \(\small\beta _{0}\) y \(\small\beta _{1}\)

Nuestro modelo de regresión lineal simple del precio en función de una función de área y la suma del error es \(precio= 635.3498 - 32464.72\cdot \frac{1}{area} +e\).

x <- 1/Area_contruida
y <- precio_millon
v <- summary(datos)
mod=lm(y~x)
summary(mod)
## 
## Call:
## lm(formula = y ~ x)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -36.987 -16.743  -5.023  18.547  44.379 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)    635.35      18.27   34.77  < 2e-16 ***
## x           -32464.72    1895.32  -17.13 5.84e-15 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 23.05 on 24 degrees of freedom
## Multiple R-squared:  0.9244, Adjusted R-squared:  0.9212 
## F-statistic: 293.4 on 1 and 24 DF,  p-value: 5.839e-15

Si bien \(\small\beta _{0}\) es significativa no tiene sentido interpretar un apartamento de infinitos metros cuadrados negativos, \(\small \beta _{1}\) por su parte también es significativa y podemos interpretarlo como que por cada aumento unitario en \(\small \frac{1}{área}\) el precio disminuye en promedio 32464.72 millones, lo que quiere decir que si el área aumenta el precio también lo hace de un modo reciproco.

b <- summary(mod)
n <- b$df[2]+2
bs <- data.frame(b$coefficients)
b0 <- data.frame(b$coefficients)$Estimate[1]
b1 <- data.frame(b$coefficients)$Estimate[2]
b0d <- c(c(bs$Std..Error[1]),c(bs$t.value[1]),c(bs$Pr...t..[1]))
b1d <- c(c(bs$Std..Error[2]),c(bs$t.value[2]),c(bs$Pr...t..[2]))

sigb0 <- function(a){
  ta <- qt(a/2,n-2,lower.tail = F)
  if(abs(b0d[2])>ta)return('significativa')
  if(abs(b0d[2])<=ta)return('no significativa')
}
sigb1 <- function(a){
  ta <- qt(a/2,n-2,lower.tail = F)
  if(abs(b1d[2])>ta)return('significativa')
  if(abs(b1d[2])<=ta)return('no significativa')
}
#nos regresa si los valores $\beta _{0}$ y $\beta _{1}$  son significativos
sigb0(0.05)
## [1] "significativa"
sigb1(0.05)
## [1] "significativa"

Calculo de intervalo de confianza e interpretación para \(\small\beta _{1}\)

un intervalo de confianza al 95% para el coeficiente \(\small\beta _{1}\) es \(\small (-36376.47, -28552.97)\) esto quiere decir que el verdadero \(\small\beta _{1}\) esta entre esos valores con un 95% de confianza, 0 no esta dentro de este intervalo por eso podemos decir con confianza que el valor es diferente de 0 y por lo tanto significativa.

icb1 <- function(a){
  ta <- qt(a/2,n-2,lower.tail = F)
  ib1 <- c(b1-(ta*b1d[1]),b1+(ta*b1d[1]))
  return(ib1)
}
icb1(0.05)# devuelve un intervalo con confianza a para  b1
## [1] -36376.47 -28552.97

Haciendo una prueba de significancia con t para nuestro \(\small\beta _{1}\) nuestra hipótesis inicial es \(\small \beta _{1}\) es igual a 0 y la alternativa es que \(\small \beta _{1}\) es diferente de 0 y tenemos que nuestro \(\small t _{0}\) es -17.13 y nuestro t es 2.06, ahora bien nuestra hipótesis inicial es rechazada, así que nuestra \(\small \beta _{1}\) es diferente de 0 y por lo tanto significativa.

b1d[2]
## [1] -17.12887
qt(0.05/2,n-2,lower.tail = F)
## [1] 2.063899
sigb1 <- function(a){
  ta <- qt(a/2,n-2,lower.tail = F)
  if(abs(b1d[2])>ta)return('significativa')
  if(abs(b1d[2])<=ta)return('no significativa')
}
sigb1(0.05)
## [1] "significativa"

Calculo e interpretación del indicador de bondad y ajuste \(\small R^{2}\)

el \(\small R^{2}\) es igual a 0.924 lo que significa que el modelo explica el 92.44% de la variabilidad del precio de las viviendas, además de que el modelo tiene un muy buen ajuste a su recta de regreción y alto puede hacer predicciones útiles.

b$r.squared
## [1] 0.9243852

Precio promedio estimado para un apartamento de 110 metros cuadrados y consideraciones adicionales

Para poder hacer el cálculo del precio promedio estimado para este apartamento primero tenemos que saber sí la \(x\) es interpolada, y ya que está dentro del rango de las \(x\), si podemos proceder a hacer cálculos, no hay que olvidar que nuestro modelo cambio y por lo tanto la \(x\) también.

xo <- 1/110
c(xo>min(x),xo<max(x))
## [1] TRUE TRUE

El precio promedio estimado para un apartamento de 110 metros cuadrados esta en un intervalo de 330.8332 millones a 349.5988 millones, por lo tanto un apartamento con un precio de 200 millones es una gran oferta en esta zona, una consideración adicional a tener en en cuenta es que la razón por la que tenga tan buen precio puede estar relacionada con el estado del apartamento, su poca deseabilidad en el mercado, o que hayan sucedido cosas terribles en el lugar.

ta <- qt(0.05/2,n-2,lower.tail = F)
yi <- predict(mod, newdata = data.frame(x=c(xo)), se.fit = T)
yo <- yi$fit
ys <- yi$se.fit
c(yo-(ta*ys),yo+(ta*ys))
##        1        1 
## 330.8332 349.5988

Función para calcular el intervalo de confianza de un \(\small \beta _{1}\) dados X y Y

Una función para calcular el intervalo de confianza de un \(\small \beta _{1}\) en un \(\small 100(1-a)%\) dados dos vectores X y Y es:

icb1 <- function(x,y,a){
  mod=lm(y~x)
  b <- summary(mod)
  n <- b$df[2]+2
  b1 <- data.frame(b$coefficients)$Estimate[2]
  b1d <- c(c(bs$Std..Error[2]),c(bs$t.value[2]),c(bs$Pr...t..[2]))
  ta <- qt(a/2,n-2,lower.tail = F)
  ib1 <- c(b1-(ta*b1d[1]),b1+(ta*b1d[1]))
  return(ib1)
}
icb1(x,y,0.05)
## [1] -36376.47 -28552.97