En la tabla se observa que el promedio de precios de los autos mazda2 en la pagina olx Colombia es de $ 52.000.000 con una desviación estandar de $ 22.700.000, los caules tienen un promedio de modelo de 2020 y un kilometraje de 72.600 Km con una variación muy alta de 48.9000 KM, se observa en la figura 1 que la mayoria de los autos se agrupan entre precios de 22 a 95 Millones, en cuanto al modelo la mayor concentracion estan entre 2008 y 2022 siento este el modelo más predominante, en la figura 3 se observa que la mayor cantidad de autos son 0 Km y los demás no superan los 280.000 km, en cuanto a la relación del precio y el tipo de transmisión figura 7 los autos con transmisión automatica tienden ha ser más costosos que los mecanicos. Los municipios donde hay más ofertas de autos son Cali, medellin y manizales.
Finalmente la Figura 5 nos muestra que existe una relación lineal exponencial entre precio y modelo entre más reciente es el modelo tiene myor precio y en la figura 6 entre precio y kilometraje tambien existe una relación exponnecial inversa a medida que el kilometraje aunmenta disminuye el precio de venta.
# Cargar los datos
bdolxcar =read_excel("C:/Mesa/Informe regresión múltiple/olxcar.xlsx",
col_types = c("text", "numeric", "numeric",
"numeric", "text", "text", "text",
"text"))
bdolxcar=bdolxcar[,2:8]
bdolxcar$tranmision=factor(bdolxcar$tranmision,levels = c("Automática","Automática Secuencial","Mecánica"))
head(bdolxcar)
## # A tibble: 6 × 7
## precio modelo kilometraje tranmision ciudad Munic…¹ Depar…²
## <dbl> <dbl> <dbl> <fct> <chr> <chr> <chr>
## 1 136900000 2010 88000 Automática Astorga Medell… Antioq…
## 2 38000000 2013 66000 Automática San Carlos Tunjue… Bogotá
## 3 64500000 2018 74000 Automática Secuencial La Base Cali Valle …
## 4 88650000 2022 0 Automática Ciudad Cam… Cali Valle …
## 5 88650000 2022 0 Automática Ciudad Cam… Cali Valle …
## 6 88650000 2022 0 Automática Ciudad Cam… Cali Valle …
## # … with abbreviated variable names ¹Municipio, ²Departamento
## Análisis univariado y bivariado
attach(bdolxcar)
summary(bdolxcar)
## precio modelo kilometraje
## Min. : 6800000 Min. :1995 Min. : 0
## 1st Qu.: 35500000 1st Qu.:2012 1st Qu.: 34650
## Median : 42300000 Median :2015 Median : 79000
## Mean : 52012161 Mean :2015 Mean : 72564
## 3rd Qu.: 66000000 3rd Qu.:2019 3rd Qu.:103000
## Max. :169800000 Max. :2022 Max. :280000
## tranmision ciudad Municipio
## Automática :130 Length:273 Length:273
## Automática Secuencial: 20 Class :character Class :character
## Mecánica :121 Mode :character Mode :character
## NA's : 2
##
##
## Departamento
## Length:273
## Class :character
## Mode :character
##
##
##
table1(~precio+modelo+kilometraje,data=bdolxcar)
| Overall (N=273) |
|
|---|---|
| precio | |
| Mean (SD) | 52000000 (22700000) |
| Median [Min, Max] | 42300000 [6800000, 170000000] |
| modelo | |
| Mean (SD) | 2020 (4.55) |
| Median [Min, Max] | 2020 [2000, 2020] |
| kilometraje | |
| Mean (SD) | 72600 (48900) |
| Median [Min, Max] | 79000 [0, 280000] |
g1=ggplot(bdolxcar,aes(x=precio))+geom_histogram()+theme_bw()
g2=ggplot(bdolxcar,aes(x=modelo))+geom_histogram()+theme(axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5))
g3=ggplot(bdolxcar,aes(x=kilometraje))+geom_histogram()+theme_bw()
g4 =ggplot(bdolxcar, aes(Municipio))+geom_bar()+ggtitle("Figura 4. Hist de municipios")+theme(axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5))
fig1 =subplot(ggplotly(g1), ggplotly(g2), ggplotly(g3), titleX = TRUE,titleY = TRUE) %>% layout(title = " Fig 1. Hist de precios | Fig 2. Hist de modelo | Fig 3. Hist de Kilometraje")
fig1
ggplotly(g4)
g5=ggplot(data = bdolxcar,mapping = aes(x=modelo,y=precio))+geom_point()+geom_smooth()+ggtitle("Figura 5. Relación precio VS modelo")+theme_bw()
g6=ggplot(data = bdolxcar,mapping = aes(x= kilometraje,y=precio))+geom_point()+geom_smooth()+ggtitle("Figura 6. Relación precio VS kilometraje")+theme_bw()
g5
g6
g7=ggplot(bdolxcar,aes(x=tranmision,y=precio,fill=tranmision))+geom_boxplot()+ggtitle("Figura 7. Relación precio VS transmision")+theme_bw()
ggplotly(g7)
Se realiza el modelo con todas las variables y se evidencia que los p-value del modelo, transmisión y kilometraje son significativos siendo este ultimo el menos de los tres, de igual forma el municipio arrojo que no aporta mucho al modelo.
El \(R^2\) ajustado es medio alto (0,61) y significa que es capaz de explicar el modelo en un 61%
##bdolxcar2=subset(bdolxcar, tranmision !='Automática Secuencial')
mod_autos=lm(precio~modelo+kilometraje+tranmision+Municipio,data=bdolxcar)
summary(mod_autos)
##
## Call:
## lm(formula = precio ~ modelo + kilometraje + tranmision + Municipio,
## data = bdolxcar)
##
## Residuals:
## Min 1Q Median 3Q Max
## -67474245 -4122186 -63399 2330645 97884343
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -6.101e+09 7.248e+08 -8.417 5.13e-15 ***
## modelo 3.057e+06 3.586e+05 8.523 2.59e-15 ***
## kilometraje -6.426e+01 3.185e+01 -2.018 0.044862 *
## tranmisionAutomática Secuencial 2.863e+06 3.678e+06 0.778 0.437291
## tranmisionMecánica -8.021e+06 2.203e+06 -3.641 0.000339 ***
## MunicipioArmenia 3.638e+06 6.884e+06 0.529 0.597654
## MunicipioBarranquilla -1.245e+06 6.990e+06 -0.178 0.858857
## MunicipioBarrios Unidos -8.863e+06 9.639e+06 -0.920 0.358847
## MunicipioBello 8.749e+05 7.322e+06 0.119 0.905000
## MunicipioBosa -8.111e+05 1.511e+07 -0.054 0.957241
## MunicipioBucaramanga -8.688e+05 7.348e+06 -0.118 0.905991
## MunicipioCaldas -6.560e+06 1.515e+07 -0.433 0.665409
## MunicipioCali 3.844e+06 5.377e+06 0.715 0.475381
## MunicipioCesar 3.601e+06 1.513e+07 0.238 0.812141
## MunicipioChapinero -5.473e+06 9.733e+06 -0.562 0.574495
## MunicipioChía 1.955e+07 9.754e+06 2.004 0.046265 *
## MunicipioCórdoba -9.539e+06 1.513e+07 -0.630 0.529134
## MunicipioCúcuta -2.260e+06 9.832e+06 -0.230 0.818425
## MunicipioCundinamarca 7.203e+06 9.648e+06 0.747 0.456099
## MunicipioEngativá 1.664e+05 7.012e+06 0.024 0.981090
## MunicipioEnvigado 4.928e+06 8.218e+06 0.600 0.549363
## MunicipioFloridablanca 1.600e+06 7.684e+06 0.208 0.835280
## MunicipioFontibón 6.582e+06 1.506e+07 0.437 0.662466
## MunicipioFunza 2.750e+06 1.514e+07 0.182 0.856013
## MunicipioGirón -7.651e+05 1.508e+07 -0.051 0.959582
## MunicipioHuila -4.964e+06 1.536e+07 -0.323 0.746797
## MunicipioIbagué 2.213e+06 9.693e+06 0.228 0.819651
## MunicipioItagüí 7.658e+05 8.849e+06 0.087 0.931112
## MunicipioKennedy -5.747e+06 8.430e+06 -0.682 0.496140
## MunicipioLos Mártires -8.764e+06 1.133e+07 -0.773 0.440118
## MunicipioManizales 7.828e+05 6.212e+06 0.126 0.899836
## MunicipioMedellín 3.949e+06 5.535e+06 0.713 0.476333
## MunicipioMeta -1.026e+07 1.522e+07 -0.675 0.500681
## MunicipioNeiva 1.306e+06 7.679e+06 0.170 0.865136
## MunicipioPereira 3.635e+05 7.104e+06 0.051 0.959238
## MunicipioPiedecuesta 1.938e+05 1.506e+07 0.013 0.989742
## MunicipioPopayán 2.683e+06 1.512e+07 0.177 0.859311
## MunicipioPuente Aranda 1.706e+05 1.508e+07 0.011 0.990984
## MunicipioRafael Uribe Uribe -1.148e+05 1.507e+07 -0.008 0.993928
## MunicipioSabaneta -3.675e+06 1.517e+07 -0.242 0.808791
## MunicipioSan Cristóbal -5.835e+06 1.124e+07 -0.519 0.604204
## MunicipioSan Juan de Pasto 1.508e+06 1.129e+07 0.134 0.893838
## MunicipioSanta Fé 1.338e+06 1.516e+07 0.088 0.929783
## MunicipioSanta Marta 3.566e+06 1.530e+07 0.233 0.815929
## MunicipioSantander 3.151e+06 1.124e+07 0.280 0.779521
## MunicipioSuba 3.823e+06 7.677e+06 0.498 0.619048
## MunicipioSucre -1.193e+07 1.555e+07 -0.767 0.443873
## MunicipioTeusaquillo 1.369e+06 1.507e+07 0.091 0.927686
## MunicipioTunjuelito -9.806e+06 1.516e+07 -0.647 0.518495
## MunicipioUsaquén 1.964e+04 1.130e+07 0.002 0.998616
## MunicipioValle del Cauca -1.896e+06 8.843e+06 -0.214 0.830473
## MunicipioValledupar 3.396e+05 1.526e+07 0.022 0.982259
## MunicipioVillavicencio -1.317e+06 8.182e+06 -0.161 0.872241
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 14170000 on 218 degrees of freedom
## (2 observations deleted due to missingness)
## Multiple R-squared: 0.6862, Adjusted R-squared: 0.6114
## F-statistic: 9.17 on 52 and 218 DF, p-value: < 2.2e-16
par(mfrow=c(2,2))
plot(mod_autos)
En este caso se van a emplear la estrategia de stepwise mixto y el valor matemático empleado para determinar la calidad del modelo va a ser Akaike(AIC), con esto se reafirma la no significancia del municipio y finalmente el modelo queda \(precio_auto\)= \(-6.607e+09\)+ \(3.308e+06*modelo-4.856e+01*kilometraje+2.279e+06*tranmisionAutomática Secuencial-7.448e+06*tranmisionMecánica\), donde los valores \(\beta\) reafirman la evidencia exploratoria entre mnayor sea el modelo y el kilometraje menor el precio, los autos de transmision automatica tienen mayor precio que los de trasmision mecanica.
El modelo con las dos variables predictorias tiene un R2 alto (0,8253), lo que significa que es capaz de explicarlo en un 82,53% de la variablidad observada del peso, el p-value de las variables diametro y altura son significativos.
Validación del modelo
Media cero: Se cumple por defecto.
Varianza Constante: Se observa en la grafica 1 de residuales vs ajustados que el comportamiento es aleatorio aunque los valores se agrupan en la recta
Normalidad: Se observa en la grafica 2 que los datos se ajustan bien a la linea de normalidad en el qqplot
Independencia: 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.
step(object = mod_autos, direction = "both", trace = 1)
## Start: AIC=8971.82
## precio ~ modelo + kilometraje + tranmision + Municipio
##
## Df Sum of Sq RSS AIC
## - Municipio 48 3.7719e+15 4.7524e+16 8898.2
## <none> 4.3752e+16 8971.8
## - kilometraje 1 8.1696e+14 4.4569e+16 8974.8
## - tranmision 2 3.0537e+15 4.6806e+16 8986.1
## - modelo 1 1.4579e+16 5.8332e+16 9047.8
##
## Step: AIC=8898.23
## precio ~ modelo + kilometraje + tranmision
##
## Df Sum of Sq RSS AIC
## <none> 4.7524e+16 8898.2
## - kilometraje 1 6.2559e+14 4.8150e+16 8899.8
## - tranmision 2 3.8093e+15 5.1334e+16 8915.1
## + Municipio 48 3.7719e+15 4.3752e+16 8971.8
## - modelo 1 2.3988e+16 7.1512e+16 9007.0
##
## Call:
## lm(formula = precio ~ modelo + kilometraje + tranmision, data = bdolxcar)
##
## Coefficients:
## (Intercept) modelo
## -6.607e+09 3.308e+06
## kilometraje tranmisionAutomática Secuencial
## -4.856e+01 2.279e+06
## tranmisionMecánica
## -7.448e+06
mod_autos_final=lm(formula = precio ~ modelo + kilometraje + tranmision, data = bdolxcar)
summary(mod_autos_final)
##
## Call:
## lm(formula = precio ~ modelo + kilometraje + tranmision, data = bdolxcar)
##
## Residuals:
## Min 1Q Median 3Q Max
## -66583782 -4676144 -1293276 2609282 99806216
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -6.607e+09 5.768e+08 -11.455 < 2e-16 ***
## modelo 3.308e+06 2.855e+05 11.587 < 2e-16 ***
## kilometraje -4.856e+01 2.595e+01 -1.871 0.0624 .
## tranmisionAutomática Secuencial 2.279e+06 3.263e+06 0.698 0.4856
## tranmisionMecánica -7.448e+06 1.725e+06 -4.319 2.22e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 13370000 on 266 degrees of freedom
## (2 observations deleted due to missingness)
## Multiple R-squared: 0.6592, Adjusted R-squared: 0.6541
## F-statistic: 128.6 on 4 and 266 DF, p-value: < 2.2e-16
par(mfrow=c(2,2))
plot(mod_autos_final)
validar_modelo=function(datos,n){
var_predict=array(NA,n)
for (i in 1:n) {
datos_training=datos[-i,]
datos_test=datos[i,]
mod=lm(formula = precio ~ modelo + kilometraje + tranmision, data = datos_training)
var_predict[i]=predict(mod,list((modelo=datos_test$modelo),(kilometraje=datos_test$kilometraje),(tranmision=datos_test$tranmision)))
}
return(var_predict)
}
El modelo ajustado haciendo la validación cruzada con la técnica MAE se obtiene un 14,74% de probabilidad de error de predicción, lo que constituye un regular modelo predictivo y el RMSE 29,317% la proporción del error versus el promedio del precio observado
precio_predict=validar_modelo(bdolxcar,dim(bdolxcar)[1])
res=data.frame(bdolxcar,precio_predict)
##head(res)
MAE=mean(abs(res$precio-na.omit(res$precio_predict)))
print(paste("MAE= ",MAE/mean(res$precio)*100))
## [1] "MAE= 14.7418563122248"
RMSE=sqrt(mean((res$precio-na.omit(res$precio_predict))^2))
print(paste("RMSE",RMSE/mean(res$precio)*100))
## [1] "RMSE 29.3174676370685"
Discutir potenciales usos de modelo como herramienta practica (como monetizar los resultados de este modelo)
Dos tipos de moluscos A y B fueron sometidos a tres concentraciones distintas de agua de mar (100%, 75%, 50%) y se observo el consumo de oxigeno midiendo la proporcion de O2 por unidad peso seco del molusco.
Realice un analisis exploratorio que permita conocer como es el consumo de oxigeno en las distintas concentraciones de agua de mar. y si estas conclusiones son las mismas para cada tipo de molusco
require(ggplot2)
require(dplyr)
load("C:/Mesa/Informe regresión múltiple/moluscos.RData")
BaseMoluscos = data.frame(BD_moluscos)
BaseMoluscos$c_agua <- as.factor(BaseMoluscos$c_agua)
head(BaseMoluscos)
## c_agua molusco cons_o
## 1 100 A 7.16
## 2 100 A 8.26
## 3 100 A 6.78
## 4 100 A 14.00
## 5 100 A 13.60
## 6 100 A 11.10
moluscoA = filter(BaseMoluscos, BaseMoluscos$molusco == "A")
moluscoB = filter(BaseMoluscos, BaseMoluscos$molusco == "B")
summary(moluscoA$cons_o)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 5.20 7.18 9.74 10.00 11.28 18.80
summary(moluscoB$cons_o)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 1.800 5.685 8.060 8.609 10.700 17.700
Se aplica summary a la variable consumo de oxígeno (cons_o), para el caso de los moluscos tipo A tenemos una media de 9.74 mientras que su mínimo es de 5.20 y su máximo de 18.80. Para el caso de los moluscos tipo B tenemos una media de 8.060, un mínimo de 1.8000 y un máximo de 17.700
attach(BaseMoluscos)
summary(BaseMoluscos)
## c_agua molusco cons_o
## 50 :16 Length:48 Min. : 1.800
## 75 :16 Class :character 1st Qu.: 6.312
## 100:16 Mode :character Median : 9.700
## Mean : 9.305
## 3rd Qu.:11.232
## Max. :18.800
ggplot(BaseMoluscos,aes(x=cons_o))+geom_histogram()+labs(x="Consumo de Oxigeno",y="Count")+ggtitle("Consumo de Oxigeno por Tipo de Molusco")+facet_grid(~molusco)
En cuanto al consumo de oxígeno independientemente de la cantidad de concentración de agua, podemos observar que es muy similar para ambos tipos de moluscos, sin embargo, para los moluscos de Tipo A podemos ver que su consumo mínimo parte desde 5 mientras que para los moluscos de tipo B aproximadamente de 2
g3=ggplot(data = BaseMoluscos,mapping = aes(x=c_agua,y=cons_o))+geom_point()+geom_smooth()+ggtitle("Relacion Concretacion de agua VS Consumo de oxigeno")+theme_bw()
g3
g4 <- ggplot(BaseMoluscos, aes(x = c_agua, y = cons_o)) +
geom_point() + facet_grid(~molusco) +theme_bw() +ggtitle("Relacion Tipo de molusco VS Consumo de Oxgeno")+
labs(x = "Concentracion agua", y = "Consumo de Oxigeno")
g4
g1=ggplot(BaseMoluscos,aes(x= c_agua, y=cons_o, fill=molusco))+
geom_boxplot()+
theme_bw()+
geom_point(position=position_jitterdodge(),alpha=0.4)+
ggtitle("Consumo de Oxigeno por tipo de Molusco en diferentes consetracciones de Agua")+
labs(x="Consentracion de Agua",y="Consumo de Oxigeno")
g1
Podemos observar que el consumo de oxígeno es más alto cuando los moluscos se encuentran sometidos a una concentración de 50% de agua de mar y el consumo disminuye cuando la concentración de agua es del 75%.
interaction.plot(
x.factor = BaseMoluscos$c_agua,
trace.factor = BaseMoluscos$molusco,
response = BaseMoluscos$cons_o,
fun = median,
ylab = "Consumo de Oxigeno",
xlab = "Concetracion de Agua",
trace.label = "Tipo de Molusco",
lwd = 1
)
Con base a la gráfica podemos identificar el consumo de oxígeno por cada tipo de molusco dependiendo de la concentración de agua de mar a la que esten expuestos.
BaseMoluscos=na.omit(BaseMoluscos)
BaseMoluscos$molusco=as.character(BaseMoluscos$molusco)
g2=ggplot(BaseMoluscos,aes(x=c_agua,y=cons_o,color=molusco))+geom_point()+theme_bw()+geom_smooth()+ ggtitle("Consumo de Oxigeno por tipo de Molusco en diferentes consetracciones de Agua")+
labs(x="Consentracion de Agua",y="Consumo de Oxigeno")
g2
Con base a esta última gráfica podemos observar que independientemente de la concentración de agua a la que esten expuestos los tipos de moluscos el molusco de tipo A fue quien obtuvo mayor consumo de oxígeno en los 3 escenarios a los que fueron sometidos.
ModMoluscos = lm(formula = cons_o ~ c_agua + molusco, data=BaseMoluscos)
summary(ModMoluscos)
##
## Call:
## lm(formula = cons_o ~ c_agua + molusco, data = BaseMoluscos)
##
## Residuals:
## Min 1Q Median 3Q Max
## -5.1750 -1.9877 -0.7019 2.1244 6.1450
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 12.9463 0.8521 15.193 < 2e-16 ***
## c_agua75 -5.2581 1.0436 -5.038 8.49e-06 ***
## c_agua100 -3.5794 1.0436 -3.430 0.00132 **
## moluscoB -1.3913 0.8521 -1.633 0.10966
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.952 on 44 degrees of freedom
## Multiple R-squared: 0.3986, Adjusted R-squared: 0.3575
## F-statistic: 9.719 on 3 and 44 DF, p-value: 4.866e-05
Concentracion de agua de 75% (c_agua75) = 8.49e-06 Concentracion de agua de 100% (c_agua100) = 0.00132 moluscoB = 0.10966
Análisis: Como evidenciamos en los resultados del summary el valor que se obtuvo para p-value es 0.05, lo que sugeriría que existe una relación lineal entre las variables, pero solo para los casos donde las concentraciones de agua son 75% y 100%, es decir que el modelo se ve afectado solo cuando hay un cambio en las concentraciones de agua, mientras que si se cambia el tipo de molusco este no tendría un cambio
anova(ModMoluscos)
## Analysis of Variance Table
##
## Response: cons_o
## Df Sum Sq Mean Sq F value Pr(>F)
## c_agua 2 230.82 115.408 13.2457 3.14e-05 ***
## molusco 1 23.23 23.227 2.6658 0.1097
## Residuals 44 383.37 8.713
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Al analizar los resultados obtenidos aplicando post anova, podemos darnos cuenta que la variable moluscos, como se indicó previamente no es significativa dentro del modelo, ya que su p-value es de 3.14e-05
por otra parte, podemos evidenciar que la variable c_agua es la más representativa, ya que si la concentración de agua varia, los resultados del modelo también se verán afectados. p-value es de 0.1097
require(agricolae)
Compara= LSD.test(ModMoluscos,c("c_agua","molusco"))
Compara
## $statistics
## MSerror Df Mean CV t.value LSD
## 8.712897 44 9.304792 31.72303 2.015368 2.974442
##
## $parameters
## test p.ajusted name.t ntr alpha
## Fisher-LSD none c_agua:molusco 6 0.05
##
## $means
## cons_o std r LCL UCL Min Max Q25 Q50 Q75
## 100:A 9.93625 2.747976 8 7.833002 12.039498 6.78 14.00 7.9850 9.295 11.7250
## 100:B 7.40625 2.844076 8 5.303002 9.509498 3.68 11.60 5.7225 6.140 10.1000
## 50:A 12.17500 3.090178 8 10.071752 14.278248 9.74 18.80 10.3100 11.110 12.5000
## 50:B 12.32625 3.517909 8 10.223002 14.429498 6.38 17.70 10.0575 12.850 14.5000
## 75:A 7.89000 2.739578 8 5.786752 9.993248 5.20 13.20 6.0775 7.180 8.8925
## 75:B 6.09500 2.739108 8 3.991752 8.198248 1.80 9.96 4.8300 5.595 7.3425
##
## $comparison
## NULL
##
## $groups
## cons_o groups
## 50:B 12.32625 a
## 50:A 12.17500 a
## 100:A 9.93625 ab
## 75:A 7.89000 bc
## 100:B 7.40625 bc
## 75:B 6.09500 c
##
## attr(,"class")
## [1] "group"
Usamos la libreria agricolae la cual ordena los valores de p de acuerdo con las medias del grupo,
ModMoluscos
##
## Call:
## lm(formula = cons_o ~ c_agua + molusco, data = BaseMoluscos)
##
## Coefficients:
## (Intercept) c_agua75 c_agua100 moluscoB
## 12.946 -5.258 -3.579 -1.391
par(mfrow=c(2,2))
plot(ModMoluscos)
Análisis Supuestos
- Linealidad: Tiene correlacion Lineal.
- Normalidad: No se cumple el supuesto de normalidad nuestro p-value = 0.04342 por lo que es menor a 0.05
shapiro.test(ModMoluscos$residuals)
##
## Shapiro-Wilk normality test
##
## data: ModMoluscos$residuals
## W = 0.9509, p-value = 0.04342
- Homocedatidad (Varianza Constante): No se evidencia Homocedaticidad por lo que se cosidera un modelo heterocedatico el cual no tiene una varianza constante y cuya P corresponde a 0.41514
library(car)
ncvTest(ModMoluscos)
## Non-constant Variance Score Test
## Variance formula: ~ fitted.values
## Chisquare = 0.6640401, Df = 1, p = 0.41514
- Independencia: los errores son independientes p-value = 0.2163
lmtest::dwtest(ModMoluscos)
##
## Durbin-Watson test
##
## data: ModMoluscos
## DW = 1.8991, p-value = 0.2163
## alternative hypothesis: true autocorrelation is greater than 0
Estudiar la relación entre ciertas características del suelo y la producción de biomasa (gr) de una planta forrajera natural se obtuvieron 45 muestras en diferentes ambientes, y en cada muestra se estimó la biomasa (respuesta Y) y se registraron las características (covariables X) del suelo en el que crecía (pH, Salinidad, Zinc y Potasio).
Del análisis preliminar se pueden extraer las siguientes conclusiones:
Las variables que tienen una mayor relación lineal con la Biomasa son: PH (r= 0.928), Zinc (r= -0.781)
Tambien se oberva que PH y Zinc están fuerte mente correlacionados (r= -0,720), al igual de sanialidaad con ZINC por lo que se sugiere dejar solo Zinc para el analisis
load("C:/Mesa/Informe regresión múltiple/Salinidad.RData")
head(Salinidad)
## Biomasa pH Salinidad Zinc Potasio
## 1 765.280 5.00 33 16.4524 1441.67
## 2 954.017 4.70 35 13.9852 1299.19
## 3 827.686 4.20 32 15.3276 1154.27
## 4 755.072 4.40 30 17.3128 1045.15
## 5 896.176 5.55 33 22.3312 521.62
## 6 1422.836 5.50 33 12.2778 1273.02
library(corrplot)
round(cor(x = Salinidad, method = "pearson"), 3)
## Biomasa pH Salinidad Zinc Potasio
## Biomasa 1.000 0.928 -0.067 -0.781 -0.073
## pH 0.928 1.000 -0.045 -0.720 0.032
## Salinidad -0.067 -0.045 1.000 -0.427 -0.020
## Zinc -0.781 -0.720 -0.427 1.000 0.079
## Potasio -0.073 0.032 -0.020 0.079 1.000
multi.hist(x = Salinidad, dcol = c("blue", "red"), dlty = c("dotted", "solid"),
main = "")
ggpairs(Salinidad, lower = list(continuous = "smooth"),
diag = list(continuous = "barDiag"), axisLabels = "none")
Se realiza el metodo mixto con todas las variables, como predictoras y realizando la selección con los mejores predictores con el calculo del R2 ajustado para ver el porcentaje de explicación del mismo y lo valores de p-value si son significativos tando del modelo, como las variables predictoras y finalmente se relaiza la selección de los mejores predictores con la medición Akaike(AIC).
El modelo con las cuatro variables predictorias tiene un R2 ajustado alto (0,9154), lo que significa que es capaz de explicarlo en un 91,54% de la variablidad observada de la biomasa, el p-value de las variables PH,Sanilidad y Zinc son significativos, la variable potación su p-value no arrojo significancia sim embargo en la prueba de selección AIC se mantiene como predictora del modelo
Validación del modelo
Media cero: Se cumple por defecto.
Varianza Constante: Se observa en la grafica 1 de residuales vs ajustados que el comportamiento es aleatorio.
Normalidad: Se observa en la grafica 2 que los datos se ajustan bien a la linea de normalidad en el qqplot
Independencia: 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.
Modelo: \(Biomasa= 1492.808+262.883 *pH -33.500*Salinidad-28.973*Zinc-0.115*Potasio\)
mod_sanilidad=lm(Biomasa~ .,data = Salinidad)
summary(mod_sanilidad)
##
## Call:
## lm(formula = Biomasa ~ ., data = Salinidad)
##
## Residuals:
## Min 1Q Median 3Q Max
## -293.98 -88.83 -9.48 88.20 387.27
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1492.8076 453.6013 3.291 0.002091 **
## pH 262.8829 33.7304 7.794 1.51e-09 ***
## Salinidad -33.4997 8.6525 -3.872 0.000391 ***
## Zinc -28.9727 5.6643 -5.115 8.20e-06 ***
## Potasio -0.1150 0.0819 -1.404 0.167979
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 158.9 on 40 degrees of freedom
## Multiple R-squared: 0.9231, Adjusted R-squared: 0.9154
## F-statistic: 120 on 4 and 40 DF, p-value: < 2.2e-16
step(object = mod_sanilidad, direction = "both", trace = 1)
## Start: AIC=460.84
## Biomasa ~ pH + Salinidad + Zinc + Potasio
##
## Df Sum of Sq RSS AIC
## <none> 1009974 460.84
## - Potasio 1 49785 1059759 461.01
## - Salinidad 1 378486 1388460 473.17
## - Zinc 1 660588 1670562 481.49
## - pH 1 1533665 2543639 500.41
##
## Call:
## lm(formula = Biomasa ~ pH + Salinidad + Zinc + Potasio, data = Salinidad)
##
## Coefficients:
## (Intercept) pH Salinidad Zinc Potasio
## 1492.808 262.883 -33.500 -28.973 -0.115
par(mfrow = c(2, 2))
plot(mod_sanilidad)