Regresion Lineal Multiple
Marco Javier Peñaloza Pérez
2022-09-16
Taller #2
Pregunta 1. Dos tipos de moluscos A y B fueron sometidos a tres concentraciones distintas de agua de mar (100%, 75% y 50%) y se observó el consumo de oxígeno midiendo la proporción de O2 por unidad de peso seco del molusco.
Punto a) Realice un análisis exploratorio que permita conocer como es el consumo de oxígeno en las distintas concentraciones de agua de mar. y si estas conclusiones son las mismas para cada tipo de molusco.
Solución
Paso 1: Cargar la base de datos:
load("C:/Users/mjpenalozap/Desktop/moluscos.RData")
BD_moluscos$c_agua = as.character(BD_moluscos$c_agua)
BD_moluscos$c_agua <- factor(BD_moluscos$c_agua, levels=unique(BD_moluscos$c_agua))
head(BD_moluscos)## 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.10dim(BD_moluscos)## [1] 48 3Se observa que la base de datos tiene 48 observaciones con tres variables: Concentración del Agua (Cualitativa), Tipo de Moluscos (Cualitativa) y Consumo de Oxígeno (Cuantitativa).
Paso 2: Resumen numérico.
La siguiente tabla presenta un resumen sobre los distintos valores que presentan cada una de las variables:
library(table1)##
## Attaching package: 'table1'## The following objects are masked from 'package:base':
##
## units, units<-table1(~c_agua+cons_o+molusco,data = BD_moluscos)| Overall (N=48) |
|
|---|---|
| c_agua | |
| 100 | 16 (33.3%) |
| 75 | 16 (33.3%) |
| 50 | 16 (33.3%) |
| cons_o | |
| Mean (SD) | 9.30 (3.68) |
| Median [Min, Max] | 9.70 [1.80, 18.8] |
| molusco | |
| A | 24 (50.0%) |
| B | 24 (50.0%) |
summary(BD_moluscos)## c_agua molusco cons_o
## 100:16 Length:48 Min. : 1.800
## 75 :16 Class :character 1st Qu.: 6.312
## 50 :16 Mode :character Median : 9.700
## Mean : 9.305
## 3rd Qu.:11.232
## Max. :18.800Se observa que los moluscos tienen igual número de tipos A y B iguales a 24, el número total ha sido igualmente distribuidos entre los tres niveles de concentración de agua de mar (50%, 75% y 100%). También se evidencia el resumen de 5 valores para la variable cuantitativa Consumo de Oxígeno, la cual tiene sus valores mínimo y máximo en 1.8 y 18.8 respectivamente, los cuartiles son \(Q_1=6.31\), \(Q_2=9.7\) y \(Q_3=11.2\), finalmente la media es de 9.3, todos estos datos son tomados de manera general.
Ahora se clasifican los datos según el tipo de Molusco:
library(RcmdrMisc)## Loading required package: car## Loading required package: carData## Loading required package: sandwichnumSummary(BD_moluscos$cons_o, statistics=c("mean", "sd", "quantiles","cv"), type=c("2", "1", "3"), quantiles=c(0, .25, .5, .75, 1), BD_moluscos$molusco)## mean sd cv 0% 25% 50% 75% 100% data:n
## A 10.000417 3.268661 0.3268525 5.2 7.180 9.74 11.2825 18.8 24
## B 8.609167 4.002435 0.4649039 1.8 5.685 8.06 10.7000 17.7 24Se puede observar que los moluscos tipo A presentan mayor promedio que los tipo B, incluso con una dispersión más baja, también se aprecian los valores de los diferentes cuartiles.
Ahora se clasifican los datos según los niveles de concentración del agua de mar:
numSummary(BD_moluscos$cons_o, statistics=c("mean", "sd", "quantiles","cv"), type=c("2", "1", "3"), quantiles=c(0, .25, .5, .75, 1), BD_moluscos$c_agua)## mean sd cv 0% 25% 50% 75% 100% data:n
## 100 8.67125 3.000940 0.3460792 3.68 6.140 8.595 10.5750 14.0 16
## 75 6.99250 2.804093 0.4010144 1.80 5.200 6.430 8.7675 13.2 16
## 50 12.25062 3.199643 0.2611820 6.38 10.085 11.455 14.5000 18.8 16Se puede resaltar que los moluscos expuestos a un nivel de agua salada de 50% presentan un mayor promedio que los demás, incluso con una dispersión menor del 26%, se observa que tienen mayor consumo de oxigeno, se aprecia también los diferentes valores de los cuartiles para cada nivel.
Paso 3: Realizar diagrama de cajas.
Primero se realizan los diagramas considerando dos variables:
library(ggplot2)
caja1=ggplot(BD_moluscos,aes(x=c_agua, y=cons_o), fill=class)+geom_boxplot(alpha=0.3)+theme_bw()+ xlab("Concentración de Agua") + ylab("Consumo de Oxígeno") + ggtitle("Consumo de Oxígeno según la concentración de agua de mar") + theme(plot.title = element_text(hjust = 0.5))
caja1
caja2=ggplot(BD_moluscos,aes(x=molusco, y=cons_o))+geom_boxplot()+theme_bw()+xlab("Tipo de Molusco") + ylab("Consumo de Oxígeno") + ggtitle("Consumo de Oxígeno por Tipo de Molusco") + theme(plot.title = element_text(hjust = 0.5)) + stat_summary()
caja2## No summary function supplied, defaulting to `mean_se()`
Se realiza el diagrama considerando simultáneamente todas las variables involucrada:
library(ggplot2)
caja=ggplot(BD_moluscos,aes(x=c_agua, y=cons_o, fill=molusco))+geom_boxplot()+theme_bw()+xlab("Concentración de Agua") + ylab("Consumo de Oxígeno") + ggtitle("Diagrama de Caja") + theme(plot.title = element_text(hjust = 0.5))
caja+ geom_boxplot() + scale_fill_brewer(palette = "Greens")
El Boxplot indica que el mayor consumo de oxígeno se da a un nivel de concentración de 50% de agua de mar y el menor a 75%. Los moluscos tipo A tienen mayor consumo de oxígeno para concentraciones de 100% y 75%, mientras que para el 50% de concentración los moluscos tipo B presentan mayor consumo.
Punto b) Estime el modelo de diseño de experimentos el cual permita evaluar el efecto de la concentración de agua de mar y los tipos de molusco sobre el consumo de oxigeno. Interprete los coeficientes del modelo, el valor p y realice un post anova de considerarlo necesario para los factores.
Solución: Parte 1: Construcción del modelo:
modelo1= lm((BD_moluscos$cons_o)~BD_moluscos$molusco+BD_moluscos$c_agua)
summary(modelo1)##
## Call:
## lm(formula = (BD_moluscos$cons_o) ~ BD_moluscos$molusco + BD_moluscos$c_agua)
##
## 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) 9.3669 0.8521 10.993 3.33e-14 ***
## BD_moluscos$moluscoB -1.3913 0.8521 -1.633 0.10966
## BD_moluscos$c_agua75 -1.6788 1.0436 -1.609 0.11486
## BD_moluscos$c_agua50 3.5794 1.0436 3.430 0.00132 **
## ---
## 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-05Se establece la ecuación del valor estimado del consumo de oxígeno:
Consumo de oxígeno = 9.36-1.39(Molusco B)-1.67(Concentración 75%)+3.5794(Concentración 50%)
El modelo estimado presenta el intercepto \(β_0=9.36\), lo cual indica que en ausencia de las variables predictoras Tipos de Moluscos y Concentración de agua de mar, el consumo de oxígeno es de 9.36, el cual presenta un p-valor que hace que su incidencia sea significativa sobre la variable dependiente. El coeficiente \(β_1=-1.39\) indica la disminución que se da en dicho factor en caso de ser un Molusco tipo B, lo cual no ocurre con los moluscos tipo A, debido al modelo estimado se dice que tiene poca incidencia sobre la variable consumo de oxígeno. El coeficiente \(β_2=-1.67\) tampoco tiene relevancia en el modelo, corresponde a la dismunicón en la consumo de oxígeno en casi de tener una concentración de 75% en agua de mar, finalmente, el coeficiente \(β_3=3.57\), es el aumento que se da en la consumo de oxígeno cuano se tiene una concentración de 100%, según el modelo, el coeficiente si tiene significancia para la ecuación de regresión.
Parte 2: Validación de Supuestos.
Errores Residuales
ei = modelo1$residuals
round(mean(ei),4)## [1] 0Se concluye que los errores tiene media igual a cero.
Prueba de normalidad para los residuales
shapiro.test(ei)##
## Shapiro-Wilk normality test
##
## data: ei
## W = 0.9509, p-value = 0.04342Se concluye que los residuales no siguen una distribución normal debido al p-valor que está por debajo del 5%.
Varianza de los residuales
lmtest::bptest(modelo1)##
## studentized Breusch-Pagan test
##
## data: modelo1
## BP = 1.5424, df = 3, p-value = 0.6725Se concluye que la varianza es constante.
Gráficos de Residuales
par(mfrow=c(1,2))
plot(modelo1, which=1)
plot(modelo1, which=2)
Test Anova
anova_del_modelo= aov(modelo1)
summary(anova_del_modelo)## Df Sum Sq Mean Sq F value Pr(>F)
## BD_moluscos$molusco 1 23.2 23.23 2.666 0.11
## BD_moluscos$c_agua 2 230.8 115.41 13.246 3.14e-05 ***
## Residuals 44 383.4 8.71
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Se concluye que l variable Concentración del agua si es significativa para el modelo, mientras que el tipo de molusco no.
Análisis PostAnova
anova(modelo1)## Analysis of Variance Table
##
## Response: (BD_moluscos$cons_o)
## Df Sum Sq Mean Sq F value Pr(>F)
## BD_moluscos$molusco 1 23.23 23.227 2.6658 0.1097
## BD_moluscos$c_agua 2 230.82 115.408 13.2457 3.14e-05 ***
## Residuals 44 383.37 8.713
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1attach(BD_moluscos)
library(agricolae)
comparar_moluscos= LSD.test(modelo1,c("c_agua","molusco"))
comparar_moluscos## NULLSe concluye que para la variable dependiente las concentraciones de 50 y 100 se obtienen resultados similares, mientras que con 50 se consiguen mayores consumos de oxigeno, las concentraciones de agua de 50 se obtiene mayor consumo de oxígeno para los moluscos tipo A y tipo B. También podemos observar que con concentraciones de agua de 75 ambos tipos de moluscos bajan su consumo de oxígeno.
Pregunta 2. Para 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).
Punto a) Realice un análisis de correlaciones que permita identificar de manera bivariada las relaciones entre las covariables y la respuesta (incluir coeficiente de correlación e interpretaciones).
Solución
Paso 1: Cargar la base de datos:
load("C:/Users/mjpenalozap/Desktop/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.02dim(Salinidad)## [1] 45 5Se tiene una base de datos con 45 observaciones y variables cuantitativas.
Paso 2: Ánalisi Exploratorio
grafico1=ggplot(Salinidad,aes(x=Salinidad,y=Biomasa)) + geom_point() + geom_smooth()+ ggtitle("Salinidad VS Biomasa")+theme(plot.title = element_text(hjust = 0.5))
grafico1## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
grafico2=ggplot(Salinidad,aes(x=pH,y=Biomasa)) + geom_point() + geom_smooth() + ggtitle("Biomasa vs PH")+theme(plot.title = element_text(hjust = 0.5))
grafico2## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
grafico3=ggplot(Salinidad,aes(x=Potasio,y=Biomasa)) + geom_point() + geom_smooth()+ ggtitle("Biomasa vs Potasio")+theme(plot.title = element_text(hjust = 0.5))
grafico3## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
grafico4=ggplot(Salinidad,aes(x=Zinc,y=Biomasa)) + geom_point() + geom_smooth()+ ggtitle("Biomasa vs Zinc")+theme(plot.title = element_text(hjust = 0.5))
grafico4## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
Paso 3: Cuadro de correlaciones entre las variables.
library(foreign)
library(Factoshiny)## Loading required package: FactoMineR## Loading required package: shiny## Loading required package: FactoInvestigatelibrary(factoextra)## Welcome! Want to learn more? See two factoextra-related books at https://goo.gl/ve3WBa##
## Attaching package: 'factoextra'## The following object is masked from 'package:agricolae':
##
## hcutlibrary(corrplot)## corrplot 0.92 loadedlibrary(biplotbootGUI)
library(dynBiplotGUI)##
## *** ------------
## *** dynBiplotGUI ***
##
## Available languages:
## Spanish-"es", English-"en", French-"fr", Portuguese-"pt".
## Use: dynBiplot("en") to launch the interface.library(psych)##
## Attaching package: 'psych'## The following objects are masked from 'package:ggplot2':
##
## %+%, alpha## The following object is masked from 'package:RcmdrMisc':
##
## reliability## The following object is masked from 'package:car':
##
## logitpairs.panels(Salinidad)
mcor = cor(Salinidad)
data.frame(mcor)## Biomasa pH Salinidad Zinc Potasio
## Biomasa 1.00000000 0.92810235 -0.06657756 -0.78146249 -0.07319518
## pH 0.92810235 1.00000000 -0.04458851 -0.72046995 0.03236212
## Salinidad -0.06657756 -0.04458851 1.00000000 -0.42663388 -0.01963288
## Zinc -0.78146249 -0.72046995 -0.42663388 1.00000000 0.07877268
## Potasio -0.07319518 0.03236212 -0.01963288 0.07877268 1.00000000Punto b) Estime el modelo de regresión lineal múltiple para explicar la biomasa en función de las covariables e interprete el valor p, los coeficientes de las variables significativas y el coeficiente R2.
Solución
Paso 1: Estimación del modelo.
attach(Salinidad)## The following object is masked _by_ .GlobalEnv:
##
## Salinidadmodelo2= lm(Biomasa~pH+Salinidad$Salinidad+Zinc+Potasio)
summary(modelo2)##
## Call:
## lm(formula = Biomasa ~ pH + Salinidad$Salinidad + Zinc + Potasio)
##
## 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$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-16Paso 2: Validación de Supuestos.
Media de los residuales
e2 = modelo2$residuals
round(mean(e2),4)## [1] 0Prueba de normalidad de los residuales
shapiro.test(e2)##
## Shapiro-Wilk normality test
##
## data: e2
## W = 0.96586, p-value = 0.2036Constancia de la varianza
lmtest::bptest(modelo2)##
## studentized Breusch-Pagan test
##
## data: modelo2
## BP = 5.4945, df = 4, p-value = 0.2402Validación Gráfica
par(mfrow=c(1,2))
plot(modelo2, which=1)
plot(modelo2, which=2)
anova(modelo2)## Analysis of Variance Table
##
## Response: Biomasa
## Df Sum Sq Mean Sq F value Pr(>F)
## pH 1 11310631 11310631 447.9573 < 2.2e-16 ***
## Salinidad$Salinidad 1 8352 8352 0.3308 0.5684
## Zinc 1 752176 752176 29.7899 2.725e-06 ***
## Potasio 1 49785 49785 1.9718 0.1680
## Residuals 40 1009974 25249
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Conclusiones:
La ecuación del modelo estimado es: \(biomasa=1492.80+262.88Ph−33.49Salinidad−28.97Zinc-0.11Potasio\)
El intercepto equivale a \(β_0=1492.80\) que corresponde al valor de la biomasa en ausencia de las variables explicativas.
El coeficiente \(β_1=262.88\) corresponde al aumento de la biomasa por acad aumento unitario en el Ph.
El coeficiente \(β_2=-33.49\) corresponde a la disminución de la biomasa por cada aumento unitario en la salinidad.
El coeficiente \(β_3=-28.97\) corrresponde a la disminución de la biomasa por cada aumento unitario en el zinc.
El coeficiente \(β_4=-0.11\) corrresponde a la disminución de la biomasa por cada aumento unitario en el Potasio. . Según el modelo, todos los coeficientes son estadísticamente significativos, por tener un p-valor por debajo del 5%. (Excepto el de Potasio)
El modelo presenta un buen ajuste, ya que el coeficiente de determinación estima una proporción de explicación del 92.31%.
No se cumple con el supuesto de normalidad de los residuales.
El análisis nde Anova indica que las variables que tienen mayor influencia sobre la biomasa son Ph y Zinc.