1 . Conjunto de datos

La base de datos Mississippi proviene del paquete SASmixed en R, el cual contiene conjuntos de datos utilizados para la modelización de efectos mixtos en estudios ambientales, biomédicos y de otras disciplinas.

Este conjunto de datos en particular fue recopilado para analizar la concentración de nitrógeno en el río Misisipi, un sistema fluvial de gran importancia en los Estados Unidos. La contaminación por nitrógeno en los ríos puede provenir de diversas fuentes, incluyendo actividades agrícolas, descargas industriales y aguas residuales. Su monitoreo es clave para la gestión de la calidad del agua y la prevención de problemas ambientales como la eutrofización.

Los datos están estructurados en 37 observaciones y tres variables, que representan factores ambientales y mediciones de la calidad del agua en diferentes puntos de muestreo.

De manera que, la estructura de los Datos, está dada por:

Variable Tipo de Dato Descripción
influent Factor ordenado Punto de muestreo en el río (niveles: 3 < 5 < 2 < 1 < 4 < 6). Indica diferentes ubicaciones donde se midieron las concentraciones de nitrógeno.
y Numérica Concentración de nitrógeno en el agua en cada punto de muestreo.
Type Factor Clasificación del tipo de muestra tomada en el río, con tres niveles: 1, 2, 3. Puede representar diferentes fuentes de agua o métodos de muestreo. 
# Cargar paquetes necesarios
# install.packages("SASmixed")
library(SASmixed)
data(Mississippi)
head(Mississippi)
##   influent  y Type
## 1        1 21    2
## 2        1 27    2
## 3        1 29    2
## 4        1 17    2
## 5        1 19    2
## 6        1 12    2
library(ggplot2)
library(lme4)
## Loading required package: Matrix
library(dplyr)
## 
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
str(Mississippi)
## 'data.frame':    37 obs. of  3 variables:
##  $ influent: Factor w/ 6 levels "1","2","3","4",..: 1 1 1 1 1 1 1 1 1 2 ...
##  $ y       : num  21 27 29 17 19 12 29 20 20 21 ...
##  $ Type    : Factor w/ 3 levels "1","2","3": 2 2 2 2 2 2 2 2 2 2 ...
##  - attr(*, "ginfo")=List of 7
##   ..$ formula     :Class 'formula'  language y ~ 1 | influent
##   .. .. ..- attr(*, ".Environment")=<environment: R_GlobalEnv> 
##   ..$ order.groups: logi TRUE
##   ..$ FUN         :function (x)  
##   ..$ outer       :Class 'formula'  language ~Type
##   .. .. ..- attr(*, ".Environment")=<environment: R_GlobalEnv> 
##   ..$ inner       : NULL
##   ..$ labels      :List of 1
##   .. ..$ y: chr "Nitrogen concentration in Mississippi River"
##   ..$ units       :List of 1
##   .. ..$ y: chr "(ppm)"
data("Mississippi")

2 . Análisis exploratorio de datos

En primer lugar se observará la distribución de las variables, la relación entre ellas y detectar posibles patrones o valores atípicos.

Distribución de la Variable Dependiente (y):visualizar cómo se distribuyen las concentraciones de nitrógeno:

ggplot(Mississippi, aes(x = y)) +
  geom_histogram(bins = 10, fill = "blue", alpha = 0.6, color = "black") +
  theme_minimal() +
  labs(title = "Distribución de Concentraciones de Nitrógeno", x = "y", y = "Frecuencia")

summary(Mississippi$y)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    2.00   14.00   20.00   20.97   28.00   42.00

La concentración de nitrógeno en el río Misisipi varía entre un mínimo de 2 y un máximo de 42. El valor promedio registrado es de aproximadamente 20.97, mientras que la mediana es 20, indicando una distribución equilibrada de los datos.

2.1 . Boxplot de Niveles de Nitrógeno (y) según el afluente

Los diagramas de caja y bigotes muyestran que hay diferencias entre los grupos, sugiere que el afluente o tipo de afluente podrían ser factores significativos en un modelo. Se observan diferencias en cuanto al los niveles de nitrogeno en el tipo 3.

str(Mississippi)
## 'data.frame':    37 obs. of  3 variables:
##  $ influent: Factor w/ 6 levels "1","2","3","4",..: 1 1 1 1 1 1 1 1 1 2 ...
##  $ y       : num  21 27 29 17 19 12 29 20 20 21 ...
##  $ Type    : Factor w/ 3 levels "1","2","3": 2 2 2 2 2 2 2 2 2 2 ...
##  - attr(*, "ginfo")=List of 7
##   ..$ formula     :Class 'formula'  language y ~ 1 | influent
##   .. .. ..- attr(*, ".Environment")=<environment: R_GlobalEnv> 
##   ..$ order.groups: logi TRUE
##   ..$ FUN         :function (x)  
##   ..$ outer       :Class 'formula'  language ~Type
##   .. .. ..- attr(*, ".Environment")=<environment: R_GlobalEnv> 
##   ..$ inner       : NULL
##   ..$ labels      :List of 1
##   .. ..$ y: chr "Nitrogen concentration in Mississippi River"
##   ..$ units       :List of 1
##   .. ..$ y: chr "(ppm)"
ggplot(Mississippi, aes(x = influent, y = y, fill = influent)) +
  geom_boxplot() +
  theme_minimal() +
  labs(title = "Distribución de los niveles de nitrógeno y según el afluente")

# Boxplot de y según Type
ggplot(Mississippi, aes(x = Type, y = y, fill = Type)) +
  geom_boxplot() +
  theme_minimal() +
  labs(title = "Distribución de los niveles de nitrógeno según el tipo de Afluente")

2.2 . Distribuición los niveles de influent dentro de cada Type

Muestra cómo están distribuidos los niveles de influent dentro de cada Type.

ggplot(Mississippi, aes(x = influent, fill = Type)) +
  geom_bar(position = "dodge") +
  theme_minimal() +
  labs(title = "Frecuencia de influent según Type")

2.3 . Gráfica “Relación entre afluente, Tipo e y(concentración de nitrogéno)

La gráfica “Relación entre influent, Type y y”: Permite identificar si hay interacción entre los factores en su relación con y.

ggplot(Mississippi, aes(x = influent, y = y, color = Type)) +
  geom_jitter(width = 0.2) +
  theme_minimal() +
  labs(title = "Relación entre influent, Type y y")

2.4 . Efecto en y(concentración de nitrógeno) según el afluente

Si se sospecha que el afluente(influent) tiene un efecto en y(concentración de nitrógeno), se puede observar las tendencias:

ggplot(Mississippi, aes(x = as.numeric(influent), y = y)) +
  geom_point() +
  geom_smooth(method = "lm", se = FALSE, color = "red") +
  theme_minimal() +
  labs(title = "Tendencia entre influente y", x = "Influent", y = "Concentración de Nitrógeno")
## `geom_smooth()` using formula = 'y ~ x'

3 . Algunos modelos

Dado que el conjunto de datos contienen información sobre concentraciones de nitrógeno en el río Misisipi con las variables “influent” (factor ordenado), “y” (concentración) y “Type” (factor con tres niveles), se puede considerar algunos enfoques en modelos de efectos fijos, aleatorios y mixtos.

3.1 . Modelo de Efectos Fijos

En este caso, se considera que todas las variables categóricas (factores) tienen efectos sistemáticos y no aleatorios sobre la variable respuesta y(niveles de niitrogeno).

Un modelo de regresión lineal con efectos fijos sería:

\[y_{ij} = \beta_0 + \beta_1 \text{influent}_i + \beta_2 \text{Type}_j + \varepsilon_{ij}\] donde:

  • \(y_{ij}\) es la concentración de nitrógeno en la observación j en el influente i.
  • \(\beta_0\) es el intercepto general.
  • \(\beta_1 influent_i\) representa el efecto fijo del influente en la concentración de nitrógeno.
  • \(\beta_2 Type_j\) representa el efecto fijo del tipo de muestra \((Type)\).
  • \(\varepsilon_{ij} \sim N(0, \sigma^2)\) es el término de error aleatorio, que se asume normalmente distribuido con media 0 y varianza \(\sigma^2\).
#Mississippi$Type <- relevel(Mississippi$Type, ref = "Type3")  # Usa Type3 como referencia , PARA CAMABIAR LA CATEGORÍA DE REFERENCIA
lm_modelo <- lm(y ~ influent + Type, data = Mississippi)
summary(lm_modelo)
## 
## Call:
## lm(formula = y ~ influent + Type, data = Mississippi)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -11.8571  -3.5000  -0.5556   4.6000  13.6000 
## 
## Coefficients: (2 not defined because of singularities)
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   21.556      2.175   9.911 3.97e-11 ***
## influent2     -7.698      3.288  -2.341 0.025827 *  
## influent3     -4.756      3.639  -1.307 0.200928    
## influent4      2.944      3.439   0.856 0.398447    
## influent5     -7.156      3.639  -1.966 0.058292 .  
## influent6     14.844      3.639   4.079 0.000293 ***
## Type2             NA         NA      NA       NA    
## Type3             NA         NA      NA       NA    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 6.525 on 31 degrees of freedom
## Multiple R-squared:  0.5933, Adjusted R-squared:  0.5277 
## F-statistic: 9.044 on 5 and 31 DF,  p-value: 2.202e-05

3.2 . Modelo de Efectos Aleatorios

Si consideramos que influent y/o Type representan muestras de una población más amplia y no niveles fijos, se pueden modelar como efectos aleatorios.

Un modelo con “influent” como efecto aleatorio sería:

\[y_{ij} = \beta_0 + \gamma_i + \varepsilon_{ij}\]

donde:

  • \(y_{ij}\) es la concentración de nitrógeno en la observación j dentro del influente i.
  • \(\beta_0\) es el intercepto general.
  • \(\gamma_i \sim N(0, \sigma^2_{\gamma})\) es un efecto aleatorio para cada nivel de influent, distribuido normalmente con media 0 y varianza \(\sigma^2_{\gamma}\).
  • \(\varepsilon_{ij} \sim N(0, \sigma^2)\) es el error aleatorio, también normalmente distribuido con media 0 y varianza \(\sigma^2\).

Implementación en R con lme4:

library(lme4)
efale_modelo <- lmer(y ~ 1 + (1 | influent), data = Mississippi)
summary(efale_modelo)
## Linear mixed model fit by REML ['lmerMod']
## Formula: y ~ 1 + (1 | influent)
##    Data: Mississippi
## 
## REML criterion at convergence: 252.4
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -1.91444 -0.53646 -0.03217  0.83714  1.95824 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  influent (Intercept) 63.32    7.958   
##  Residual             42.66    6.531   
## Number of obs: 37, groups:  influent, 6
## 
## Fixed effects:
##             Estimate Std. Error t value
## (Intercept)   21.223      3.429   6.189

3.3 . Modelo Mixto (Efectos Fijos y Aleatorios)

Un modelo mixto considera efectos fijos y aleatorios al mismo tiempo.

Si tratamos Type como efecto fijo y influent como efecto aleatorio:

\[y_{ijk} = \beta_0 + \beta_1 \text{Type}j + \gamma_i + \varepsilon_{ijk}\]

donde:

  • \(y_{ijk}\) es la concentración de nitrógeno en la observación \(k\) dentro del tipo \(j\) y el afluente \(i\).
  • \(\beta_0\) es el intercepto general.
  • \(\beta_1 Type_j\) representa el efecto fijo del tipo de muestra (Type).
  • \(\gamma_i \sim N(0, \sigma^2_{\gamma})\) es el efecto aleatorio de influent, que sigue una distribución normal con media 0 y varianza \(\sigma^2_{\gamma}\).
  • \(\varepsilon_{ijk} \sim N(0, \sigma^2)\) es el término de error aleatorio, también normalmente distribuido.
modelo_mixto <- lmer(y ~ Type + (1 | influent), data = Mississippi)
summary(modelo_mixto)
## Linear mixed model fit by REML ['lmerMod']
## Formula: y ~ Type + (1 | influent)
##    Data: Mississippi
## 
## REML criterion at convergence: 234.5
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.0877 -0.6392 -0.0257  0.7055  2.0191 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  influent (Intercept) 14.97    3.869   
##  Residual             42.51    6.520   
## Number of obs: 37, groups:  influent, 6
## 
## Fixed effects:
##             Estimate Std. Error t value
## (Intercept)   15.600      3.426   4.554
## Type2          4.338      4.324   1.003
## Type3         20.800      5.934   3.505
## 
## Correlation of Fixed Effects:
##       (Intr) Type2 
## Type2 -0.792       
## Type3 -0.577  0.457

3.4 . Modelo anidado

Si se considera al tipo de afluente(Type) que representa un grupo dentro de afluente(influent). Un modelo anidado se podría usar cuando una variable categórica está contenida dentro de otra jerárquicamente, por ejemplo, si Type representa diferentes tipos de mediciones dentro de cada influent.

Un modelo anidado con efectos mixtos se podría representar como:

\[y_{ijk} = \beta_0 + \beta_1 \cdot influent_i + u_i + v_{ij} + \varepsilon_{ijk}\]

donde:

  • \(y_{ijk}\): Concentración de nitrógeno para la observación k, en el tipo j dentro del punto de muestreo i.
  • \(\beta_0\): Intercepto general del modelo.
  • \(\beta_1 \cdot influent_i\): Efecto fijo de influent (ubicación de muestreo en el río).
  • \(u_i \sim N(0, \sigma_u^2)\): Efecto aleatorio del punto de muestreo (influent).
  • \(v_{ij} \sim N(0, \sigma_v^2)\): Efecto aleatorio de Type dentro de cada influent.
  • \(\varepsilon_{ijk} \sim N(0, \sigma^2)\): Error aleatorio residual.

En el caso de usar el paquete nlme, que permite definir modelos mixtos de la siguiente manera:

  • \(y \sim influent\): Modela la concentración de nitrógeno en función del factor influent.
  • \(random = ~1 | influent/Type\): \(Type\) está anidado dentro de \(influent\), lo que significa que la variabilidad en \(Type\) está dentro de cada \(influent\).
  • \(~1\) indica un intercepto aleatorio, permitiendo que cada \(influent\) tenga su propio nivel base.
library(nlme)
## 
## Attaching package: 'nlme'
## The following object is masked from 'package:dplyr':
## 
##     collapse
## The following object is masked from 'package:lme4':
## 
##     lmList
modelo_anidado <- lme(y ~ influent, random = ~1 | influent/Type, data = Mississippi)
summary(modelo_anidado)
## Warning in pt(-abs(tVal), fDF): NaNs produced
## Linear mixed-effects model fit by REML
##   Data: Mississippi 
##        AIC      BIC    logLik
##   233.0256 245.9315 -107.5128
## 
## Random effects:
##  Formula: ~1 | influent
##         (Intercept)
## StdDev:     7.00665
## 
##  Formula: ~1 | Type %in% influent
##         (Intercept) Residual
## StdDev:     7.00665 6.524839
## 
## Fixed effects:  y ~ influent 
##                 Value Std.Error DF    t-value p-value
## (Intercept) 21.555556  10.14479 31  2.1247916  0.0417
## influent2   -7.698413  14.39392  0 -0.5348379     NaN
## influent3   -4.755556  14.47818  0 -0.3284637     NaN
## influent4    2.944444  14.42909  0  0.2040631     NaN
## influent5   -7.155556  14.47818  0 -0.4942304     NaN
## influent6   14.844444  14.47818  0  1.0252978     NaN
##  Correlation: 
##           (Intr) infln2 infln3 infln4 infln5
## influent2 -0.705                            
## influent3 -0.701  0.494                     
## influent4 -0.703  0.496  0.493              
## influent5 -0.701  0.494  0.491  0.493       
## influent6 -0.701  0.494  0.491  0.493  0.491
## 
## Standardized Within-Group Residuals:
##         Min          Q1         Med          Q3         Max 
## -1.81723140 -0.53641168 -0.08514471  0.70499821  2.08434252 
## 
## Number of Observations: 37
## Number of Groups: 
##           influent Type %in% influent 
##                  6                  6

De manera similiar, el modelo anterior se puede ajustar con la librería lme4, en donde:

  • \(y ~ influent\): Modelo el efecto fijo de influent sobre la concentración de nitrógeno. _ \((1 | influent/Type)\):\(influent/Type\) indica que Type está anidado dentro de influent.Esto modela efectos aleatorios para influent y para Type dentro de influent
library(lme4)
modelo_anidado2 <- lmer(y ~ influent + (1 | influent:Type), data = Mississippi)
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## unable to evaluate scaled gradient
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Hessian is numerically singular: parameters are not uniquely determined
summary(modelo_anidado2)
## Linear mixed model fit by REML ['lmerMod']
## Formula: y ~ influent + (1 | influent:Type)
##    Data: Mississippi
## 
## REML criterion at convergence: 215
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -1.81723 -0.53641 -0.08514  0.70500  2.08434 
## 
## Random effects:
##  Groups        Name        Variance Std.Dev.
##  influent:Type (Intercept)  3.881   1.970   
##  Residual                  42.574   6.525   
## Number of obs: 37, groups:  influent:Type, 6
## 
## Fixed effects:
##             Estimate Std. Error t value
## (Intercept)   21.556      2.935   7.345
## influent2     -7.698      4.310  -1.786
## influent3     -4.756      4.583  -1.038
## influent4      2.944      4.426   0.665
## influent5     -7.156      4.583  -1.561
## influent6     14.844      4.583   3.239
## 
## Correlation of Fixed Effects:
##           (Intr) infln2 infln3 infln4 infln5
## influent2 -0.681                            
## influent3 -0.640  0.436                     
## influent4 -0.663  0.451  0.425              
## influent5 -0.640  0.436  0.410  0.425       
## influent6 -0.640  0.436  0.410  0.425  0.410
## optimizer (nloptwrap) convergence code: 0 (OK)
## unable to evaluate scaled gradient
##  Hessian is numerically singular: parameters are not uniquely determined

Las advertencias en la salida, indican que hay problemas de identificación en los parámetros del modelo. Esto puede deberse a:

  • Multicolinealidad: Algunas variables están altamente correlacionadas.
  • Datos insuficientes: Hay pocas observaciones en cada grupo de influent:Type.
  • Modelo demasiado complejo: Demasiados efectos aleatorios para la cantidad de datos.

En cuanto a los efectos aleatorios:

Varianza de influent:Type = 3.881 → La variabilidad entre niveles de influent:Type es baja.

  • Varianza residual = 42.574, indica que hay una gran cantidad de variabilidad no explicada por el modelo.
  • Desviación estándar de influent:Type = 1.970, en este caso cuánto varían las medias de cada Type dentro de influent.

En cuanto a los efectos fijos:

Intercepto (36.400): Representa la concentración promedio de nitrógeno en la categoría de referencia de influent.

  • Efectos de in bbfluent:
  • influent1 reduce la concentración en -14.844 unidades respecto al intercepto.
  • influent2, influent3, influent5 tienen efectos negativos significativos, indicando que esos lugares tienen menores concentraciones.

En cuanto a la evaluación de ajuste del modelo:

  • Normalidad en los errores, se cumple.
qqnorm(resid(modelo_anidado))
qqline(resid(modelo_anidado), col = "red")

- Homocedasticidad (varianza constante):

plot(fitted(modelo_anidado), resid(modelo_anidado))
abline(h = 0, col = "red")

- Comparación con un Modelo Más Simple:

modelo_simple <- lmer(y ~ influent + (1 | influent), data = Mississippi)
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## unable to evaluate scaled gradient
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Hessian is numerically singular: parameters are not uniquely determined
anova(modelo_simple, modelo_anidado)
## Warning in anova.merMod(modelo_simple, modelo_anidado): additional arguments
## ignored
## Analysis of Variance Table
##          npar Sum Sq Mean Sq F value
## influent    5 1291.6  258.33  6.0678

En este caso el p-valor es alto, el modelo anidado no es significativamente mejor y podemos simplificarlo.

4 Referencias

  • Littell, R. C., Milliken, G. A., Stroup, W. W., Wolfinger, R. D., & Schabenberger, O. (2006). SAS for Mixed Models (2nd ed.). SAS Institute Inc.

  • Pinheiro, J., & Bates, D. (2003). SASmixed: Mixed Effects Models in S and S-PLUS. R package version 0.2-3. Disponible en: https://cran.r-project.org/package=SASmixed

citation("SASmixed")
## To cite package 'SASmixed' in publications use:
## 
##   Littell Ob, Milliken, Stroup, Wolfinger, Bates mbD, Maechler M,
##   Bolker B, Walker S (2014). _SASmixed: Data sets from "SAS System for
##   Mixed Models"_. R package version 1.0-4,
##   <https://CRAN.R-project.org/package=SASmixed>.
## 
## A BibTeX entry for LaTeX users is
## 
##   @Manual{,
##     title = {SASmixed: Data sets from "SAS System for Mixed Models"},
##     author = {Original by Littell and {Milliken} and {Stroup} and {Wolfinger} and modifications by Douglas Bates and Martin Maechler and Ben Bolker and Steven Walker},
##     year = {2014},
##     note = {R package version 1.0-4},
##     url = {https://CRAN.R-project.org/package=SASmixed},
##   }
## 
## ATTENTION: This citation information has been auto-generated from the
## package DESCRIPTION file and may need manual editing, see
## 'help("citation")'.