El presente documento contiene dos ejercicios de simulación de regresiones lineales multiples con diferentes tematicas que se presentan a continuación:

1.Moluscos

Dos tipos de moluscos A y B fueron sometidos a tres concentraciones distintas de agua de mar (100%, 75% y 50%) y observó el consumo de oxígeno midiendo la proporción de O2 por unidad de peso seco del molusco.

load("moluscos.RData")

BD_moluscos$c_agua <- factor(BD_moluscos$c_agua)
BD_moluscos$molusco <- factor(BD_moluscos$molusco)
BD_moluscos
# A tibble: 48 × 3
   c_agua molusco cons_o
   <fct>  <fct>    <dbl>
 1 100    A         7.16
 2 100    A         8.26
 3 100    A         6.78
 4 100    A        14   
 5 100    A        13.6 
 6 100    A        11.1 
 7 100    A         8.93
 8 100    A         9.66
 9 100    B         6.14
10 100    B         6.14
# … with 38 more rows

1.a. Análisis exploratorio

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.

range_tipo <- ggplot(BD_moluscos, mapping=aes(), main="Distribución Producción") + geom_boxplot(aes(y=cons_o, x=molusco), fill='#67B7D1') + 
  labs(title ="Rangos de concentración de oxigeno por tipo de molusco", x="Tipo de molusco", y="Concentración Oxigeno")

ggplotly(range_tipo)

Nuestra primera grafica nos permite evidencia que el rango de los moluscos tipo B es significamente mayor que rango tipo A, aún así la media de concentración de oxigeno de tipo A es mayor que la del tipo B.

range_agua <- ggplot(BD_moluscos, mapping=aes(y=cons_o, x=c_agua)) + geom_boxplot(varwidth = TRUE) + labs(title ="Rangos de concentración de oxigeno por concentración de agua de mar", x="Concentración de agua de mar", y="Concentración Oxigeno")
ggplotly(range_agua)

No se evidencia una relación de concentración de oxigeno por concentración de agua de mar, debido a que de 100% a 75% se ve un decrecimiento del oxigeno mientras que del 75% o del 100% si se ve un aumento respecto al 50%. El 50% supera significativamente la concentración de oxigeno de un 75% al 100%, para el caso del 75% de concentración de agua del mar, el 50% evidencia casi el doble de oxigeno.

histo_tipo <- ggplot(data=BD_moluscos, aes(x=cons_o, fill=molusco)) +
    geom_histogram(binwidth=3, position="dodge") + labs(title ="Frecuencia niveles de oxigeno por tipo", x="Concentración de oxigeno",
                                                                      y="Conteo")

ggplotly(histo_tipo)

Entre los niveles de concentración de oxigeno entre 0 a 6.75, se observa que existe mayor frecuencia de moluscos tipo B, en adelante se evidencia un sesgo a la izquierda donde los moluscos tipo A se encuentran en mayor frecuencia.

histo_agua <- ggplot(BD_moluscos, aes(x=cons_o, fill=molusco)) +
    geom_histogram(position="identity", alpha=0.5, binwidth=3)+
    facet_grid(. ~ c_agua) + labs(title ="Frecuencia niveles de oxigeno por tipo por concentración de agua de mar", x="Concentración de oxigeno",
                                                                      y="Conteo")

ggplotly(histo_agua)

En todas las concentraciones de agua de mar se evidencia que el molusco tipo A esta logrando mayor concentración de oxigeno que los tipo B. Adicionalmente podemos observar que la mayor cantidad de moluscos concentrados se encuentran en la concentración de agua 75%. Y la mayor cantidad de oxgieno se encuentra en el 50%, sin embargo, el 100% parece mostrar un buen comportamiento también en concentración de oxigeno pero realmente es dado por los moluscos tipo A.

La conclusión de esta exploración inicial muestra como ha mostrado que a menor concentración de agua de mar, el molusco tipo B esta presentando mayor niveles de oxigeno, mientras que a mayor niveles de concentración de agua de mar el moluscto tipo A esta presentando mayores niveles de oxigeno. Por otro lado, los niveles de concentración de oxigeno son mucho más estables en los tipo A que los tipo B, el tipo tiene una reacción más fuerte al nivel de agua de mar.

Correlación agua de mar vs oxigeno tipo A

filter1 <- filter(BD_moluscos,BD_moluscos$molusco=="A")
cor(filter1$cons_o, strtoi(filter1$c_agua))
[1] -0.2856287

Correlación agua de mar vs oxigeno tipo B

filter2 <- filter(BD_moluscos,BD_moluscos$molusco=="B")
cor(filter2$cons_o, strtoi(filter2$c_agua))
[1] -0.5126334

Asi lo confirma la correlación de Person, quien indica que hay una baja correlación negativa en los moluscos tipo A, es decir, a mayor concentración de agua de mar, menor oxigeno se esta obteniendo, sin embargo, esta mostrando una menor correlación. Mientras que en los moluscos tipo B, esta correlación llega a un -51%.

1.b.Experimentos

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.

Modelo usado

yij = MU + Bi + Lambdaj + (BLambda)ij + Eij

model1 <- lm(BD_moluscos$cons_o ~ BD_moluscos$c_agua + BD_moluscos$molusco + (BD_moluscos$c_agua*BD_moluscos$molusco))

summary(model1)

Call:
lm(formula = BD_moluscos$cons_o ~ BD_moluscos$c_agua + BD_moluscos$molusco + 
    (BD_moluscos$c_agua * BD_moluscos$molusco))

Residuals:
   Min     1Q Median     3Q    Max 
-5.946 -1.736 -0.710  2.237  6.625 

Coefficients:
                                           Estimate Std. Error t value Pr(>|t|)
(Intercept)                                 12.1750     1.0466  11.633 1.02e-14
BD_moluscos$c_agua75                        -4.2850     1.4800  -2.895  0.00599
BD_moluscos$c_agua100                       -2.2387     1.4800  -1.513  0.13787
BD_moluscos$moluscoB                         0.1513     1.4800   0.102  0.91909
BD_moluscos$c_agua75:BD_moluscos$moluscoB   -1.9462     2.0931  -0.930  0.35777
BD_moluscos$c_agua100:BD_moluscos$moluscoB  -2.6813     2.0931  -1.281  0.20722
                                              
(Intercept)                                ***
BD_moluscos$c_agua75                       ** 
BD_moluscos$c_agua100                         
BD_moluscos$moluscoB                          
BD_moluscos$c_agua75:BD_moluscos$moluscoB     
BD_moluscos$c_agua100:BD_moluscos$moluscoB    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 2.96 on 42 degrees of freedom
Multiple R-squared:  0.4226,    Adjusted R-squared:  0.3539 
F-statistic: 6.149 on 5 and 42 DF,  p-value: 0.0002324

Modelo Anova

anova1 <- aov(lm(BD_moluscos$cons_o ~ BD_moluscos$c_agua + BD_moluscos$molusco + (BD_moluscos$c_agua*BD_moluscos$molusco)))
summary(anova1)
                                       Df Sum Sq Mean Sq F value   Pr(>F)    
BD_moluscos$c_agua                      2  230.8  115.41  13.171 3.63e-05 ***
BD_moluscos$molusco                     1   23.2   23.23   2.651    0.111    
BD_moluscos$c_agua:BD_moluscos$molusco  2   15.4    7.68   0.876    0.424    
Residuals                              42  368.0    8.76                     
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
plot(anova1)

etaSquared(anova1)
                                           eta.sq eta.sq.part
BD_moluscos$c_agua                     0.36211526  0.38544681
BD_moluscos$molusco                    0.03643950  0.05936774
BD_moluscos$c_agua:BD_moluscos$molusco 0.02409168  0.04005632

Comparacion del modelo Anova con LSD - Least significant difference

TukeyHSD(anova1)
  Tukey multiple comparisons of means
    95% family-wise confidence level

Fit: aov(formula = lm(BD_moluscos$cons_o ~ BD_moluscos$c_agua + BD_moluscos$molusco + (BD_moluscos$c_agua * BD_moluscos$molusco)))

$`BD_moluscos$c_agua`
            diff        lwr       upr     p adj
75-50  -5.258125 -7.8007168 -2.715533 0.0000289
100-50 -3.579375 -6.1219668 -1.036783 0.0039342
100-75  1.678750 -0.8638418  4.221342 0.2551550

$`BD_moluscos$molusco`
        diff       lwr       upr     p adj
B-A -1.39125 -3.115713 0.3332134 0.1109762

$`BD_moluscos$c_agua:BD_moluscos$molusco`
                diff          lwr        upr     p adj
75:A-50:A   -4.28500  -8.70331211  0.1333121 0.0619482
100:A-50:A  -2.23875  -6.65706211  2.1795621 0.6582697
50:B-50:A    0.15125  -4.26706211  4.5695621 0.9999983
75:B-50:A   -6.08000 -10.49831211 -1.6616879 0.0023423
100:B-50:A  -4.76875  -9.18706211 -0.3504379 0.0277338
100:A-75:A   2.04625  -2.37206211  6.4645621 0.7369714
50:B-75:A    4.43625   0.01793789  8.8545621 0.0485578
75:B-75:A   -1.79500  -6.21331211  2.6233121 0.8281801
100:B-75:A  -0.48375  -4.90206211  3.9345621 0.9994669
50:B-100:A   2.39000  -2.02831211  6.8083121 0.5936863
75:B-100:A  -3.84125  -8.25956211  0.5770621 0.1209401
100:B-100:A -2.53000  -6.94831211  1.8883121 0.5334845
75:B-50:B   -6.23125 -10.64956211 -1.8129379 0.0017250
100:B-50:B  -4.92000  -9.33831211 -0.5016879 0.0212746
100:B-75:B   1.31125  -3.10706211  5.7295621 0.9478135
plot(TukeyHSD(anova1))

df<-df.residual(anova1)
MSerror<-deviance(anova1)/df
out <- with(BD_moluscos,LSD.test(cons_o,c_agua,df,MSerror))
out
$statistics
   MSerror Df     Mean      CV  t.value      LSD
  8.762171 42 9.304792 31.8126 2.018082 2.112028

$parameters
        test p.ajusted name.t ntr alpha
  Fisher-LSD      none c_agua   3  0.05

$means
      cons_o      std  r       LCL       UCL  Min  Max    Q25    Q50     Q75
100  8.67125 3.000940 16  7.177821 10.164679 3.68 14.0  6.140  8.595 10.5750
50  12.25062 3.199643 16 10.757196 13.744054 6.38 18.8 10.085 11.455 14.5000
75   6.99250 2.804093 16  5.499071  8.485929 1.80 13.2  5.200  6.430  8.7675

$comparison
NULL

$groups
      cons_o groups
50  12.25062      a
100  8.67125      b
75   6.99250      b

attr(,"class")
[1] "group"
plot(out,variation="IQR")

Ahora graficamos las las variables de la base de datos de moluscos

ggplot () +
aes (x = BD_moluscos$c_agua, y = BD_moluscos$cons_o, color = BD_moluscos$molusco) +
 geom_line (aes (grupo = BD_moluscos$molusco))  +
 geom_point () +
 labs (title = ' Puntuación del examen por horas estudiadas y técnica de estudio ',
     color = ' Técnica ',
     x = ' Horas dedicadas a estudiar ',
     y = ' Puntuación media del examen ') +
 theme_minimal ()

anova_pareado <- Anova(model1, idata = BD_moluscos(cons_o),
                       idesign = ~ cons_o, type = "III")
summary(anova_pareado, multivariate = F)
     Sum Sq                Df          F value              Pr(>F)       
 Min.   :   0.0915   Min.   : 1.0   Min.   :  0.01044   Min.   :0.00000  
 1st Qu.:  15.3563   1st Qu.: 1.0   1st Qu.:  0.65982   1st Qu.:0.01639  
 Median :  73.4943   Median : 2.0   Median :  2.53506   Median :0.22282  
 Mean   : 328.5597   Mean   : 9.6   Mean   : 35.10437   Mean   :0.34118  
 3rd Qu.: 368.0112   3rd Qu.: 2.0   3rd Qu.: 36.97961   3rd Qu.:0.54762  
 Max.   :1185.8450   Max.   :42.0   Max.   :135.33690   Max.   :0.91909  
                                    NA's   :1           NA's   :1

Se puede evidenciar que Si quieremos conocer el cambio de oxigeno, en los moluscos debemos centrarnos en la concentracion de agua debido a que se esta mostrando una diferencia en la poblacion por medio de esta variable

Finalmente graficamos la interacion

interaction.plot(x.factor = BD_moluscos$c_agua, # variable del eje x 
                 trace.factor = BD_moluscos$molusco, # variable para líneas 
                 response = BD_moluscos$cons_o, # variable del eje y 
                 fun = median, # métrica para trazar
                 ylab = "concentracion de oxigeno",
                 xlab = "concentracion de agua de mar",
                 col = c ("red", "blue"),
                 lty = 1, # tipo de línea 
                 lwd = 2, # ancho de línea
                 trace.label = "Tipo Molusco")

### si un molusco tiene que tener mas concentracion oxigeno, su concentracion de agua de mar debe ser igual a 50, ### si un molusco tiene que tener menor concentracion oxigeno, su concentracion de agua debe ser igual a 75

shapiro.test(anova1$res)

    Shapiro-Wilk normality test

data:  anova1$res
W = 0.95824, p-value = 0.08571
shapiro.test(anova1$res)

    Shapiro-Wilk normality test

data:  anova1$res
W = 0.95824, p-value = 0.08571
#resi2=
#bartlett.test(anova1$res)

Se observa que se cumple es supuesto de normalida debido a que no cumple la hipotesis nula que es equivalente a que no es normal

2.Salinidad

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)

.

load(file="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

2.a. Análisis de correlación

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).

corrplot(cor(Salinidad),# Matriz de correlación
         method = "number", # Método para el gráfico de correlación
         type = "lower",    # Estilo del gráfico (también "upper" y "lower")
         diag = TRUE,      # Si TRUE (por defecto), añade la diagonal
         tl.col = "black", # Color de las etiquetas
         bg = "white",# Color de fondo
         is.corr=T,
         insig = "label_sig",
         title = "",       # Título
         col = NULL)

Las correlaciones más fuertes entre estas covariables y Biomasa son pH y Zinc, y existe una baja correlación con salinidad y postasio. Sin embargo, se debe realizar un análisis de multicolinealidad debido a que pH y Zinc también tienen una correlación negativa fuerte y el Zinc y la Salinidad una relación de casi 43%. Se debe analizar las covariables que se incluiran en el modelo. Zinc es una variable que esta correlaccionada con varias covariables

pairs(~Biomasa  + pH + Salinidad + Zinc + Potasio, data = Salinidad)

A pesar de que Zinc presenta una correlación fuerte se evidencia un poco más de dispersión que pH, salinidad evidencia fuerte aleatoriedad.

2.b. Regresión lineal multiple

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.

summary(salinidadfullmodel <- lm(Biomasa ~ ., data = Salinidad))

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

El modelo con todas sus covariables tiene muy buen ajuste y variables significativas excepto por Potasio, sin embargo, revisaremos si tenemos alguna mejora.

slm1 <- stepAIC(salinidadfullmodel, trace=TRUE, direction="backward")
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
slm1$anova
Stepwise Model Path 
Analysis of Deviance Table

Initial Model:
Biomasa ~ pH + Salinidad + Zinc + Potasio

Final Model:
Biomasa ~ pH + Salinidad + Zinc + Potasio


  Step Df Deviance Resid. Df Resid. Dev      AIC
1                         40    1009974 460.8448

Según lo anterior el modelo no puede realizar mejoras al cambiar sus covariables. Revisemos si se cumplen los supuestos.

#Distribución normal

residuales_sml1 <- slm1$residuals
shapiro.test(residuales_sml1)

    Shapiro-Wilk normality test

data:  residuales_sml1
W = 0.96586, p-value = 0.2036

Cumplimos con el supuesto de normalidad porque nuestro valor p es 0.20, indica que no hay significancia para rechaza la hipotesis nula equivalente a normal.

#Varianza constante

lmtest::bptest(slm1) 

    studentized Breusch-Pagan test

data:  slm1
BP = 5.4945, df = 4, p-value = 0.2402

Cumplimos el supuesto de homocedasticidad porque nuestro valor p es 0.24, indica que hay noy significancia para rechaza la hipotesis nula equivalente a varianza constante, esto es validado por la gráfica inicial donde no se nota una apertura del canal de los puntos.

#Errores independientes

lmtest::dwtest(slm1) 

    Durbin-Watson test

data:  slm1
DW = 1.6647, p-value = 0.06483
alternative hypothesis: true autocorrelation is greater than 0

Hasta el momento parece ser un modelo perfecto, pero la relación que existe entre la covariable Zinc y pH persiste, por tal razón y para no incurrir a errores en la interpretación del efecto que tiene cada covariable procederemos a eliminar la variable con menor correlación que es Zinc.

summary(salinidadnomulti <- lm(Biomasa ~ pH + Salinidad + Potasio, data = Salinidad))

Call:
lm(formula = Biomasa ~ pH + Salinidad + Potasio, data = Salinidad)

Residuals:
    Min      1Q  Median      3Q     Max 
-608.83  -93.37   -6.08  112.37  380.96 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept) -511.8776   290.0872  -1.765   0.0851 .  
pH           405.0152    24.2890  16.675   <2e-16 ***
Salinidad     -3.9855     8.1904  -0.487   0.6291    
Potasio       -0.1906     0.1023  -1.863   0.0697 .  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 201.9 on 41 degrees of freedom
Multiple R-squared:  0.8728,    Adjusted R-squared:  0.8635 
F-statistic: 93.76 on 3 and 41 DF,  p-value: < 2.2e-16

Nuestro modelo tiene un menor rendimiento, sin embargo, aseguramos que el modelo sea más confiable, pero ahora nuestras variables salinidad y Potasio no tienen tanta significacia o explican el modelo.

summary(salinidaduno <- lm(Biomasa ~ pH, data = Salinidad))

Call:
lm(formula = Biomasa ~ pH, data = Salinidad)

Residuals:
    Min      1Q  Median      3Q     Max 
-566.28  -89.26  -19.42  142.42  413.28 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  -780.18     117.99  -6.612  4.7e-08 ***
pH            404.08      24.72  16.346  < 2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 205.7 on 43 degrees of freedom
Multiple R-squared:  0.8614,    Adjusted R-squared:  0.8582 
F-statistic: 267.2 on 1 and 43 DF,  p-value: < 2.2e-16

Solo existe una mejora de 0.01% en el R2, si se mantiene el modelo simple con solo 1 covariable el pH, se puede obtener un rendimiento similar. Ahora verifiquemos los supuestos:

#Distribución normal

residuales_sml <- salinidaduno$residuals
shapiro.test(residuales_sml)

    Shapiro-Wilk normality test

data:  residuales_sml
W = 0.9681, p-value = 0.2468

Cumplimos con el supuesto de normalidad porque nuestro valor p es 0.24, indica que no hay significancia para rechaza la hipotesis nula equivalente a normal.

#Varianza constante

lmtest::bptest(salinidaduno) 

    studentized Breusch-Pagan test

data:  salinidaduno
BP = 0.63672, df = 1, p-value = 0.4249

Cumplimos el supuesto de homocedasticidad porque nuestro valor p es 0.43, indica que hay noy significancia para rechaza la hipotesis nula equivalente a varianza constante, esto es validado por la gráfica inicial donde no se nota una apertura del canal de los puntos.

#Errores independientes

lmtest::dwtest(salinidaduno) 

    Durbin-Watson test

data:  salinidaduno
DW = 1.2166, p-value = 0.002001
alternative hypothesis: true autocorrelation is greater than 0