Análisis estadístico de datos biológicos en R

Seminario

Autor

José Francisco Calvo Sendín
jfcalvo@um.es | https:/webs.um.es/jfcalvo

Departamento

Ecología e Hidrología, Facultad de Biología

Fecha de publicación

25 de septiembre de 2023

Fecha de actualización

7 de noviembre de 2023

Preparación

R es un potente software estadístico y gráfico que puede descargarse gratuitamente desde la página https://cran.r-project.org/. A través de la web de la red CRAN (The Comprehensive R Archive Network) puede encontrarse la documentación necesaria para iniciarse y profundizar en el uso de R, así como de sus bibliotecas o paquetes de funciones, creadas por una amplísima comunidad de programadores e investigadores de todos los ámbitos de las ciencias.

Para realizar los ejercicios propuestos en el seminario se necesitan varios paquetes de R instalados por defecto (MASS, mgcv, nlme y nnet), y otros que requieren una instalación adicional (effects, aods3, lme4, car y MuMIn). Para instalarlos usaremos:

install.packages( c( "effects", "aods3", "lme4", "car", "MuMIn" ) )

Importación de datos

Los datos que se utilizarán en el seminario, tanto simulados como reales, están importados y guardados en un archivo de datos de R (aedbr.RData), disponible en un servidor de la Universidad de Murcia. Una vez iniciado R debemos cargar este archivo con la función load, indicando la dirección URL:

load( url( "https://webs.um.es/jfcalvo/aedbr.RData" ) )

Usando ls() podemos ver los datos y funciones cargados. Utiliza la función info para obtener información sobre los objetos; por ejemplo: info( "aliens ) o info( "MonteCarlo" ). Teclea info() para más información.

Como alternativa podemos descargar el archivo en nuestro ordenador accediendo con un navegador de internet a https://webs.um.es/jfcalvo/aedbr.RData. A continuación, una vez iniciado R, es conveniente cambiar el directorio de trabajo a dicha carpeta (> Archivo > Cambiar dir...), de forma que podamos acceder fácilmente al archivo descargado: > Archivo > Cargar área de trabajo... También podemos clicar sobre el archivo para que se abra directamente en el programa.

Para aprender a importar (y exportar) datos en R, en la web del proyecto CRAN podemos encontrar una extensa guía sobre los paquetes y funciones que facilitan estas tareas.

Existen numerosos métodos de importación de archivos de datos en R, dependiendo de su formato. Entre los más habituales figuran los archivos de texto sin formato o texto plano (ASCII, American Standard Code for Information Interchange), los archivos con formato csv (datos separados por comas, o por puntos y comas), y los archivos de hojas de cálculo como Excel.

La importación de archivos de texto en R se puede realizar con la función read.table. Un aspecto muy importante para tener en cuenta es el separador de decimales de nuestros datos. Aunque en español se utiliza la coma, lo habitual es usar puntos (así es en R). Por tanto, necesitamos generalmente usar la función con el argumento dec = ",". Si la primera columna de datos no lleva cabecera, la función interpreta que corresponde al nombre de las filas, y la importación se realiza adecuadamente; en caso contrario tendríamos que introducir el argumento header = TRUE.

Para leer archivos csv se utilizan las funciones read.csv (cuando el separador de datos es la coma), y read.csv2 (cuando el separador de datos es el punto y coma).

En el caso de archivos Excel, la importación requiere de paquetes adicionales, como xlsx, que cuenta con la función read.xlsx.

ANOVA, regresión lineal y ANCOVA

Comenzaremos realizando análisis de la varianza (ANOVA) utilizando los datos simulados contenidos en el objeto ejemplo:

ejemplo

Podemos conocer la estructura de este objeto y las características de sus variables con:

str( ejemplo )

Observamos que las variables no numéricas están reconocidas como variables tipo “factor” o tipo “carácter”. Aunque en muchos casos las funciones de R convierten automáticamente las variables de texto en factores, esta transformación es necesaria en diversos tipos de análisis y funciones, por lo que generalmente interesa hacer una conversión permanente con la función factor. Por ejemplo:

factor( ejemplo$ind ) -> ejemplo$ind

Conviene siempre hacer inicialmente una representación gráfica de los datos:

plot( ejemplo$cov, ejemplo$y, col = ejemplo$trat )

Para evitar tener que referirnos constantemente al objeto ejemplo al nombrar sus variables, podríamos utilizar la función attach, aunque no es una opción muy recomendable. Es mejor utilizar el argumento data, admitido por la mayoría de las funciones:

plot( y ~ cov, data = ejemplo, col = trat, pch = 19 )

La función boxplot resulta muy útil para representar los datos en función de variables cualitativas:

boxplot( y ~ trat, data = ejemplo )

Para analizar el efecto de la variable cualitativa realizaremos un análisis de la varianza (ANOVA) con la función aov (y summary, para mostrar los resultados con más detalle):

summary( aov( y ~ trat, data = ejemplo ) )

El análisis anterior es equivalente a una prueba de la \(t\) de Student:

t.test( y ~ trat, data = ejemplo, var.equal = TRUE )

Como alternativa se pueden utilizar pruebas no paramétricas. Por ejemplo:

wilcox.test( y ~ trat, data = ejemplo )
kruskal.test( y ~ trat, data = ejemplo )

Otra alternativa es un test de aleatorización, que realizaremos con la función MonteCarlo, incluida en aedbr.RData:

MonteCarlo( ejemplo$y, ejemplo$trat, sim = 1000 )

Con la variable independiente cuantitativa haremos una regresión (modelo lineal, con errores de distribución normal o gaussiana):

lm( y ~  cov, data = ejemplo ) -> mod1
summary( mod1 )
plot( y ~ cov, data = ejemplo )
abline( 14.098, 3.089 )

Una representación gráfica adecuada de modelo debe incluir el intervalo de confianza. Para ello utilizaremos la función predictorEffects del paquete effects:

library( effects )
plot( predictorEffects( mod1 ) )

Si queremos probar el ajuste a un modelo unimodal (parábola), solo tenemos que añadir en la ecuación el término cuadrático de la variable cuantitativa (estadísticamente no significativo en este caso):

summary( lm( y ~ cov + I( cov^2 ), data = ejemplo ) )

También podemos usar un modelo lineal para hacer un análisis equivalente al ANOVA:

summary( lm( y ~ trat, data = ejemplo ) )

Los modelos lineales nos permiten además combinar variables de diferente naturaleza. Por ejemplo, un análisis de la covarianza (ANCOVA):

lm( y ~ trat + cov, data = ejemplo ) -> mod2
summary( mod2 )
plot( y ~ cov, col = trat, data = ejemplo, pch = 19 )
abline( 20.9283, 3.0894 )
abline( 21.4583 -13.66, 3.0894, col = 2 )

Interpretemos un modelo más complejo que incluya términos de interacción entre las variables:

lm( y ~ trat * cov, data = ejemplo ) -> mod3
summary( mod3 )
abline( 30.55, 1.34, lwd = 2 )
abline( 30.55 - 32.9033, 1.34 + 3.4988, lwd = 2, col = 2 )

La representación de este último modelo, con los intervalos de confianza, la haríamos con:

plot( predictorEffects( mod3 ) )
plot( predictorEffects( mod3, ~ trat, xlevels = list( cov = c( 2, 5, 9) ) ) )
plot( predictorEffects( mod3, ~ cov ), multiline = TRUE, confint = list( style = "bands" ) )

También podemos utilizar la función aov para hacer un ANCOVA. Sin embargo, los estadísticos que proporciona aov en este caso no son apropiados y debe utilizarse junto a otra función adicional (drop1):

drop1( aov( y ~ trat * cov, data = ejemplo ), ~., test = "F" )

Una función muy útil es anova ( no debe confundirse con aov), que permite comparar dos o más modelos:

anova( mod1, mod2, mod3 )

ANOVA de dos factores

Analizaremos un ejemplo con datos procedentes de un diseño de bloques aleatorios. Se trata de la respuesta de crecimiento de la hierba algodonera (Eriophorum angustifolium) a cuatro tratamientos de fertilización en cinco localidades (bloques) en la tundra de Alaska (fuente: Krebs, 1999, pág. 358).

cottongrass
aov( growth ~ trat + loc, data = cottongrass ) -> modB
summary( modB )

Las comparaciones múltiples las podemos hacer con el test de Tukey:

TukeyHSD( modB )
plot( TukeyHSD( modB, "trat" ) )

El test de Tukey también es aplicable a análisis de ANOVA de un único factor con más de dos tratamientos. Por ejemplo, examinemos el resultado de no incluir el efecto bloque en el análisis:

x11()
plot( TukeyHSD( aov( growth ~ trat, data = cottongrass ) ) )

Modelos lineales generalizados

Utilizaremos una matriz de datos simulados sobre la ocupación y reproducción de una rapaz territorial:

territorios
str( territorios )

En esta tabla aparecen símbolos NA (Not Available) en la variable pollos que representan valores perdidos (missing values), y que se corresponden con los casos en los que no hay ocupación. Con este ejemplo se explica bien la diferencia entre un 0 (la pareja no se reproduce o fracasa en la reproducción) y un NA ( no hay valor de reproducción porque el territorio no está ocupado ).

La ocupación, expresada con datos binarios (1/0), la analizaremos con modelos de regresión logística (distribución binomial de errores), usando la función glm:

glm( ocup ~ alt, data = territorios, family = binomial ) -> modG1
summary( modG1 )

Comprueba cómo se obtiene la curva ajustada utilizando la ecuación del modelo logístico y la función plogis, que realiza la transformación logística inversa:

\[\ln \left( \dfrac{ p }{ 1-p } \right)= b_0 + b_1 x \quad \rightarrow \quad p = \dfrac{ e^{ b_0 + b_1 x } }{ 1 - e^{ b_0 + b_1 x } }\]

plot( ocup ~ alt, data = territorios )
exp( -8.869028 + 0.010527 * ( 700:1100 ) ) -> Ey
lines( 700:1100, Ey / ( 1 + Ey ) )

Para mayor comodidad puede usarse la función plogis, que realiza la transformación logística inversa:

plot( ocup ~ alt, data = territorios )
plogis( -8.869028 + 0.010527 * ( 700:1100 ) ) -> p
lines( 700:1100, p )

La representación con los intervalos de confianza de este último modelo, requiere un nuevo argumento (type = "response"):

plot( predictorEffects( modG1 ), type = "response", ylim = c( 0, 1 ) )

Por su parte, la reproducción (productividad, expresada como número de pollos) la analizaremos con modelos de regresión discreta o de Poisson (datos de conteos con distribución de errores de Poisson):

glm( pollos ~ alt, data = territorios, family = poisson ) -> modG2
summary( modG2 )

Comprueba cómo se obtiene la curva ajustada utilizando la ecuación del modelo de Poisson:

\[\ln{E y} = b_0 + b_1 x \quad \rightarrow \quad E y = e^{ b_0 + b_1 x }\]

plot( pollos ~ alt, data = territorios )
exp( 0.8613412 -0.0008336 * ( 700:1100 ) ) -> Ey
lines( 700:1100, Ey )
plot( predictorEffects( modG2 ), type = "response", ylim = c( 0, 2 ) )

Interpretemos ahora los modelos con variables cualitativas:

glm( pollos ~ orient, data = territorios, family = poisson ) -> modG3
summary( modG3 )

Utilizaremos también la función anova para testar la significación estadística de cada variable/factor:

anova( modG3, test = "Chisq" )

Las diferencias entre los grupos de la variable cualitativa pueden analizarse modificando la modalidad de referencia con la función relevel. Por ejemplo:

glm( pollos ~ relevel( orient, 2 ), data = territorios, family = poisson ) -> modG3b
summary( modG3b )

En estos modelos con variables independientes cualitativas puede resultar interesante aplicar la opción de eliminar la estimación de la constante ( ~ -1 ), lo que facilita la interpretación de los coeficientes de cada modalidad (aunque no puedan evaluarse las diferencias estadísticas entre ellas):

glm( pollos ~ -1 + orient, data = territorios, family = poisson ) -> modG3c
summary( modG3c )

Finalmente, el análisis “completo” de la ocupación y la reproducción lo obtendríamos introduciendo los términos de interacción entre las dos variables ambientales:

glm( ocup ~ alt * orient, data = territorios, family = binomial ) -> modG4
anova( modG4, test = "Chisq" )
plot( predictorEffects( modG4, ~ alt ), type = "response", multiline = TRUE, confint = list( style = "bands" ) )
glm( pollos ~ alt * orient, data = territorios, family = poisson ) -> modG5
anova( modG5, test = "Chisq" )
plot( predictorEffects( modG5, ~ alt ), type = "response", multiline = TRUE, confint = list( style = "bands" ) )

Modelo “clásico” de regresión logística

La regresión logística se aplica habitualmente a datos procedentes de experimentos en los que la respuesta biológica se expresa de manera binaria, y que se presentan pareados como número de “éxitos” y número de “fracasos” (o número total de “intentos”). Usaremos como ejemplo los datos de un experimento de germinación de la planta halófila Halocnemum strobilaceum. Se indica el número de éxitos (semillas que germinan) y fracasos (semillas que no germinan) en réplicas sometidas a diferentes tratamientos de presión osmótica (fuente: Pujol 2000).

halocnemum
glm( cbind( germ, nogerm ) ~ pres, data = halocnemum, family = binomial ) -> modG6
summary( modG6 )
plot( germ/25 ~ pres, data = halocnemum, pch = 19 )
lines( modG6$fit ~ pres, data = halocnemum, lwd = 2 )
plot( predictorEffects( modG6 ), type = "response", ylim = c( 0, 1.1 ) )

Diagnósticos

Una vez realizados nuestros análisis conviene evaluar las asunciones de los modelos (linealidad, homocedasticidad, normalidad), así como detectar la existencia de puntos u observaciones atípicas (ouliers), observaciones palanca (leverage points) y observaciones influyentes. La función plot aplicada al objeto que contiene el modelo proporciona gráficas diagnósticas sobre las asunciones, y la función influence.measures proporciona información (mediante el cálculo de diversas medidas) sobre el efecto de cada observación, señalando las que tienen una influencia desproporcionada en el ajuste del modelo. Para el modelo de Halocnemum, por ejemplo:

plot( modG6 )
influence.measures( modG6 )

En el caso de modelos de regresión lineal (distribución de errores gaussiana) en los que existan ouliers u observaciones influyentes se pueden usar como alternativa modelos de regresión robusta (por ejemplo, función rlm del paquete MASS).

La asunción de independencia de los datos (cada observación es independiente de las demás) solo se asegura con el diseño de un adecuado procedimiento experimental o de muestreo, y no puede ser testado con diagramas diagnósticos. En los ejemplos de análisis utilizados hasta ahora en el seminario venimos asumiendo la independencia de las observaciones, pero esta circunstancia no es válida en los datos de los objetos ejemplo y territorio. En el primer caso se trata de medidas repetidas (no independientes) realizadas sobre cuatro individuos, y en el segundo las observaciones se agrupan en territorios. Esta cuestión se abordará de manera apropiada en la sección dedicada a los modelos mixtos.

Sobredispersión

La sobredispersión ocurre cuando los datos presentan una mayor varianza de lo esperado según un determinado modelo de distribución. Esto implica que (desvianza residual / grados de libertad) > 1. Podemos hacer tests de sobredispersión aplicando la función gof del paquete aods3:

library( aods3 )

Siguiendo con el ejemplo de la germinación de Halocnemum, comprobamos como los dos tests que proporciona la función gof son significativos:

gof( modG6 )

Con modelos logísticos puede usarse la familia quasibinomial, que estima idénticos coeficientes, pero que corrige los valores de probabilidad. En este caso, sin embargo, la interpretación de los resultados no cambia:

glm( cbind( germ, nogerm ) ~ pres, data = halocnemum, family = quasibinomial ) -> modQB
summary( modQB )

El problema de la sobredispersión es más frecuente con datos de Poisson. En estos casos, para regresiones discretas existe también la familia quasipoisson. Usaremos como ejemplo los datos correspondientes a conteos de individuos de curruca cabecinegra (Curruca melanocephala) en 58 itinerarios de censo de 1 km en la Región de Murcia. Analizaremos la abundancia de esta especie en función de la altitud (m s. n. m.):

curruca
glm( count ~ alt, data = curruca, family = poisson ) -> modP
summary( modP )
gof( modP )

En este caso, el modelo quasipoisson modifica sustancialmente el valor de \(P\) del coeficiente correspondiente a la variable altitud:

glm( count ~ alt, data = curruca, family = quasipoisson ) -> modQP
summary( modQP )

No obstante, la mejor alternativa para modelar datos de conteos con sobredispersión es la regresión negativa binomial. Aplicaremos la función glm.nb incluida en el paquete MASS:

library( MASS )
glm.nb( count ~ alt, data = curruca ) -> modNB
summary( modNB )

Con el modelo negativo binomial cambian tanto los coeficientes como el valor de probabilidad. Además, en este caso, el test de sobredispersión resulta no significativo, por lo que podría considerarse el modelo más adecuado para los datos de abundancia de curruca:

gof( modNB )
plot( predictorEffects( modNB ), type = "response" )

Otra alternativa para abordar el problema de la sobredispersión consiste, por ejemplo, en especificar la familia negativa binomial en la función glm, introduciendo el parámetro () (theta) estimado con el modelo de regresión negativa binomial:

summary( glm( count ~ alt, data = curruca, family = negative.binomial( modNB$theta ) ) )

Otro problema que surge cuando se analizan datos de conteos es el exceso de ceros, en casos en los que dicho exceso es generado por un proceso distinto del que genera los valores del conteo. Para estos datos existen diversos tipos de modelos de regresión, como zero-inflated, hurdle y zero-truncated, con variantes para las familias Poisson y negativa binomial. Los paquetes pscl y VGAM disponen de funciones que permiten realizar estos análisis.

Puede encontrarse más información sobre estos y otros tipos de regresión (ordinal, multinomial, probit …) en: https://stats.oarc.ucla.edu/other/dae/

Modelos GAM

Los GAM (generalized additive models) permiten el ajuste de los datos a curvas complejas. Como ejemplo de aplicación usaremos de nuevo los datos de la curruca. Emplearemos la función gam del paquete mgcv, especificando la familia y, opcionalmente, un número de nodos para la curva. Por ejemplo:

library( mgcv )
gam( count ~ s( alt, k = 5 ), data = curruca, family = poisson ) -> modGP
summary( modGP )

La curva estimada puede representarse sobre los datos originales:

plot( count ~ alt, data = curruca, pch = 19 )
lines( modGP$fit[ order( curruca$alt ) ] ~ sort( curruca$alt ), lwd = 2, col = 2 )

Esta curva puede compararse con las de los modelos obtenidos en la sección anterior:

lines( modP$fit[ order( curruca$alt ) ] ~ sort( curruca$alt ), lwd = 2, col = 3 )
lines( modNB$fit[ order( curruca$alt ) ] ~ sort( curruca$alt ), lwd = 2, col = 4 )

El modelo GAM sugiere un mejor ajuste de los datos a una curva unimodal, que podría haberse contemplado en el ejercicio de la sección anterior usando como fórmula la expresión alt + I( alt^2 ).

Por último, ejecutaremos un GAM con la familia negativa binomial:

gam( count ~ s( alt, k = 5 ), data = curruca, family = nb ) -> modGNB
lines( modGNB$fit[ order( curruca$alt ) ] ~ sort( curruca$alt ), lwd = 2, col = 5 )

La comparación entre los diferentes modelos, con objeto de determinar cuál proporciona un mejor ajuste a los datos, puede realizarse con la función model.sel del paquete MuMIn (sección “Selección de modelos e inferencia multimodelo”).

Los smoothing splines son una alternativa a los modelos GAM:

smooth.spline( curruca$alt, curruca$count, nk = 5 ) -> modS
lines( modS, lwd = 2, col = 6 )

Regresión no lineal

El ajuste de datos a ecuaciones no lineales puede realizarse mediante la función nls. Como ejemplo, ajustaremos la ecuación de Michaelis-Menten a los datos del objeto puromycin, correspondiente a un estudio sobre la velocidad de una reacción enzimática en células tratadas con puromicina (fuente: paquete de R nlstools):

puromycin
plot( rate ~ conc, data = puromycin )
nls( rate ~ Vmax * conc / ( Km + conc ), data = puromycin, start = c( Vmax = 200, Km  = 0.05 ) ) -> modN
summary( modN )
lines( puromycin$conc, 212.7 * puromycin$conc / ( 0.06412 + puromycin$conc ) )

Modelos mixtos

Los modelos de regresión mixtos permiten analizar de manera apropiada la falta de independencia en datos que se agrupan en función de individuos o unidades experimentales o de muestreo y que corresponden habitualmente a medidas u observaciones repetidas. En un modelo mixto se asume que elemento de agrupación, denominado factor (o efecto) aleatorio es una muestra aleatoria de una población mayor. Por otra parte, las variables potencialmente explicativas de la respuesta biológica objeto de la investigación se denominan factores (o efectos) fijos.

Retomaremos los datos del objeto ejemplo, analizados al inicio del seminario obviando la falta de independencia, y que serán analizados aquí con modelos más apropiados a su estructura. En esta matriz variable ind representa el factor aleatorio y las variables trat y cov se consideran efectos fijos.

ejemplo

Necesitaremos cargar los paquetes nlme y lme4 (este último para modelos mixtos generalizados).

library( nlme )
library( lme4 )

Visualicemos inicialmente la estructura agrupada de los datos en una gráfica:

plot( groupedData( y ~ cov | ind, data = ejemplo ) )

Usaremos la función lme para ejecutar un modelo mixto sencillo:

lme( y ~ cov, random = ~ 1 | ind, data = ejemplo ) -> modM1
summary( modM1 )
modM1$fit

El análisis proporciona una única respuesta al efecto fijo para los cuatro individuos, y otra distinta para cada uno de ellos (efecto aleatorio), aunque en este caso todas con la misma pendiente:

plot( modM1$fit[,2] ~ cov, data = ejemplo, col = ind )
points( modM1$fit[,1] ~ cov, data = ejemplo, pch = 19 )

Un modelo más complejo podría contemplar la correlación entre constante y pendiente:

lme( y ~ cov, random = ~ 1 + cov | ind, data = ejemplo ) -> modM2
summary( modM2 )
plot( modM2$fit[,2] ~ cov, data = ejemplo, col = ind )
points( modM2$fit[,1] ~ cov, data = ejemplo, pch = 19 )

Con la función anova podemos comprobar que el segundo modelo mejora considerablemente el ajuste:

anova( modM1, modM2 )

Para ejecutar modelos lineales mixtos generalizados usaremos las funciones glmer y glmer.nb del paquete lme4. Usaremos de nuevo como ejemplo los datos del objeto territorios, considerando la variable terr como factor aleatorio (1 | terr) y las variables ambientales (alt y orient) como factores fijos. Observa que con glmer la forma de especificar los factores aleatorios es diferente a la de lme:

territorios
glmer( ocup ~ alt + ( 1 | terr ), data = territorios, family = binomial ) -> modM3
summary( modM3 )

En este caso la varianza entre territorios es 0, por lo que el modelo mixto apenas difiere del “clásico” (modG1, sección “Modelos lineales generalizados”):

summary( modG1 )
anova( modM3, modG1 )

Aunque la función anova permite comparar modelos obtenidos con la función glmer (entre sí y con otros tipos de modelos), cuando se aplica sobre un único modelo no proporciona los valores de \(P\) de los coeficientes. En su lugar, puede usarse la función Anova del paquete car.

library( car )
Anova( modM3 )

La función Anova es necesaria también para obtener los valores de \(P\) cuando se utiliza la función lmer del paquete lme4 (equivalente a la función lme del paquete nlme para ajustar modelos lineales mixtos).

Para ejecutar modelos lineales mixtos generalizados para la familia negativa binomial usaremos glmer.nb, en este caso con la variable de reproducción. Comprobaremos que en este caso tampoco existen grandes diferencias con el modelo “clásico”:

glmer.nb( pollos ~ orient + ( 1 | terr ), data = territorios ) -> modM4
summary( modM4 )
summary( modG3 )
anova( modM4, modG3 )

Modelos anidados

Son un tipo particular de modelos mixtos, para datos con una estructura jerárquica. Utilizaremos como ejemplo los datos del objeto glycogen, consistentes en medidas de glucógeno en plasma de ratas sometidas a diferentes tratamientos: dos medidas independientes en 3 preparaciones por rata (6 ratas) y 3 tratamientos (fuente: Sokal y Rohlf 1995, pág. 289). Se trata de un diseño de ANOVA anidado ( Modelo I: tratamiento fijo y niveles subordinados aleatorios ). Es importante señalar que las variables deben ser de tipo factor.

glycogen
aov( gly ~ trat + Error( rat / prep ), data = glycogen ) -> modA
summary( modA )

En el ejemplo anterior, el único test de interés es el del tratamiento, que es el que proporciona la función. No obstante, si queremos obtener los resultados de los otros dos, podemos calcular los valores de \(F\) dividiendo la varianza (Mean Sq) de cada factor por la de su subordinado:

  1. Ratas dentro de tratamientos (\(F\) y \(P\)):

    265.9 / 49.5 
1 - pf( 265.9 / 49.5, 3, 12 )

  2. Preparaciones dentro de ratas (\(F\) y \(P\)):

    49.5 / 21.17
1 - pf( 49.5 / 21.17, 12, 18)

Cuando existen efectos aleatorios interesa conocer los componentes de varianza. Se pueden calcular con los valores del ANOVA anterior, pero es más sencillo utilizar un modelo de regresión mixto, que sirve además como alternativa a la función aov, proporcionando los mismos resultados. Hay que indicar la estructura anidada con el argumento random = ~ 1 | rat/prep:

lme( gly ~ trat, random = ~ 1 | rat/prep, data = glycogen ) -> modM5
anova( modM5 )
VarCorr( modM5 )

En este caso, para expresar los componentes de varianza como porcentajes sobre el total de varianza:

c( 36.06482, 14.16667, 21.16667 ) * 100 / sum( 36.06482, 14.16667, 21.16667 )

Análisis del ejemplo del objeto cottongrass (diseño de bloques aleatorios) con un modelo mixto:

lme( growth ~ trat, random = ~ 1 | loc, data = cottongrass ) -> modM6
anova( modM6 )

Comprobemos que, para el factor tratamiento, proporciona el mismo resultado que el ANOVA realizado anteriormente:

anova( modB )

Análisis split-plot

El diseño split-plot es un caso especial de diseño anidado. Analizaremos como ejemplo los datos de cosecha de tres variedades de avena (expresadas como 1/4 lb por sub-plot, cada uno de 1/80 acre. Hay seis bloques (I-VI) y cuatro tratamientos de fertilización con nitrógeno. El diseño split-plot en este caso presenta una estructura de plots (tres parcelas de cultivo, una con cada variedad) anidados en bloques (distintas fincas), y sub-plots (los cuatro tratamientos de fertilización) anidados en plots (fuente: Venables y Ripley 2002, pág. 282):

avena
lme( yield ~ variety * factor( nitro ), random = ~ 1 | block/variety, data = avena ) -> modM6
anova( modM6 )

El mismo análisis se puede realizar con la función aov:

summary( aov( yield ~ variety * factor( nitro ) + Error( block/variety ), data = avena ) )

Análisis de datos pareados

Es un caso particular de modelo mixto en el que los datos se emparejan por observaciones. Utilizaremos como ejemplo los datos del objeto sewage, que corresponden de un estudio en el que se analiza la densidad (transformada logarítmicamente) de coliformes por ml en aguas sometidas a dos métodos de cloración a lo largo de 8 días de tratamiento (fuente: paquete de R PairedData). El método clásico de análisis en estos casos es un test de la \(t\) (especificando el argumento paired = TRUE), pero también pueden analizarse con un modelo mixto, con el mismo resultado:

sewage
t.test( coliform ~ method, data = sewage, paired = TRUE )
lme( coliform ~ method, random = ~ 1 | day, data = sewage ) -> modM7
summary( modM7 )

Modelos multinomiales

Los modelos multinomiales son una generalización de los modelos binomiales para analizar variables de respuesta con más de dos modalidades. Necesitaremos cargar el paquete nnet para usar la función multinom.

library( nnet )

Como ejemplo usaremos una tabla de datos (multiestado) derivada de la matriz del objeto territorios. En este caso la variable de respuesta (estado) refleja tres modalidades para cada territorio: vacío (0), ocupado sin éxito (1) y ocupado con éxito (2), que relacionaremos con la altitud. Obviaremos también aquí la estructura agrupada en territorios.

multiestado
multinom( estado ~ alt, data = multiestado ) -> modM8
summary( modM8 )

Como la función multinom no proporciona valores de \(P\), necesitaremos comparar con la función anova el modelo de interés con el modelo nulo (éxito ~ 1), o bien usar directamente la función Anova del paquete car:

multinom( estado ~ 1, data = multiestado ) -> modM0
anova( modM8, modM0 )
Anova( modM8 )

La interpretación de los resultados de un modelo multinomial a través de sus coeficientes es compleja, pero podemos recurrir a su representación gráfica:

plot( predictorEffects( modM8 ), multiline = TRUE, confint = list( style = "bands" ) )
x11()
plot( predictorEffects( modM8 ), axes = list( y = list( style = "stacked" ) ) )

Tablas de contingencia: chi-cuadrado, regresión logística y análisis loglineal

El análisis de datos categóricos, frecuentemente organizados en forma de tabla de contingencia, presenta una elevada complejidad, con numerosas alternativas de análisis que responden a diferentes enfoques y estructuras de los datos. Como ejemplo en esta sección usaremos una tabla de datos (reprod) derivada de la matriz del objeto territorios, con objeto de analizar el éxito en la reproducción (pollos > 0) en función de la orientación. Obviaremos también en este caso la estructura agrupada en territorios asumiendo independencia en las observaciones.

Comenzaremos elaborando una tabla de contingencia a partir de la matriz de datos:

reprod
table( reprod ) -> reprodT
reprodT

Las pruebas clásicas de análisis son el test de \(\chi^2\) y el test exacto de Fisher:

chisq.test( reprodT )
fisher.test( reprodT )

La alternativa es una regresión logística con la tabla original (reprod):

glm( éxito ~ orientación, data = reprod, family = binomial ) -> modT1
Anova( modT1 )

También una regresión logística con la tabla modificada (estructura “éxitos/fracasos”):

glm( reprodT ~ rownames( reprodT ), family = binomial ) -> modT2
Anova( modT2 )

Modificaremos ahora el formato de la tabla para expresar los datos como frecuencias:

data.frame( reprodT ) -> reprodF
reprodF

De esta forma podemos usar una regresión logística usando la opción de ponderación (weights = Freq):

glm( éxito ~ orientación, weights = Freq, data = reprodF, family = binomial ) -> modT3
Anova( modT3 )

Otra alternativa de análisis son los modelos loglineales: se utiliza una regresión de Poisson con la tabla de frecuencias (donde Freq es la variable dependiente) para obtener el denominado “modelo saturado”. El único test de interés aquí es el correspondiente al término de interacción:

glm( Freq ~ éxito * orientación, data = reprodF, family = poisson ) -> modT4
Anova( modT4 )

Cálculo de frecuencias estimadas de cada celda con modelos loglineales:

El modelo saturado (modT4) predice la frecuencia de cada celda de la tabla reprodT:

summary( modT4 )
predict( modT4, type = "response" )

En cualquier caso, utilizando los coeficientes podemos calcular las frecuencias estimadas. Para el caso del modelo modT4, por ejemplo, la frecuencia de éxitos en orientación Este se obtiene con:

exp( 1.9459 -0.3365 )

La frecuencia de fracasos en orientación Sur:

exp( 1.9459 + 0.1335 )

Y la frecuencia de éxitos en orientación Norte:

exp( 1.9459 -0.3365 -24.2485 + 25.5835 )

Por último, analizaremos los datos de la tabla original con un modelo multinomial:

multinom( éxito ~ orientación, data = reprod ) -> modT5
summary( modT5 )
Anova( modT5 )

Tablas de contingencia con más de un factor

Usaremos los datos del objeto aliens, que presenta, en forma de tabla de frecuencias, información sobre 744 especies invasoras de ecosistemas de agua dulce europeos, caracterizadas en función de variables cualitativas (fuente: Nunes et al. 2015). Analizaremos si el carácter de “alto impacto” (variable hi) de las especies está relacionado con la vía de invasión (pathway) y con el reino al que pertenecen (kingdom).

aliens
xtabs( freq ~ hi + pathway + kingdom, data = aliens )

Realizaremos un análisis loglineal con las frecuencias como variable dependiente:

glm( freq ~ hi * pathway * kingdom, data = aliens, family = poisson ) -> modT6
Anova( modT6)

Para la interpretación de los resultados solo interesan los tests correspondientes a los términos de interacción con la variable dependiente (hi:pathway, hi:kingdom y hi:pathway:kingdom).

Selección de modelos e inferencia multimodelo

Estos procedimientos constituyen un nuevo paradigma para el análisis estadístico. En síntesis, un procedimiento de selección de modelos evalúa el ajuste a los datos de un conjunto de modelos “candidatos”, y realiza la estimación de parámetros en función de las estimaciones ponderadas de dicho parámetro en los diferentes modelos. Esta evaluación se realiza generalmente mediante el criterio de información de Akaike (AIC), o su versión corregida (AICc). Se busca la “mejor” explicación posible para los datos observados (best approximating model) y, bajo esta perspectiva, el uso de los valores de \(P\) de la estadística frecuentista no tiene sentido.

Para el ejemplo de esta sección utilizaremos datos de concentraciones de calcio (mg/100 ml) en plasma sanguíneo de aves de ambos sexos sometidas a tres tratamientos hormonales (fuente: Zar 2014, pág. 265). Necesitaremos cargar el paquete MuMIn.

library( MuMIn )
calcio

El procedimiento de selección puede realizarse con cualquier tipo de modelo. Usaremos aquí la función glm para realizar regresiones lineales. Si no se especifica, la función elige por defecto la familia gaussiana, y los resultados son equivalentes a los que se obtendrían con la función lm.

Utilizaremos cinco modelos candidatos distintos, incluido el modelo nulo:

glm( Ca ~ 1, data = calcio ) -> mc0
glm( Ca ~ hormone, data = calcio ) -> mc1
glm( Ca ~ sex, data = calcio ) -> mc2
glm( Ca ~ hormone + sex, data = calcio ) -> mc3
glm( Ca ~ hormone * sex, data = calcio ) -> mc4

Con el conjunto de modelos “candidatos” obtendremos la tabla:

model.sel( mc0, mc1, mc2, mc3, mc4 ) -> ms
ms
coef( ms )

Con la función de inferencia multimodelo (model.avg) estimaremos los diferentes coeficientes:

summary( model.avg( ms ) )

La función proporciona dos tablas de coeficientes, una en la que se consideran todos los modelos candidatos (full average), y otra en la que se consideran solo los modelos donde aparece la variable en cuestión (conditional average). También resulta muy interesante la información sobre la “importancia relativa” de las variables (suma de los pesos de Akaike de los modelos en los que aparece cada variable).

sw( model.avg( ms ) )

Los valores estimados por el modelo promediado se pueden obtener con la función predict:

predict( model.avg( ms ), se.fit = TRUE )

Referencias

  • Krebs CJ ( 1999 ) Ecological Methodology. 2ª ed. Benjamin/Cummings, Menlo Park, CA.

  • Nunes AL, Tricarico E, Panov VE, Cardoso AC, Katsanevakis S ( 2015 ) Pathways and gateways of freshwater invasions in Europe. Aquatic Invasions 10: 359-370.

  • Pujol JA ( 2000 ) Ecología de la germinación de semillas y del crecimiento de plántulas de especies halófitas del sureste ibérico. Tesis Doctoral. Universidad de Murcia.

  • Sokal RR, Rohlf FJ ( 1995 ) Biometry. 3ª ed. Freeman, New York.

  • Venables WN, Ripley BD ( 2002 ) Modern Applied Statistics with S. 4ª ed. Springer, New York.

  • Zar JH ( 2014 ) Biostatistical Analysis. 5ª ed. Pearson, Harlow, UK.

Descripción de los datos y funciones

Utiliza la función info. Por ejemplo: info( "aliens" ) o info( "MonteCarlo" ).