1. Predicción de los precios de las acciones. Analizar el comportamiento de los precios de las Acciones de Ecopetrol según la variación del precio del barril de petróleo WTI producido en Colombia. Se tienen los siguientes precios

prec_ecop=c(1090, 1170, 1160, 1230, 1155, 1165, 1205, 1170, 1150, 1130, 1110, 1105, 
1085, 1060, 1035, 1015, 955, 961) 
  
prec_petr=c(35.62,36.31,37.35,34.95,34.53,35.81,36.14,37.50,37.80,36.81,37.87,37.04,
36.76,35.97,33.97,33.27,31.41,30.44)

datos_ecopetrol=data.frame(prec_ecop,prec_petr)

head(datos_ecopetrol)
prec_ecop prec_petr
1090 35.62
1170 36.31
1160 37.35
1230 34.95
1155 34.53
1165 35.81

a.Proponga un modelo de regresión lineal simple que permita predecir el valor de las Acciones de Ecopetrol con base en el Precio del barril de petróleo en Colombia. Indique la ecuación de regresión y el valor del R2.

Por lo cual la ecuación es la siguiente:

\[\begin{equation*} \widehat{E[y|_{x_i}]} = \widehat{y_{i}}=\widehat{\beta_{0}}+\widehat{\beta_{1}}x_{i} \end{equation*}\]

\[Y(PrecioAccionEcopetrol) = 177.768 + 26.192*X(PrecioBarrilPetroleo)\] \[ R^2 = 0.4692 \]

modelo_ecopetrol=lm(prec_ecop~prec_petr)

summary(modelo_ecopetrol)
## 
## Call:
## lm(formula = prec_ecop ~ prec_petr)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -59.90 -40.74 -15.94  33.40 136.82 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)   
## (Intercept)  177.768    232.828   0.764  0.45627   
## prec_petr     26.192      6.542   4.004  0.00102 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 57.13 on 16 degrees of freedom
## Multiple R-squared:  0.5005, Adjusted R-squared:  0.4692 
## F-statistic: 16.03 on 1 and 16 DF,  p-value: 0.001024

b.Pruebe la significancia del modelo propuesto en “a)” plantee las hipótesis respectivas y use el concepto de Valor _p para tomar la decisión sobre las hipótesis. Use α = 0.05

##Significancia de Intercepto b0

\[H_{0} : \beta_{0} = 0\] \[H_{1} : \beta_{0} \neq 0 \] Luego en este caso el valor p para el b0 es de 0.45627, por lo tanto no se rechaza la hipotesis H0 ya que el valor p es mayor que el alfa de 0.05

##Significancia de Pendiente b1

\[H_{0} : \beta_{1} = 0\] \[H_{1} : \beta_{1} \neq 0 \] Para la pendiente el valor p fue de 0.00102, por lo que se rechaza la H0 al ser nuestro valor p menor que el alfa de 0.05

require(ggplot2)
require(plotly)

grafica_modelo_ecopetrol=ggplot(data = datos_ecopetrol,aes(y=prec_ecop, x=prec_petr)) + geom_point() + geom_smooth(method = "loess",formula = "y ~ x")

ggplotly(grafica_modelo_ecopetrol)

c.Interprete los coeficientes del modelo propuesto en “a)”

Interpretación: Para empezar el b0 nos dice que si no existe precio del barril de petroleo o este es cero, entonces el precio de una accion de ecopetrol seria igual a $177.768 aproximadamente. Por otro lado, Lo que indica la pendiente b1, es cual es el aumento del precio de la acción por cada incremento en un peso del barril de petroleo, es decir, que si el precio del barril aumenta un 1% entonces, en el precio de las acciones de ecopetrol el efecto sera que los precios aumenten 26.192 millones aproximadamente por cada peso adicional. Por otro lado, revisando el valor p nos dice que el coeficiente b1 si es significativo ya que nos presenta un valor de cero. Además, el r2 nos dice que se explica en un 46.92%% los precios de cada acción.

d.Haga un análisis de los residuos. ¿Qué supuesto no se cumple?

Supuestos:

1.Errores del modelo tienen media cero en el modelo. Se cumple.

2.Errores del modelo tienen varianza Constante. Se cumple.

3.Errores del modelo se distribuyen normal. No se cumple.

4.Errores del modelo son independientes. No se cumple.

### 1.Para el supuesto de media igual a cero para el modelo, podemos ver que el supuesto si se cumple ya que el valor p es 1 y es mayor que el alfa de 0.05 diciendonos que no se rechaza la H0.


##h0= Media igual a 0

##h1= Media diferente a 0


summary(modelo_ecopetrol$residuals)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##  -59.90  -40.74  -15.94    0.00   33.40  136.82
t.test(modelo_ecopetrol$residuals, mu = 0)
## 
##  One Sample t-test
## 
## data:  modelo_ecopetrol$residuals
## t = -4.2309e-16, df = 17, p-value = 1
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  -27.56364  27.56364
## sample estimates:
##     mean of x 
## -5.527407e-15
### 2. Para el supuesto de varianza constante, podemos ver que tenemos un valor p de 0.9813, el cual es mayor que el alfa de 0.05 mostrandonos que no se rechaza la h0, por lo tanto el supuesto se cumple.

##h0  : no existe heteroscedasticidad (homoscedasticidad).
##h1 : existe heteriscedastucidad.

library(lmtest)

lmtest::gqtest(modelo_ecopetrol)
## 
##  Goldfeld-Quandt test
## 
## data:  modelo_ecopetrol
## GQ = 0.17924, df1 = 7, df2 = 7, p-value = 0.9813
## alternative hypothesis: variance increases from segment 1 to 2
### 3. Para el supuesto de normalidad se realizo la prueba de Shapiro sobre los residuales del modelo y se logro ver que el valor p es 0.04276, el cual es menor que el alfa de 0.05. Esto, nos indica que se rechaza la H0 mostrando que el supuesto de normalidad no se cumple.


shapiro.test(modelo_ecopetrol$residuals)
## 
##  Shapiro-Wilk normality test
## 
## data:  modelo_ecopetrol$residuals
## W = 0.89259, p-value = 0.04276
### 4. Para el supuesto de independencia podemos ver que el valor p fue de 0.0004666 menor al valor del alfa de 0.05, por lo que rechazamos la h0, lo cual nos dice que el supuesto no se cumple 

##H0 igual a 0. Errores independientes. 
##H1 diferente 0 Errores no independientes.

lmtest::dwtest(modelo_ecopetrol)
## 
##  Durbin-Watson test
## 
## data:  modelo_ecopetrol
## DW = 0.74504, p-value = 0.0004666
## alternative hypothesis: true autocorrelation is greater than 0

e.Concluya sobre la validez del modelo propuesto en a)

Respuesta: El modelo no es valido ya que dos de los cuatro supuestos no se cumplen. Ademas, el R2 nos da un valor de 46,92% el cual no explica de una manera significativa la relacion entre el precio de la accion de ecopetrol y el precio por barril.

3. Con base en los datos de precios de vivienda de la actividad en clase realizar un informe que contenga los siguientes puntos utilizando R y RMarkdown (publicar el informe final en Rpubs presentando código, resultados e interpretaciones).

a.Realice un filtro a la base de datos e incluya solo las ofertas de apartamentos, de la zona norte de la ciudad con precios inferiores a los 500 millones de pesos y áreas menores a 300 mt2. Presente los primeros 3 registros de la base y algunas tablas que comprueben la consulta. (Adicional un mapa con los puntos de la base, discutir si todos los puntos se ubican en la zona norte o se presentan valores en otras zonas, por que?).

Respuesta: No todos los puntos se encuentran en la zona norte y esto se puede dar si existe algun error en la base de datos en la cual se haya puesto un barrio de una zona diferente como Zona Norte o que existan datos nulos.

library(readxl)
datos_vivienda = read_excel("C:/Users/Luisa/Downloads/Datos_Vivienda (1).xlsx")


ID=1:dim(datos_vivienda)[1]

datos_vivienda=data.frame(ID,datos_vivienda)

pos=which(datos_vivienda$Tipo=="Apartamento" & datos_vivienda$Zona=="Zona Norte" & datos_vivienda$precio_millon < 500 & datos_vivienda$Area_contruida < 300)


datos_sub=datos_vivienda[pos,]
head(datos_sub,3)
ID Zona piso Estrato precio_millon Area_contruida parqueaderos Banos Habitaciones Tipo Barrio cordenada_longitud Cordenada_latitud
31 31 Zona Norte 2 3 135 56 1 1 3 Apartamento torres de comfandi -76.46745 3.40763
71 71 Zona Norte NA 3 78 54 2 1 3 Apartamento chiminangos -76.47820 3.44898
89 89 Zona Norte NA 5 340 106 2 2 3 Apartamento la flora -76.48200 3.43500
sapply(datos_sub, class)
##                 ID               Zona               piso            Estrato 
##          "integer"        "character"        "character"          "numeric" 
##      precio_millon     Area_contruida       parqueaderos              Banos 
##          "numeric"          "numeric"        "character"          "numeric" 
##       Habitaciones               Tipo             Barrio cordenada_longitud 
##          "numeric"        "character"        "character"          "numeric" 
##  Cordenada_latitud 
##          "numeric"
# column to numeric
datos_sub <- transform(datos_sub, parqueaderos = as.numeric(parqueaderos))

#Reemplazar NA
datos_sub[is.na(datos_sub)] <- 0

require(leaflet)

leaflet() %>% addCircleMarkers(lng = datos_sub$cordenada_longitud,lat = datos_sub$Cordenada_latitud,radius = 0.3,color = "black",label = datos_sub$ID) %>% addTiles()

b.Realice un análisis exploratorio de datos enfocado en la correlación entre la variable respuesta (precio del apartamento) en función del área construida, estrato y si tiene parqueadero. Use gráficos interactivos con plotly e interprete los resultados.

Interpretación: Segun lo visto en las graficas podemos ver que la relacion con mas sentido es la que existe entre el precio de la vivienda y el area construida de la misma. Por otro lado, se logro ver de que la relacion de los estratos y parqueaderos son un poco diferentes y estas no se presentande manera lineal. Es importante recalcar, que en la vida real si una vivienda se encuentra en un lugar de estrato alto, con mayor cantidad de parqueaderos y con mayor cantidad de area construida se puede en cierto modo afirmar que la vivienda sera mucho mas costosa.

require(table1)
table1(~Area_contruida+precio_millon+Estrato+parqueaderos,data = datos_sub)
Overall
(N=1077)
Area_contruida
Mean (SD) 85.9 (34.0)
Median [Min, Max] 76.0 [35.0, 287]
precio_millon
Mean (SD) 234 (110)
Median [Min, Max] 220 [65.0, 495]
Estrato
Mean (SD) 4.19 (0.934)
Median [Min, Max] 4.00 [3.00, 6.00]
parqueaderos
Mean (SD) 0.889 (0.698)
Median [Min, Max] 1.00 [0, 4.00]
#Diferentes correlaciones entre las variables

library(GGally)
library(dplyr)
ggpairs(select_if(datos_sub, is.numeric), lower = list(continuous = "smooth"),diag = list(continuous = "barDiag"), axisLabels = "none")

require(ggplot2)
require(plotly)

g1=ggplot(data = datos_sub,aes(y=precio_millon,x=Area_contruida)) + geom_point() + geom_smooth()

ggplotly(g1)
correlacion_vivienda = cor(datos_sub$Area_contruida, datos_sub$precio_millon)
correlacion_vivienda
## [1] 0.6937608
#Interpretación: Podemos ver que la correlacion nos dice que existen una relacion positiva entre el precio  y el area construida, es decir que, mientras el area construida sea mayor el precio sera mas alto.
require(ggplot2)
require(plotly)

g2=ggplot(data = datos_sub,aes(y=precio_millon,x=datos_sub$Estrato)) + geom_point() + geom_smooth()

ggplotly(g2)
## Warning: Use of `datos_sub$Estrato` is discouraged. Use `Estrato` instead.
## Use of `datos_sub$Estrato` is discouraged. Use `Estrato` instead.
## `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'
## Warning: Computation failed in `stat_smooth()`:
## x has insufficient unique values to support 10 knots: reduce k.
correlacion_vivienda_e = cor(datos_sub$Estrato, datos_sub$precio_millon)
correlacion_vivienda_e
## [1] 0.8230851
#Interpretación: para empezar la grafica no nos muestra relacion lineal, pero al realizar la funcion de correlacion nos dice que el estrato y el precio de una vivienda si estan correlacionadas ya que si el estrato sube, entonces el precio de la vivienda sera mayor.
require(ggplot2)
require(plotly)

g3=ggplot(data = datos_sub,aes(y=precio_millon,x=datos_sub$parqueaderos)) + geom_point() + geom_smooth()

ggplotly(g3)
## Warning: Use of `datos_sub$parqueaderos` is discouraged. Use `parqueaderos` instead.
## Use of `datos_sub$parqueaderos` is discouraged. Use `parqueaderos` instead.
## `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'
## Warning: Computation failed in `stat_smooth()`:
## x has insufficient unique values to support 10 knots: reduce k.
correlacion_vivienda_p = cor(datos_sub$parqueaderos, datos_sub$precio_millon)
correlacion_vivienda_p
## [1] 0.5472535
#Interpretación: para empezar la grafica no nos muestra relacion lineal, pero al realizar la funcion de correlacion nos dice que la cantidad de parqueaderos y el precio de una vivienda si estan correlacionadas, pero no de manera muy fuerte ya que el resultado es de 0.54.
library(plotly)

plot_ly(x=datos_sub$Area_contruida,
        y=datos_sub$precio_millon,
        z=datos_sub$Estrato,
        size=.5)
#Interpretacion: En cuanto a este grafico, podemos ver que ciertas relaciones si tienen sentido, ya que se muestra que si el estrato es mayor y tiene mayor area contruida, la vivienda sera mucho mas costosa.
  1. Estime un modelo de regresión lineal múltiple con las variables del punto anterior e interprete los coeficientes si son estadísticamente significativos. Las interpretaciones deber están contextualizadas y discutir si los resultados son lógicos. Adicionalmente interprete el coeficiente R2 y discuta el ajuste del modelo e implicaciones (que podrían hacer para mejorarlo).

Interpretación: Como se encontro en el punto anterior y por logica se puede decir que si son logicos los resultados que se mostraron ya que si se explican las diferentes variables en el modelo y esto se ve con el r2 que fue de 76%. Ademas, de como se decia en el punto anterior la logica de usar las variables tiene sentido ya que eso es lo que se espera de una vivienda.Se podria realizar un ajuste al modelo para mostrar mejores datos, pero creo que lo que nos esta resultando tiene sentido, ya que tambien es importante usar la experiencia para la toma de decisiones.

#Significancia: Coeficientes significativos

  1. Lo que indica la pendiente b1, es cual es el aumento del precio por cada incremento en una unidad del area construida, es decir, que si el el area construida aumenta un 1% entonces, en las viviendas el efecto sera que los precios aumenten 0.95 millones aproximadamente por cada metro cuadrado adicional. Por otro lado, revisando el valor p nos dice que el coeficiente b1 si es significativo ya que nos presenta un valor de cero. Además, el r2 nos dice que se explica en un 76% aproximandamente los precios de una vivienda.

  2. Lo que indica la pendiente b2, es cual es el aumento del precio por cada incremento en el Estrato, es decir, que si el el estrato sube entonces, en las viviendas el efecto sera que los precios aumenten 68.99 millones aproximadamente. Por otro lado, revisando el valor p nos dice que el coeficiente b1 si es significativo ya que nos presenta un valor de cero. Además, el r2 nos dice que se explica en un 76% aproximandamente los precios de una vivienda.

  3. Lo que indica la pendiente b3, es cual es el aumento del precio por cada incremento en la cantidad de parqueaderos, es decir, que si los parqueaderos aumentan en numero entonces, en las viviendas el efecto sera que los precios aumenten 22.65 millones aproximadamente por cada parqueadero adicional. Por otro lado, revisando el valor p nos dice que el coeficiente b1 si es significativo ya que nos presenta un valor de cero. Además, el r2 nos dice que se explica en un 76% aproximandamente los precios de una vivienda.

\[Y(PrecioMillonVivienda) = -157.38220 + 0.94938*X1(AreaConstruida) + 68.99436*X2(Estrato) + 22.64906*X3(Parqueaderos) \] \[ R^2 = 0.762 \]

modelo_vivienda=lm(precio_millon~Area_contruida + Estrato + parqueaderos, data = datos_sub)

summary(modelo_vivienda)
## 
## Call:
## lm(formula = precio_millon ~ Area_contruida + Estrato + parqueaderos, 
##     data = datos_sub)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -216.571  -31.564   -1.213   27.889  224.053 
## 
## Coefficients:
##                  Estimate Std. Error t value Pr(>|t|)    
## (Intercept)    -157.38220    7.71190 -20.408  < 2e-16 ***
## Area_contruida    0.94938    0.06054  15.682  < 2e-16 ***
## Estrato          68.99436    2.26623  30.445  < 2e-16 ***
## parqueaderos     22.64906    2.73064   8.294 3.24e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 53.75 on 1073 degrees of freedom
## Multiple R-squared:  0.7629, Adjusted R-squared:  0.7623 
## F-statistic:  1151 on 3 and 1073 DF,  p-value: < 2.2e-16
  1. Realice la validación de supuestos del modelo e interprete los resultados (no es necesario corregir en caso de presentar problemas solo realizar sugerencias de que se podría hacer).
modelo_vivienda
## 
## Call:
## lm(formula = precio_millon ~ Area_contruida + Estrato + parqueaderos, 
##     data = datos_sub)
## 
## Coefficients:
##    (Intercept)  Area_contruida         Estrato    parqueaderos  
##      -157.3822          0.9494         68.9944         22.6491
par(mfrow=c(2,2))
plot(modelo_vivienda)

# 1.    Gráfico de ajuste: Me muestra el modelo ajustado. Que son los residuales contra valores ajustados. SE ESPERA QUE ESTO NO TENGA NINGUN COMPORTAMIENTO, ES DECIR, QUE SEA ALEATORIO. En este caso lo que se ve en la grafica tiene sentido ya que se muestran datos aleatorios. Es decir, que la asociacion lineal que suponemos si es posiblemente lineal.


# 2.    Gráfico de normalidad: (Normal Q-Q). Se muestra que no todos los datos están sobre la línea de la distribución normal, solo algunos se ven un poco alejados de la linea, pero estos pueden ajustarse lo suficiente para que funcione.

# 3 y 4. Son gráficos de datos atípicos. 

# Por otro lado, se debe transformar el modelo de regresion lineal, ya que  no todos los supuestos se cumplen. 
  1. Con el modelo identificado predecir el precio de un apartamento con 100 mt2, de estrato 4 y con parqueadero. ¿Si este apartamento lo están ofreciendo en 450 millones cual seria su opinión con base en el resultado del modelo considera que es una buena oferta?

Interpretacion: No considero que es una buena oferta ya que lo que se deberia pagar por un apartamento de estos es aproximadamente 236 Millones. Por otro lado, lo unico que no tuve encuenta en el modelo es que si esta oferta contiene mas de 1 parqueadero, este puede costar un poco mas, pero no llegaria a los 450 millones probablemente, ya que al reemplazarlo en el modelo el valor maximo para un apartamento con estas caracteristicas y 4 parqueaderos seria de 304 millones.

#Predicciones de modelo

predict(modelo_vivienda,newdata = list(Area_contruida=100,Estrato=4,parqueaderos=1),interval = "confidence")
##        fit      lwr      upr
## 1 236.1822 232.2449 240.1195
predict(modelo_vivienda,newdata = list(Area_contruida=100,Estrato=4,parqueaderos=2),interval = "confidence")
##        fit      lwr      upr
## 1 258.8313 251.7806 265.8819
predict(modelo_vivienda,newdata = list(Area_contruida=100,Estrato=4,parqueaderos=3),interval = "confidence")
##        fit      lwr      upr
## 1 281.4803 269.5918 293.3689
predict(modelo_vivienda,newdata = list(Area_contruida=100,Estrato=4,parqueaderos=4),interval = "confidence")
##        fit      lwr      upr
## 1 304.1294 287.0889 321.1699

f.Con las predicciones del modelo sugiera potenciales ofertas para una persona interesada en un apartamento en la zona norte con mas de 100 mt2 de área, de estrato 4, que tenga parqueadero y tenga encuentra que la persona tiene un crédito preaprobado de máximo 400 millones de pesos. Realice un análisis y presente en un mapa al menos 5 ofertas potenciales que debe discutir.

library(readxl)
datos_vivienda = read_excel("C:/Users/Luisa/Downloads/Datos_Vivienda (1).xlsx")


ID=1:dim(datos_vivienda)[1]

datos_vivienda=data.frame(ID,datos_vivienda)

pos_cliente=which(datos_vivienda$Tipo=="Apartamento" & datos_vivienda$Zona=="Zona Norte" & datos_vivienda$Barrio=="la flora" & datos_vivienda$precio_millon < 400 & datos_vivienda$precio_millon>236 & datos_vivienda$Area_contruida > 100 & datos_vivienda$Estrato==4 & datos_vivienda$parqueaderos>=1)


datos_sub_cliente=datos_vivienda[pos_cliente,]
head(datos_sub_cliente,)
ID Zona piso Estrato precio_millon Area_contruida parqueaderos Banos Habitaciones Tipo Barrio cordenada_longitud Cordenada_latitud
1264 1264 Zona Norte 4 4 380 123 1 3 3 Apartamento la flora -76.51437 3.48618
2606 2606 Zona Norte NA 4 350 130 1 2 3 Apartamento la flora -76.52100 3.49000
2632 2632 Zona Norte 1 4 290 108 1 2 3 Apartamento la flora -76.52115 3.48930
3067 3067 Zona Norte NA 4 265 125 2 3 4 Apartamento la flora -76.52353 3.48157
3189 3189 Zona Norte 2 4 380 126 2 3 4 Apartamento la flora -76.52432 3.48254
require(leaflet)

leaflet() %>% addCircleMarkers(lng = datos_sub_cliente$cordenada_longitud,lat = datos_sub_cliente$Cordenada_latitud,radius = 18,color = "black",label = datos_sub_cliente$ID) %>% addTiles()

4. Con base en los datos de arboles proponga un modelo de regresión lineal múltiple que permita predecir el peso del árbol en función de las covariables que considere importantes y seleccionándolas de acuerdo con un proceso adecuado. Tenga en cuenta realizar una evaluación de la significancia de los parámetros, interpretación y proponga un método de evaluación por medio de validación cruzada. Presente métricas apropiadas como el RMSE y MAE.

Interpretación: \[Y(PrecioMillonVivienda) = -9.1493 + 4.4257*X1(Diametro) + 0.5878*X2(Altura) \]

  1. Correlación: se encontro que la correlación de las covariables diametro y altura tienen una relacion directa fuerte en cuanto al peso.

  2. Significancia:

  1. R2: En este caso el resultado nos muestra que el diametro y la altura explican el peso en un 82% aproximadamente. \[ R^2 = 0.82 \]

  2. Supuestos: Se cumplen 3 supuestos de 4, el unico que no se cumple es el de normalidad. Lo que nos dice que si se le hacen ciertos ajustes al modelo, este realmente puede ser bueno para la prediccion.

##Base de Datos

library(readxl)
data_arboles <- read_excel("C:/Users/Luisa/Downloads/data arboles.xlsx", 
    col_types = c("text", "text", "numeric", 
        "numeric", "numeric"))

head(data_arboles,3)
finca mg peso diametro altura
FINCA_1 GENOTIPO_1 13.73 4.7 5.0
FINCA_1 GENOTIPO_1 14.58 5.3 5.6
FINCA_1 GENOTIPO_1 15.88 4.8 5.8
#Diferentes correlaciones entre las variables

library(GGally)
library(dplyr)
ggpairs(select_if(data_arboles, is.numeric), lower = list(continuous = "smooth"),diag = list(continuous = "barDiag"), axisLabels = "none")

### Validacion Cruzada

## Conjunto de 90 datos (80%(72) para modelar y 20%(18) para validar):

id_modelacion= sample(1:90, size=72)
id_modelacion
##  [1] 67 25 58 30 27 10 77 32 68 54 45 65 69 83 21 51 40 35 33  1 47 38 84 88 20
## [26] 75 14  3 12  6  2 19 63 76 89 85 18 46 11 43 79 22 62  8 82 41 61 66 23  7
## [51] 71 81  4 26 64 42 73 39  5 50 52 60 59 17 72 36 31  9 16 78 57 86
# data set con 72 datos para modelar:
modelacion_arboles= data_arboles[id_modelacion,]
head(modelacion_arboles)
finca mg peso diametro altura
FINCA_3 GENOTIPO_2 23.14 6.0 7.5
FINCA_1 GENOTIPO_1 9.93 3.5 4.3
FINCA_2 GENOTIPO_1 22.75 6.6 6.1
FINCA_1 GENOTIPO_2 8.73 4.0 5.3
FINCA_1 GENOTIPO_2 9.89 4.3 6.3
FINCA_1 GENOTIPO_2 16.62 5.9 7.1
# data set con 18 datos para validar:
validacion_arboles= data_arboles[-id_modelacion,]
head(validacion_arboles)
finca mg peso diametro altura
FINCA_1 GENOTIPO_1 7.47 2.2 3.5
FINCA_1 GENOTIPO_1 11.04 3.7 4.6
FINCA_1 GENOTIPO_1 13.61 4.6 5.5
FINCA_1 GENOTIPO_2 6.98 3.0 4.8
FINCA_1 GENOTIPO_2 14.44 5.5 7.4
FINCA_2 GENOTIPO_2 30.83 7.9 10.9
##Modelo Multiple

modelo_arboles=lm(peso~diametro + altura, data = modelacion_arboles)

summary(modelo_arboles)
## 
## Call:
## lm(formula = peso ~ diametro + altura, data = modelacion_arboles)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -5.3566 -2.3800  0.2778  1.9640  8.6928 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -7.29745    1.49266  -4.889 6.35e-06 ***
## diametro     4.60395    0.69496   6.625 6.32e-09 ***
## altura       0.08602    0.55997   0.154    0.878    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.969 on 69 degrees of freedom
## Multiple R-squared:  0.8201, Adjusted R-squared:  0.8149 
## F-statistic: 157.3 on 2 and 69 DF,  p-value: < 2.2e-16
### 1.Para el supuesto de media igual a cero para el modelo, podemos ver que el supuesto si se cumple ya que el valor p es 1 y es mayor que el alfa de 0.05 diciendonos que no se rechaza la H0.


##h0= Media igual a 0

##h1= Media diferente a 0


summary(modelo_arboles$residuals)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
## -5.3566 -2.3800  0.2778  0.0000  1.9640  8.6928
t.test(modelo_arboles$residuals, mu = 0)
## 
##  One Sample t-test
## 
## data:  modelo_arboles$residuals
## t = 5.4978e-16, df = 71, p-value = 1
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  -0.6877475  0.6877475
## sample estimates:
##  mean of x 
## 1.8963e-16
### 2. Para el supuesto de varianza constante, podemos ver que tenemos un valor p de 0.3079, el cual es mayor que el alfa de 0.05 mostrandonos que no se rechaza la h0, por lo tanto el supuesto se cumple.

##h0  : no existe heteroscedasticidad (homoscedasticidad).
##h1 : existe heteriscedastucidad.

library(lmtest)

lmtest::gqtest(modelo_arboles)
## 
##  Goldfeld-Quandt test
## 
## data:  modelo_arboles
## GQ = 1.3055, df1 = 33, df2 = 33, p-value = 0.224
## alternative hypothesis: variance increases from segment 1 to 2
### 3. Para el supuesto de normalidad se realizo la prueba de Shapiro sobre los residuales del modelo y se logro ver que el valor p es 0.02168, el cual es menor que el alfa de 0.05. Esto, nos indica que se rechaza la H0 mostrando que el supuesto de normalidad no se cumple.


shapiro.test(modelo_arboles$residuals)
## 
##  Shapiro-Wilk normality test
## 
## data:  modelo_arboles$residuals
## W = 0.96855, p-value = 0.0683
### 4. Para el supuesto de independencia podemos ver que el valor p fue de 0.416 menor al valor del alfa de 0.05, por lo que no rechazamos la h0, lo cual nos dice que el supuesto si se cumple 

##H0 igual a 0. Errores independientes. 
##H1 diferente 0 Errores no independientes.

lmtest::dwtest(modelo_arboles)
## 
##  Durbin-Watson test
## 
## data:  modelo_arboles
## DW = 1.8916, p-value = 0.3308
## alternative hypothesis: true autocorrelation is greater than 0
#Prediccion 

prediccion_modelo = predict(modelo_arboles, list(diametro=validacion_arboles$diametro, altura=validacion_arboles$altura))
prediccion_modelo
##         1         2         3         4         5         6         7         8 
##  3.132318 10.132864 14.353837  6.927306 18.660835 30.011390 18.097211 29.473574 
##         9        10        11        12        13        14        15        16 
## 23.230372 22.804387 33.582716 34.189351 15.274627 17.202229 24.280197 13.355627 
##        17        18 
## 11.070858 13.859033
## Validacion Peso Real

validacion_peso= validacion_arboles$peso


error_peso= validacion_peso - prediccion_modelo

resultado_peso=data.frame(validacion_peso, prediccion_modelo, error_peso)
resultado_peso
validacion_peso prediccion_modelo error_peso
7.47 3.132318 4.3376822
11.04 10.132864 0.9071361
13.61 14.353837 -0.7438373
6.98 6.927306 0.0526939
14.44 18.660835 -4.2208346
30.83 30.011390 0.8186101
17.01 18.097211 -1.0872113
32.69 29.473574 3.2164262
22.01 23.230372 -1.2203716
20.24 22.804386 -2.5643865
45.41 33.582716 11.8272839
47.87 34.189351 13.6806486
12.95 15.274627 -2.3246266
17.26 17.202229 0.0577709
31.16 24.280197 6.8798033
13.23 13.355626 -0.1256265
14.91 11.070858 3.8391420
18.49 13.859033 4.6309669
## MAE

## El error entre la prediccion y el valor real del peso es de 2.2557 y este modelo se equivoca en un 12% aproximadamente segun este indicador.

MAE= mean(abs(error_peso))
MAE
## [1] 3.47417
(MAE/mean(validacion_peso))*100
## [1] 16.56119
## RMSE

## Debido a que tenemos un RMSE pequeno podemos decir que el modelo es bueno.

RMSE = sqrt(mean((error_peso)^2))
RMSE
## [1] 5.132993
library(Metrics)
rmse(validacion_peso, prediccion_modelo)
## [1] 5.132993