Tutorial (ecología.funcional_PGLS): Mínimos Cuadrados Generalizados Filogenéticos (PGLS) en R

Introducción

En este tutorial aprenderemos a implementar modelos PGLS (“Phylogenetic Generalized Least Squares”, o Mínimos Cuadrados Generalizados Filogenéticos) utilizando PGLS en con diferentes aproximaciones (serie de paquetes) en R. Pero antes es bueno entender qué es un PGLS, desde el punto de vista del modelado estadístico (ver Cuadro 1). El PGLS es, básicamente, un modelo de regresión lineal mixto, generalizado (GLMM) o no, específicado para controlar la falta de independencia entre especies - que tienen una historia evolutiva en común - por medio de una matriz de varianza-covarianza construída a partir de uno o más árboles filogenéticos.

Cuadro 1: Cuando un modelo de regresión lineal, simple (SLR) o múltiple (MLR), o su análogo generalizado (GLM), incluye uno o más términos aleatorios para dar cuenta de la existencia de correlaciones entre los datos en el predictor lineal, pasa a ser un modelo mixto (LMM: “Linear Mixed Model”, GLMM: “Generalized Linear Mixed Model”, o Modelo Lineal Mixto y Modelo Lineal Generalizado Mixto, respectivamente). En Wikipedia podemos leer que para la estimación de los parámetros de los modelos mixtos los métodos de Mínimos Cuadrados Ordinarios (OLS, de “Ordinary Least Squares”) suelen ser estadísticamente ineficaces. Entonces los LMM y GLMM hacen uso de estimación por Mínimos Cuadrados Generalizados (“GLS”, por la sigla en inglés de “Generalized Least Squares”), una técnica para estimar los parámetros desconocidos de un modelo de regresión lineal cuando existe cierto grado de correlación entre los residuos de un modelo de regresión. Los GLS pueden utilizarse para estimar parámetros en modelos con una estructura jerárquica, o anidamiento (como los llamados modelos multi-nivel), ya que lidian con problemas de heterogeneidad, anidamiento, correlación temporal y espacial - y otros tipos de correlación - (ver capítulo 4 de Zuur et al. 2009: “Dealing with heterogeneity”).

Para introducir los PGLS me basaré casi textualmente en la excelente revisión de Symonds y Blomberg (2014) “A Primer on Phylogenetic Generalised Least Squares” (capítulo 5 del libro “Modern Phylogenetic Comparative Methods and Their Application in Evolutionary Biology”, editado por László Zsolt Garamszegi, ver haciendo click). Los autores resaltan que los mínimos cuadrados generalizados filogenéticos (PGLS) son uno de los principales métodos estadísticos empleados para establecer la naturaleza de la asociación evolutiva entre dos o más rasgos biológicos (por ejemplo, la relación entre la masa corporal y la esperanza de vida) teniendo en cuenta la filogenia de las especies consideradas. Por “asociación evolutiva” entendemos la evidencia de que los rasgos están asociados a lo largo del tiempo evolutivo.

Cuadro 2: Otros capítulos del libro son de lectura recomendable. Los Capítulos 2 y 3 proporcionan una mayor discusión sobre la preparación de filogenias para el análisis comparativo. incluída inferencia basada en árboles compuestos y superárboles. En el Capítulo 6 se revisan las consideraciones estadísticas sobre los supuestos asumidos por los PGLS y sus limitaciones en tópicos con el tamaño muestral. En el Capítulo 7 se cubren varios métodos para tratar la variación intraespecífica y el error de medición en el marco PGLS. En el capítulo 9 se revisan enfoques que afrontan la comprobación de hipótesis con PGLS para análisis con una variable de respuesta de carácter discreto (como la prueba de cambios concentrados, las comparaciones por pares, el método de verosimilitud de Pagel y la regresión logística filogenética. El Capítulo 12 aplica inferencia multimodelo para combinar muchas posibles resoluciones del árbol filogenético. Además el Capítulo 10 y el 12 cubren aproximaciones Bayesianas para PGLS.

La forma más sencilla de concebir la PGLS es como una regresión ponderada enfocada a estudiar la asociación entre dos o más rasgos (sopesando las relaciones filogenéticas entre las especies). En una regresión estándar (ajustada mediante “OLS”, el método de Mínimos Cuadrados Ordinarios), cada punto de datos independiente contribuye por igual a la estimación de la línea de regresión. En cambio, la PGLS “pondera a la baja” los puntos que proceden de especies con una historia filogenética común (usando GLS). Lo distintivo es que un PGLS permite incorporar en el modelo de regresión el hecho de que las especies estrechamente emparentadas tienen rasgos más parecidos debido a su ascendencia compartida y, por tanto, producen residuos más parecidos a partir de la línea de regresión de mínimos cuadrados. En otras palabras, los PGLS incorporan la estructura de covarianza esperada de estos residuos (calculada a partir del árbol filogenético), y generan estimaciones modificadas de la pendiente y el intercepto que pueden dar cuenta de la autocorrelación interespecífica debida a la filogenia.

Para los mínimos cuadrados generalizados (GLS) necesitaremos considerar entonces un elemento adicional de la ecuación de regresión, que es la matriz de varianza-covarianza que representa la estructura de covarianza esperada de los residuos de la ecuación de regresión. En el caso de OLS, el supuesto implícito es que no hay covarianza entre los residuos (es decir, todas las especies son independientes entre sí, y los residuos de especies estrechamente relacionadas no son más similares en promedio que los residuos de especies distantes). En otras palabras, sigue el supuesto de que no hay efectos filogenéticos en los residuos. Los elementos diagonales de la matriz de varianza-covarianza (la línea de valores desde arriba a la izquierda hasta abajo a la derecha) representan la varianza de los residuos, mientras que los restantes elementos no diagonales son iguales a cero, lo que significa que no hay covariación entre los residuos. Cuando se asume esta estructura de varianza-covarianza, los resultados de GLS son los mismos que los de OLS (la contribución de la matriz al cálculo de la regresión es casi nula).

Usualmente los puntos de datos de las especies no son independientes debido a su historia filogenética común. En consecuencia, los errores también pueden no ser independientes o estar autocorrelacionados (los residuos de especies estrechamente relacionadas pueden ser similares). Por lo tanto, habrá covariación en los residuos, que debemos tener en cuenta en nuestra matriz de varianza-covarianza (Los elementos no diagonales podrán ser distintos a cero). La cantidad de historia evolutiva compartida entre las especies entonces estará representada por esa covarianza esperada para los residuos de la regresión. Los elementos diagonales, que mencionamos que representan la varianza de los residuos, tomarán el valor de la longitud total de las ramas desde la raíz del árbol hasta las puntas. Esto será igual para cada celda de la diagonal si la filogenia es ultramétrica (es decir, todas las puntas están a la misma distancia de la raíz de la filogenia).

La estructura de covarianza esperada de los residuos es calculada a partir del árbol filogenético suponiendo un modelo evolutivo para los rasgos analizados. En su formulación más básica el PGLS supone que los rasgos continuos evolucionan según un proceso de “paseo aleatorio” (o “movimiento browniano”) de tal forma que el cambio en el valor de un rasgo en un periodo de tiempo dado viene dado por un número aleatorio extraído de una distribución normal con una desviación estándar dada y una media de 0 (el valor tiene la misma probabilidad de aumentar o disminuir). Según este modelo, las especies que comparten un antepasado común más reciente deberían tener (por eso hablamos de covarianza esperada) valores de rasgos más parecidos que las especies emparentadas a mayor distancia, ya que sus rasgos han tenido menos tiempo para divergir. Esto refleja que tengan valores similares para los residuos en la matriz de varianza-covarianza. Sin embargo, hay muchas situaciones en las que los rasgos son evolutivamente lábiles, en las que las especies estrechamente emparentadas no son necesariamente más parecidas en sus rasgos. Entonces, una cuestión central es si la matriz de varianza-covarianza filogenética esperada describe con exactitud la estructura de error de los datos (estructura de error observada).

Suponemos que la filogenia es exacta y que los valores de los rasgos de las especies han evolucionado mediante un modelo de movimiento browniano de evolución gradual, en el que la cantidad de cambio evolutivo a lo largo de una rama es proporcional a la longitud de la rama. Sin embargo, si la filogenia o el modelo evolutivo no son exactos y, en realidad, hay menos o ninguna covarianza filogenética en los residuos (la expectativa OLS), entonces utilizar la filogenia tal y como está estimada puede ser inadecuado. Lo que necesitamos es una forma de determinar el grado de autocorrelación filogenética en los datos. Esto se puede conseguir estimando la señal filogenética, una medida de la bondad del ajuste de la matriz de varianza-covarianza filogenética esperada respecto de la estructura de la covarianza observada. La estructura de covarianza observada es calculada a partir de los residuos del modelo como la covarianza entre un par de especies \(i,j\) para cada uno de sus rasgos \(X\) (esto es, el desvío del valor del rasgo para cada una de las especies del par respecto del valor medio del rasgo (\(X-\bar{X}\)).

La magnitud de la señal nos guiará en la decisión de corregir, o no, por la filogenia. Una medida de la señal filogenética es \(\lambda\), que es un multiplicador de los elementos no diagonales de la matriz de varianza-covarianza esperada, estimado mediante máxima verosimilitud (no por el PGLS). Un valor de \(\lambda = 0\) es coherente con la ausencia de señal filogenética en el rasgo (lleva a cero los elementos no diagonales, produciendo la matriz de covarianza no filogenética). Si no hay señal filogenética en los datos, el PGLS devolverá estimaciones idénticas a las de un análisis de regresión por OLS. Mientras que, un valor de \(\lambda = 1\) es coherente con una fuerte señal filogenética (da una matriz idéntica a la matriz de covarianza filogenética esperada bajo un modelo de evolución de movimiento browniano). Los valores intermedios de \(\lambda\) indican una señal filogenética intermedia. Para \(\lambda < 1\), se acortan las ramas internas y se alargan las ramas terminales del árbol. Los valores superiores a 1 no son válidos porque los valores no diagonales de la matriz de covarianza no pueden superar a los diagonales en GLS (las especies no pueden ser más parecidas a otras especies que a sí mismas). Si la señal filogenética es intermedia, el PGLS puede corregir la filogenia en el grado adecuado.

Es importante recordar que la señal filogenética, en el contexto global del PGLS, se cuantifica en torno a los errores residuales del modelo de regresión (no a la intensidad de la señal en la variable de respuesta o en las variables predictoras). Sin embargo, muchos de los paquetes de R que ajustan PGLS permiten estimar \(\lambda\) para rasgos individuales, además del global. En consecuencia, los valores de \(\lambda\) para la regresión PGLS pueden variar con respecto a los propios rasgos individuales. Por ejemplo, la estimación de máxima verosimilitud de \(\lambda\) para la regresión puede ser 0, y los valores de los rasgos individuales ser de 1 y 0.9 (para la variable respuesta y la predictora, respectivamente). Así, aunque hay una fuerte señal filogenética en nuestros rasgos individuales, no hay señal cuando ambos rasgos se relacionan en un modelo, y las estimaciones reales de la regresión filogenética serán idénticas a las estimaciones de la regresión por OLS. Tengan en cuenta que esto sólo se aplica si su filogenia es ultramétrica (todas las puntas están a la misma distancia de la raíz del árbol).

Formalmente no existe un criterio definido para interpretar si los valores intermedios de \(\lambda\) indican una señal filogenética “débil” o “fuerte”, ya que depende del perfil de verosimilitud de \(\lambda\) para el conjunto de datos específico. Sin embargo, se pueden utilizar pruebas de razón de verosimilitudes (\(LR\)) y calcular valores \(P\) para evaluar si el valor estimado de máxima verosimilitud de \(\lambda\) difiere significativamente de 0 ó 1. Hay quienes sostienen, sin embargo, que dichas pruebas de cociente de verosimilitud en las que el valor nulo no puede superar un determinado valor (como menos de 0 o más de 1) son inherentemente conservadoras. Además \(\lambda\) es sensible al número de especies, y puede funcionar mal como una medida de la señal filogenética en tamaños de muestra pequeños.

Existen otras dos medidas de la señal filogenética, relacionadas con \(\lambda\), que también son modificadores de la longitud de rama y se calculan mediante estimación de máxima verosimilitud. La primera de ellas, \(\delta\), es una transformación de potencia de las longitudes de rama sumadas desde la raíz hasta las puntas del árbol, y la segunda, \(\kappa\), es una transformación de potencia de las propias longitudes de rama individuales. Al igual que con \(\lambda\), ambas pueden utilizarse para inferir algo sobre el proceso evolutivo. \(\delta\) es una medida de si la evolución de los rasgos se ha acelerado (\(\delta > 1\)) o ralentizado (\(\delta < 1\)) a lo largo del tiempo evolutivo. \(\kappa\) es una medida del modo de evolución, con \(\kappa = 0\) representando un cambio evolutivo que es independiente de la longitud de rama, lo que indica un modelo puntuado de evolución. Al igual que \(\lambda\), tanto \(\delta\) como \(\kappa\) pueden aplicarse al cálculo de PGLS, aunque no se utilizan tan a menudo como \(\lambda\) en ese contexto.

Aunque el PGLS se utiliza con frecuencia para examinar la asociación entre un par de rasgos, también puede manejar múltiples variables predictoras. Para la PGLS, la variable dependiente (respuesta) suele ser una variable continua. La(s) variable(s) predictora(s) también puede(n) ser continua(s), pero el PGLS puede tratar datos ordinales pseudocontinuos y datos discretos binarios. Las variables discretas multiestado con propiedades no ordinales (por ejemplo, dieta: insectívora, herbívora, piscívora, etc.) pueden tratarse en un marco PGLS si se recodifican como caracteres binarios separados (por ejemplo, piscívora: no (0) o sí (1)). La comprobación de hipótesis con PGLS no es apropiada para análisis con una variable de respuesta de carácter discreto. Existen métodos distintos para tratar variables de respuesta discretas, como la prueba de cambios concentrados, las comparaciones por pares, el método de verosimilitud de Pagel y la regresión logística filogenética. En el capítulo 9 del libro se revisan algunos de estos enfoques.

Cuadro 3: Los dos requisitos para un análisis PGLS son un conjunto de datos de rasgos comparativos de especies y una filogenia para esas especies. La filogenia puede producirse de novo a partir del análisis filogenético de datos de secuencias de ADN, por ejemplo. También puede tomarse de una fuente ya publicada y podarse hasta las especies pertinentes, o aumentarse como una filogenia compuesta utilizando otras fuentes. Lo ideal es que la filogenia incluya la longitud de las ramas y esté totalmente resuelta. Si no está totalmente resuelta, hay que determinar si las politomías (cuando más de dos especies descienden de un nodo) representan una filogenia conocida o desconocida (es decir, el verdadero proceso evolutivo, en cuyo caso las llamamos politomías “duras”) o sólo incertidumbre sobre el verdadero patrón de relaciones (politomías “blandas”). Puede ocurrir que no se disponga de información sobre la longitud de las ramas de la filogenia, en cuyo caso se pueden establecer todas las longitudes de las ramas como iguales, o utilizar un algoritmo en el que la profundidad de cada nodo del árbol está relacionada con el número de especies hijas derivadas de ese nodo, entre otras estrategias. Una vez compilada, la filogenia debe formatearse para que pueda ser leída por el paquete informático utilizado para el análisis. Normalmente, será un archivo Nexus con el árbol almacenado presentado en ese archivo en formato Newick. La mayoría de los paquetes de análisis filogenético y manipulación de árboles pueden guardar los árboles en este formato. Varios paquetes dentro del marco estadístico R pueden derivar estimaciones PGLS muy rápida y eficientemente, incluyendo “ape” (Paradis et al. 2004), “picante” (Kembel et al. 2010), “caper” (Orme et al. 2012), “phytools” (Revell 2012), “nlme” (Pinheiro et al. 2013), y “phyreg” (Grafen 2014).

Ejemplos de relaciones lineales

En esta sección mostraré cómo implementar PGLS con diferentes paquetes de R cuando la variable dependiente varía linealmente en función de una o más variables independientes.

PGLS con gls de nlme, ape y phytools

El primer ejemplo proviene del material de un taller impartido en Ilhabela, Brasil, en 2015 Comparative methods in R, ver haciendo click. Hace uso del paquete “ape” (“Analyses of Phylogenetics and Evolution”) y del paquete “nlme” (“Nonlinear Mixed-Effects Models”, el clásico de modelos mixtos de José Pinheiro y Douglas Bates), entre otros.

rm(list=ls())  # limpia variables de la memoria.
graphics.off() # limpia gráficos.

library(ape)
# install.packages("geiger") # ejecuta instalación (si no se tiene ya) 
library(geiger)
library(nlme)
library(phytools)

En el ejemplo usan datos de los rasgos de unos lagartos del género \(Anolis\) y un árbol filogenético con sus especies (en formato “.phy”). Las funciones read.table(“myDataFile.csv”, header=TRUE, sep=“,”) o “read.csv” permiten leer un archivo delimitado por comas (csv) con el encabezados de columna. La función “read.tree” lee un árbol filogfenético con formato “Newick”. Todos estos archivos deben estar en el directorio de trabajo del/la usuario/a.

anoleData<-read.csv("./0Data/PGLS_example/anolisDataAppended.csv", row.names=1)
# Preparación de datos e Introducción del Problema:

# Visualizo estructura y matriz de datos de los rasgos de las especies de Anolis
str(anoleData)  # clase de objeto, tipo y número de observaciones de las variables
#> 'data.frame':    100 obs. of  10 variables:
#>  $ SVL                   : num  4.04 3.82 3.53 4.04 4.38 ...
#>  $ PCI_limbs             : num  -3.248 3.409 2.582 0.123 2.036 ...
#>  $ PCII_head             : num  0.372 -1.783 0.779 -1.675 -3.743 ...
#>  $ PCIII_padwidth_vs_tail: num  -1.04 2.21 -2.45 2.39 0.52 ...
#>  $ PCIV_lamella_num      : num  -2.415 0.95 0.7 0.585 2.056 ...
#>  $ awesomeness           : num  -0.24165 -0.25903 -0.46206 0.00687 -0.13231 ...
#>  $ hostility             : num  -0.173 0.127 0.178 0.145 0.465 ...
#>  $ attitude              : num  0.6444 0.296 -0.0196 0.2938 -0.5496 ...
#>  $ ecomorph              : chr  "TG" "TW" "GB" "TC" ...
#>  $ island                : chr  "Cuba" "Cuba" "Cuba" "Hispaniola" ...
# View(anoleData) # visualiza la matriz de datos
head(anoleData) # alternativa sintética de View.

anoleTree<-read.tree("./0Data/PGLS_example/anolis.phy")

El árbol filogenético puede visualizarse con la función plot():

plot(anoleTree)

Figura 1. Árbol filogenético de las especies de lagartos del género \(Anolis\)

El paquete “geiger” tiene una función para comprobar que hay coincidencia entre los nombres del árbol y el “data frame”.

name.check(anoleTree, anoleData)
#> [1] "OK"

Usando métodos gráficos del base y el paquete “ape” observamos que existe una correlación negativa (inversa) entre las variables “awesomeness” y “hostility”.

plot(anoleData[,c("awesomeness", "hostility")])

Figura 2. Relación entre dos rasgos de las especies del género \(Anolis\)

Ahora sí construyo el modelo PGLS con el árbol filogenético ajustado por un modelo evolutivo browniano en el argumento para la estructura de correlación. El código debajo muestra el Uso de PGLS, con ape (ver observación en el Cuadro 4):

# Ajuste del modelo lineal mediante mínimos cuadrados generalizados con nlme:
anoleData$species <- rownames(anoleData) # agrego columna con nombre de las especies
anoleData2 <- anoleData[, c(11, 1:10)]

pglsModel <- gls(hostility~awesomeness,
                 correlation=corBrownian(1, phy=anoleTree, form=~species),
                 data=anoleData2,
                 method="ML")

#correlation=corBrownian(phy=anoleTree, form=~species), 

summary(pglsModel)
#> Generalized least squares fit by maximum likelihood
#>   Model: hostility ~ awesomeness 
#>   Data: anoleData2 
#>        AIC     BIC    logLik
#>   190.9775 198.793 -92.48875
#> 
#> Correlation Structure: corBrownian
#>  Formula: ~species 
#>  Parameter estimate(s):
#> numeric(0)
#> 
#> Coefficients:
#>                  Value  Std.Error    t-value p-value
#> (Intercept)  0.1505613 0.26262709   0.573289  0.5678
#> awesomeness -0.9775847 0.04515861 -21.647804  0.0000
#> 
#>  Correlation: 
#>             (Intr)
#> awesomeness -0.042
#> 
#> Standardized residuals:
#>         Min          Q1         Med          Q3         Max 
#> -0.76019997 -0.39056977 -0.04941607  0.19596725  1.07373699 
#> 
#> Residual standard error: 0.8877372 
#> Degrees of freedom: 100 total; 98 residual

## Para graficar el modelo:
# Selección de columnas de rasgos
host<-anoleData2[,"hostility"]
awe<-anoleData2[,"awesomeness"]

# Incluimos nombres de especies (filas)
names(host)<-names(awe)<-rownames(anoleData2)

coef(pglsModel) # extraigo coeficientes de recta de regresión
#> (Intercept) awesomeness 
#>   0.1505613  -0.9775847
plot(host~awe)  # 
# Uso ordenada al origen y pendiente para graficar recta del PGLS
abline(a=coef(pglsModel)[1], b=coef(pglsModel)[2])

Figura 3. Relación entre dos rasgos de las especies del género \(Anolis\)

Cuadro 4: Hago notar que al replicar el código del modelo PGLS del ejemplo original (el de la página de internet de donde saqué el ejemplo), R devuelve el mensaje de advertencia: “Warning: No covariate specified, species will be taken as ordered in the data frame. To avoid this message, specify a covariate containing the species names with the ‘form’ argument”, que advierte que, como no se especificó un vector con el orden de las especies, por defecto se consideró el orden de las especies tal como aparecen en el dataframe. Para evitar este mensaje de “warning” agregué un vector con el nombre de las especies como nueva columna en el data.frame original. En el argumento “form” debemos incluir una fórmula unilateral de la forma \(\sim t\), o \(\sim t|g\), especificando la covariable \(t\) de los taxones y, opcionalmente, un factor de agrupación \(g\). La covariable (especies) para esta estructura de correlación debe ser caracter, con entradas que coincidan con las etiquetas de las puntas en el árbol filogenético. Cuando un factor de agrupación está presente en “form”, se supone que la estructura de correlación se aplica sólo a las observaciones dentro del mismo nivel de agrupación. Se supone que las observaciones con diferentes niveles de agrupación no están correlacionadas. Por defecto es \(\sim 1\), lo que corresponde a utilizar el orden de las observaciones en los datos como una covariable, y sin grupos. Esto significa que \(t\) debe contener los nombres de las especies o etiquetas de las puntas, deletreados exactamente como están en su filogenia.

La función gls del paquete “nlme” tiene la flexibilidad necesaria para incluir un predictor discreto en el ajuste del PGLS:

pglsModel2<-gls(hostility~ecomorph,
                correlation=corBrownian(1, phy=anoleTree, form=~species),
                data=anoleData2,
                method="ML")
anova(pglsModel2)

coef(pglsModel2)
#> (Intercept)  ecomorphGB   ecomorphT  ecomorphTC  ecomorphTG  ecomorphTW 
#>   0.4843515  -0.6315992  -1.0585278  -0.8558138  -0.4085610  -0.4039460 
#>   ecomorphU 
#>  -0.7021719

Incluso podemos incluir múltiples predictores:

pglsModel3<-gls(hostility~ecomorph*awesomeness,
                correlation=corBrownian(1, phy=anoleTree, form=~species),
                data=anoleData2,
                method="ML")

anova(pglsModel3)

También podemos suponer que la estructura del error sigue un modelo “OU” en lugar de un movimiento browniano:

# No converge - ¡y esto es difícil de arreglar!
# pglsModelLambda <- gls(hostility~awesomeness,
#                        correlation=corPagel(1, phy=anoleTree,
#                                             form=~species, fixed=FALSE),
#                       data=anoleData2, method="ML")
# Se trata de un problema de escala. Podemos solucionarlo rápidamente alargando la longitud de las ramas. Esto no afectará al análisis, salvo el cambio de escala de un parámetro molesto.
tempTree<-anoleTree
tempTree$edge.length<-tempTree$edge.length * 100

pglsModelLambda<-gls(hostility~awesomeness,
                     correlation=corPagel(1, phy=tempTree,
                                          form=~species, fixed=FALSE),
                     data=anoleData2,
                     method="ML")
summary(pglsModelLambda)
#> Generalized least squares fit by maximum likelihood
#>   Model: hostility ~ awesomeness 
#>   Data: anoleData2 
#>        AIC      BIC    logLik
#>   72.56056 82.98124 -32.28028
#> 
#> Correlation Structure: corPagel
#>  Formula: ~species 
#>  Parameter estimate(s):
#>     lambda 
#> -0.1585633 
#> 
#> Coefficients:
#>                  Value  Std.Error    t-value p-value
#> (Intercept)  0.0612470 0.01581847   3.871868   2e-04
#> awesomeness -0.8776519 0.03104246 -28.272628   0e+00
#> 
#>  Correlation: 
#>             (Intr)
#> awesomeness -1    
#> 
#> Standardized residuals:
#>         Min          Q1         Med          Q3         Max 
#> -1.78946302 -0.71477505  0.00309539  0.78509306  2.23215144 
#> 
#> Residual standard error: 0.3709858 
#> Degrees of freedom: 100 total; 98 residual

pglsModelOU<-gls(hostility~awesomeness,
                 correlation=corMartins(1, phy=tempTree, form=~species),
                 data=anoleData2,
                 method="ML")
summary(pglsModelOU)
#> Generalized least squares fit by maximum likelihood
#>   Model: hostility ~ awesomeness 
#>   Data: anoleData2 
#>        AIC      BIC    logLik
#>   96.63478 107.0555 -44.31739
#> 
#> Correlation Structure: corMartins
#>  Formula: ~species 
#>  Parameter estimate(s):
#>    alpha 
#> 4.441625 
#> 
#> Coefficients:
#>                  Value  Std.Error    t-value p-value
#> (Intercept)  0.1084258 0.03952884   2.742954  0.0072
#> awesomeness -0.8811632 0.03657646 -24.090988  0.0000
#> 
#>  Correlation: 
#>             (Intr)
#> awesomeness -0.269
#> 
#> Standardized residuals:
#>        Min         Q1        Med         Q3        Max 
#> -1.8664557 -0.8132899 -0.1103815  0.6474918  2.0919152 
#> 
#> Residual standard error: 0.376904 
#> Degrees of freedom: 100 total; 98 residual

Otro PGLS con gls de nlme, ade4 y ape

Aquí damos un ejemplo de análisis PGLS utilizando un árbol y datos del paquete ade4. Analizamos los datos utilizando funciones de los paquetes “ape” y “nlme”. Existen otras formas de realizar análisis PGLS en R, sobre todo utilizando la función pgls del paquete caper (vean el ejemplo siguiente). El paquete “nlme” tiene la ventaja de ser un paquete maduro de propósito general con útiles funciones de “ayuda”, y que se incluye en la distribución básica de R.

El primer paso en cualquier análisis es importar el árbol y los datos a R. Para la mayoría de los propósitos, las funciones de R “read.table” y “read.tree” del paquete “ape” serán suficientes. El paquete “ape” también tiene una función para leer archivos Nexus, “read.nexus”, que es otro formato de archivo de árbol popular. Para fines de demostración en este tutorial, utilizaremos un árbol y datos del paquete “ade4”. Considere los siguientes datos de historia de vida y árbol (incluyendo longitudes de rama) para 18 especies de lagartos:

library(ade4) # fuente de datos de ejemplo y árbol
library(ape)  # manejo de árboles
library(nlme) # modelos de regresión
data(lizards)
tree <- read.tree(text = lizards$hprA)
dat <- lizards$traits[tree$tip.label, ] # ordena los datos según el árbol
plot(tree, main = "Phylogeny for 18 Lizard Species",
     direction = "up",
     srt = -90,
     label.offset = 1)

Figura 4. Árbol filogenético de las 18 especies de lagarto analizadas

head(dat) # sólo las primeras líneas de datos

Supongamos que estamos interesados en la regresión de la talla de madurez sobre la edad de madurez. Si asumiéramos un modelo de evolución de movimiento browniano, podríamos construir la matriz de correlación filogenética implícita en el árbol y utilizarla en una llamada a la función gls del paquete nlme:

mat <- vcv(tree, corr=TRUE) # construyo la matriz
fit <- gls(matur.L ~ age.mat,
           correlation=corSymm(mat[lower.tri(mat)], fixed=TRUE),
           data=dat)

Hay que tener en cuenta varias cosas:

• La función vcv calcula una matriz de covarianza por defecto, y una matriz de correlación con el argumento corr=TRUE.

• La función gls funciona como las otras funciones de regresión en R, excepto que tiene un argumento de correlación.

• El triángulo inferior de la matriz de correlación se pasa a la función corSymm, que especifica una matriz de correlación simétrica general. Puesto que no estamos estimando ningún parámetro en esta matriz, establecemos fixed=TRUE, indicando a gls que la matriz de correlación es fija mientras dure la optimización de parámetros.

• Este análisis supone que el árbol es ultramétrico, es decir, que todas sus puntas están alineadas. En este caso es así, pero es fácil imaginar casos en los que no sea así. (por ejemplo, la inclusión de especies extinguidas o determinadas transformaciones de la longitud de las ramas). La solución es utilizar el argumento “weights” de gls:

tip.heights <- diag(mat)
fit <- gls(matur.L ~ age.mat,
           correlation = corSymm(mat[lower.tri(mat)], fixed = TRUE),
           weights = varFixed(~tip.heights), data = dat)

El paquete ape dispone de otras estructuras de correlación para su uso con gls. La más sencilla es corBrownian, que se ajusta al mismo modelo que la anterior:

fit2 <- gls(matur.L ~ age.mat,
            correlation=corBrownian(phy=tree), data=dat)
anova(fit, fit2) # los modelos son equivalentes

Usar corBrownian es más sencillo que usar corSymm, ya que todo lo que necesita decirle es el nombre del árbol. No es necesario manipular directamente la matriz de correlaciones.

Podemos hacer algunos gráficos de diagnóstico básicos para ver si el modelo se aparta de los supuestos de regresión habituales. En particular, podemos observar los residuos frente a los valores ajustados para comprobar si hay heteroscedasticidad y/o curvatura, y un gráfico cuantil-cuantil de los residuos (normalizados) para comprobar si hay desviaciones del supuesto de residuos distribuidos normalmente:

plot(fit2, resid(., type="n")~fitted(.),
     main="Normalized Residuals v Fitted Values",
     abline=c(0,0))
res <- resid(fit2, type="n")
qqnorm(res)
qqline(res)
anova(fit, fit2) # los modelos son equivalentes

Figura 5. Diagnóstico gráfico de los residuos del modelo PGLS para comprobar normalidad y homocedasticidad

Obsérvese que los gráficos tienen muy buen aspecto, salvo quizá por un valor atípico (el que tiene un residual normalizado \(> 3\) (figura 5, izq.). Podríamos volver a examinar los datos originales para ver si hubo un error o algo inusual en esa especie. Una opción es identificar a esa especie y excluirla del análisis, y volver a ejecutar el análisis sin el valor atípico (especie “Pg”, vean el código debajo):

res[which.max(res)] # especie con mayor residuo (valor atípico O "outlier")
#>       Pg 
#> 3.762375

dat3 <- dat[-which(rownames(dat) == "Pg"),]
tree3 <- drop.tip(tree, "Pg")
fit3 <- gls(matur.L ~ age.mat, correlation=corBrownian(phy=tree3), data=dat3)
summary(fit3)
#> Generalized least squares fit by REML
#>   Model: matur.L ~ age.mat 
#>   Data: dat3 
#>        AIC      BIC    logLik
#>   138.6185 140.7426 -66.30924
#> 
#> Correlation Structure: corBrownian
#>  Formula: ~1 
#>  Parameter estimate(s):
#> numeric(0)
#> 
#> Coefficients:
#>                 Value Std.Error  t-value p-value
#> (Intercept) 27.615859 21.545175 1.281765  0.2194
#> age.mat      2.902111  1.436332 2.020501  0.0616
#> 
#>  Correlation: 
#>         (Intr)
#> age.mat -0.727
#> 
#> Standardized residuals:
#>        Min         Q1        Med         Q3        Max 
#> -0.5714836 -0.3228441 -0.1222878  0.1507580  1.7015545 
#> 
#> Residual standard error: 28.94061 
#> Degrees of freedom: 17 total; 15 residual

Cuadro 5: Noten que si ustedes adaptan el código de este ejemplo para correrlo con R en sus computadoras, cada vez que ajusten el PGLS con “gls”, recibirán el mensaje de advertencia (referido en el Cuadro 4): “Warning: No covariate specified, species will be taken as ordered in the data frame. To avoid this message, specify a covariate containing the species names with the ‘form’ argument”, que advierte que, como no se especificó un vector con el orden de las especies, por defecto se consideró el orden de las especies tal como aparecen en el dataframe. En el Cuadro 4 expliqué cómo hacer para evitar este mensaje de “warning” en el ejemplo anterior. No hice las modificaciones al código necesarias para evita el mensaje de advertencia en este y los siguientes ejemplos. En este tutorial, si no ven los mensajes de advertencia para éste, y varios de los siguientes ejemplos, es simplemente porque (por una cuestión de espacio y claridad) especifiqué en el código del RMarkdown que generó el tutorial que no me devuelva los “warnings”.

Vemos en la devolución del “summary” del PGLS ajustado (summary(fit3), fit3 es el PGLS sin la especie “Pg”) que la pendiente de la relación entre la longitud en la madurez y la edad en la madurez es de 2.902 (aunque, por poco, no es significativa: \(p = 0.0616\)). Esto significa que por cada unidad de aumento en la edad de madurez, hay un aumento asociado de casi 3 unidades en la longitud de madurez. Las especies de maduración más tardía son más largas en la madurez.

Quizá nos preocupe la cantidad de “señal” filogenética de nuestros datos. Es posible que los datos observados no hayan evolucionado según un modelo de movimiento browniano. En su lugar, podemos incluir un parámetro adicional en la matriz de correlación para tener en cuenta el error en la estimación de la longitud de las ramas. Nuestra nueva matriz es entonces una función del árbol y de este nuevo parámetro: \(Matriz_{Nueva} = f(parámetro.árbol)\).

El paquete “ape” proporciona cuatro estructuras de correlación de un solo parámetro para su uso con la función “gls”. Se trata de “corPagel”, “corGrafen”, “corMartins” y “corBlomberg”. Cada estructura de correlación construye la matriz de correlación de una manera diferente. La más sencilla es “corPagel”, en la que el parámetro (\(\lambda\)) multiplica los diagonales de la matriz de correlación del movimiento browniano.

El \(\lambda\) puede oscilar entre cero (sin señal filogenética) y uno (coherente con el movimiento browniano). Los valores intermedios implican que los datos apoyan un modelo que está entre las dos situaciones anteriores. El \(\lambda\) representa la cantidad de “señal” filogenética en los datos, condicional en el árbol. En un contexto de regresión, podemos utilizar el método de máxima verosimilitud (restringida) para ajustar simultáneamente el modelo de regresión y estimar \(\lambda\), que representa la señal filogenética en los residuos. También podemos obtener un intervalo de confianza para \(\lambda\) y realizar pruebas de las hipótesis \(\lambda = 0\) y \(\lambda = 1\).

fitPagel <- gls(matur.L ~ age.mat,
                correlation = corPagel(value=0.8,
                                       phy=tree3),
                data=dat3)
intervals(fitPagel, which="var-cov")
#> Approximate 95% confidence intervals
#> 
#>  Correlation structure:
#>            lower      est.    upper
#> lambda 0.4899848 0.8989786 1.307972
#> 
#>  Residual standard error:
#>    lower     est.    upper 
#> 11.76031 21.88207 40.71533

Podemos obtener pruebas directas de las hipótesis \(\lambda = 0\) (independencia) y \(\lambda = 1\) (movimiento browniano) ajustando modelos forzados a tener \(\lambda = 0\) (independencia) y \(\lambda = 1\) (movimiento browniano) y luego comparándolos con una prueba de razón de verosimilitud \(\chi^2\).

fitPagel0 <- gls(matur.L ~ age.mat,
                 correlation = corPagel(value = 0,
                                        phy = tree3,
                                        fixed = TRUE),
                 data = dat3) # independencia (lambda = 0)
fitPagel1 <- gls(matur.L ~ age.mat,
                 correlation = corPagel(value = 1,
                                        phy = tree3,
                                        fixed = TRUE),
                 data = dat3) # movimiento Browniano (lambda = 1)

anova(fitPagel, fitPagel0)

anova(fitPagel, fitPagel1)

Vemos que las pruebas de la razón de verosimilitud \(\chi^2\) que comparan nuestro modelo de \(\lambda\) “flotante” con un modelo que tiene \(\lambda\) fijado en cero o en uno tienen valores \(p\) de 0.5576 y 0.4987, respectivamente. Por tanto, no podemos rechazar la hipótesis nula de “residuos independientes” o movimiento browniano. Entonces, ¿por qué el intervalo de confianza de \(\lambda\) no incluye el cero? La estimación correcta de \(\lambda\) es difícil con un tamaño de muestra tan pequeño. Esto puede verse ejecutando el análisis gls con diferentes valores de partida para \(\lambda\). La \(\lambda\) óptima converge a valores que son o bien 0.899 o bien menores que cero. Se puede obtener una impresión de la estimación de \(\lambda\) trazando la probabilidad del modelo para varios valores fijos de \(\lambda\) (figura 6):

lambda <- seq(0, 1, length.out = 500)
lik <- sapply(lambda, function(lambda)
  logLik(gls(matur.L ~ age.mat,
             correlation = corPagel(value = lambda, phy = tree3, fixed = TRUE), data = dat3)))

plot(lik ~ lambda,
     type = "l",
     main = expression(paste("Likelihood Plot for ", lambda)),
     ylab = "Log Likelihood", xlab = expression(lambda))
abline(v = fitPagel$modelStruct, col = "red")

Figura 6. Distribución de probabilidad de valores de lambda óptimos para el modelo PGLS ajsutado

Podemos ver que para valores iniciales de \(\lambda\) inferiores a 0.7 aproximadamente, la optimización numérica se dibujará hacia \(\lambda \le 0\). En caso contrario, la optimización se dibujará hacia \(\lambda = 0.899\). Este valor está indicado por la línea roja. Con datos para más especies, podríamos estimar mejor \(\lambda\). Esto es una advertencia sobre la dificultad de la estimación de parámetros con tamaños de muestra pequeños.

Nos gustaría construir un gráfico de la relación entre la longitud en la madurez y la edad en la madurez, e incluir la línea de tendencia (figura 7):

with(dat3, plot(matur.L ~ age.mat,
                xlab = "Age at Maturity",
                ylab = "Length at Maturity",
                main = "Length at maturity versus \nage at maturity for 17 lizard species"))
abline(fitPagel)

Figura 7. Relación lineal entre la longitud al alcanzar la madurez en función de la edad al alcanzar la madurez

Hasta ahora nos hemos centrado en la estimación de \(\lambda\) como parte de un análisis de regresión. También podemos querer conocer el valor de \(\lambda\) para un único rasgo. Esto es sólo una cuestión de ajuste de un modelo de sólo intercepto. Aquí lo hacemos para todas las variables del “data frame”:

# hacemos una matriz para guardar los resultados
mat <- matrix(NA, ncol=3, nrow=(dim(dat3)[2]),
              dimnames=list(variable=names(dat3),
                            c("Lower" , "Estimate", "Upper")))
for (i in 1:(dim(dat)[2])) { # bucle a través de cada variable
  form <- formula(paste(names(dat)[i], " ~ 1")) # construir la fórmula del modelo
  this.fit <- gls(form,
                  data=dat,
                  correlation=corPagel(0.1, phy=tree),
                  control=glsControl(opt="optim")) # ajustar el modelo
  ints <- try(intervals(this.fit)$corStruct) # obtener lambda e intervalo de confianza
# si hay un error, simplemente obtenga la estimación puntual y establezca los límites de confianza en NA
  if (inherits(ints, "try-error"))
    ints <- c(NA, coef(this.fit$modelStruct), NA)
  mat[i,] <- ints # guardar los resultados en la matriz
  } # finaliza el bucle
#> Error in intervals.gls(this.fit) : 
#>   cannot get confidence intervals on var-cov components: Non-positive definite approximate variance-covariance
#> Error in intervals.gls(this.fit) : 
#>   cannot get confidence intervals on var-cov components: Non-positive definite approximate variance-covariance
## Lambda para cada variable:
mat
#>           
#> variable       Lower  Estimate    Upper
#>   mean.L   0.9501762 0.9819832 1.013790
#>   matur.L  0.9382276 0.9779053 1.017583
#>   max.L    0.9570556 0.9853112 1.013567
#>   hatch.L  0.8832186 0.9582008 1.033183
#>   hatch.m  0.9780110 0.9929446 1.007878
#>   clutch.S 0.8927188 0.9636695 1.034620
#>   age.mat         NA 1.0005152       NA
#>   clutch.F        NA 1.0000167       NA

Hay que tener en cuenta algunas cosas:

• Construimos la fórmula del modelo “sobre la marcha” utilizando las funciones “formula” y “paste”. Esta es una técnica útil si necesita ajustar muchas regresiones diferentes utilizando distintas combinaciones de variables de respuesta y/o explicativas.

• Utilizamos una llamada a “glsControl” para cambiar el optimizador del ajuste. A veces esto puede mejorar la convergencia.

• Utilizamos la función “try” para detectar errores. A veces se puede obtener una estimación puntual para \(\lambda\), pero si la matriz de varianza-covarianza asociada no es definida para valores positivos, no será posible estimar los límites de confianza (y habrá un error). Por lo tanto, le dijimos a R que si había un error, sólo informara de la estimación puntual para \(\lambda\) y dejara los límites de confianza como “NA”, es decir, valores no disponibles o faltantes.

Las estimaciones de \(\lambda\) para todas las variables son muy cercanas a uno, lo que indica una “señal” filogenética fuerte en estas variables.

En raras ocasiones, habrá información externa a los datos relativa al valor de \(\lambda\). Es decir, \(\lambda\) puede ser conocido. Es fácil utilizar un valor preespecificado para \(\lambda\). Basta con pasarlo a la función “corPagel” con el argumento “value”, junto con el argumento “fixed=TRUE”:

new.lambda <- 0.7 # valor lambda predeterminado
fitPagel4 <- gls(matur.L ~ age.mat,
                 correlation = corPagel(value = new.lambda,
                                        fixed = TRUE,
                                        phy = tree3),
                 data = dat3)
fitPagel4
#> Generalized least squares fit by REML
#>   Model: matur.L ~ age.mat 
#>   Data: dat3 
#>   Log-restricted-likelihood: -66.15628
#> 
#> Coefficients:
#> (Intercept)     age.mat 
#>    9.886470    4.531729 
#> 
#> Correlation Structure: corPagel
#>  Formula: ~1 
#>  Parameter estimate(s):
#> lambda 
#>    0.7 
#> Degrees of freedom: 17 total; 15 residual
#> Residual standard error: 18.73424

Noten que es un error estimar primero \(\lambda\) para un rasgo y luego utilizarlo en un análisis posterior del mismo conjunto de datos. Esto se debe a que la estimación (restringida) de máxima verosimilitud de \(\lambda\) para un único rasgo puede ser una mala estimación de \(\lambda\) cuando se aplica a los residuos de un modelo más complejo, como un modelo de regresión. Tiene sentido estimar \(\lambda\) y los parámetros de regresión juntos para que la \(\lambda\) que se estima sea el estimador de máxima verosimilitud (restringido) para los datos y el modelo completos. Es útil saber cómo cambia el árbol cuando se le aplican diferentes valores de \(\lambda\). A continuación se muestran seis versiones del árbol original después de la transformación (figura 8). El código utiliza la función “vcv2phylo” para convertir de nuevo de una matriz a un árbol, y requiere el paquete “matrixcalc”:

library(matrixcalc)
source("vcv2phylo.R") # cargar el código de vcv2phylo
## Para cargarlo, antes debemos crear (copíar) el archivo vcv2phylo.R de la función en el directorio de trabajo
ll <- seq(1, 0, by = -0.2)
mats <- sapply(ll,
               function(x) corMatrix(Initialize(corPagel(x, phy = tree),
                                                data = dat)),
               simplify = FALSE)

trees <- lapply(mats, vcv2phylo)
class(trees) <- "multiPhylo"

par(mfrow = c(2, 3), mar = rep(1, 4))

for (i in 1:length(ll))
  plot(trees[[i]],
       main = substitute(paste(lambda, " = ", ll),
                         list(ll = ll[i])),
       direction = "upwards",
       show.tip.label = FALSE)

Figura 8. Versiones del árbol filogenético original a partir de diferentes valores de lambda

Podemos ver que \(\lambda = 1\) no es más que el árbol original. Valores decrecientes de \(\lambda = 1\) alargan las ramas terminales y comprimen las ramas internas (figura 8). \(\lambda = 0\) reduce todas las ramas internas a longitud cero, produciendo una filogenia en estrella. Queda como ejercicio mostrar el efecto de otras transformaciones de un solo parámetro en el árbol.

PGLS con caper

El tercer ejemplo proviene del material práctico del capítulo 7.5 (“Studying Correlated Evolution with Phylogenetic Generalized Least Squares (PGLS)”) del libro “The Comparative Method in Evolutionary Anthropology and Biology” de Charles Nunn, disponible en un “wiki” de Duke University, ver haciendo click. El capítulo hace uso del paquete “caper”. El ejemplo implementado por Natalie Cooper y Charlie Nunn muestra cómo incorporar al PGLS las tres medidas de la señal filogenética vistas (\(\lambda, \kappa, \delta\)), también denominados parámetros de escalamiento, para mejorar significativamente el ajuste de los datos al modelo y, por tanto, la estimación de la correlación de rasgos. Si ninguno de los parámetros de escalamiento se aparta de 1.0 (lo que implica un movimiento browniano), GLS devolverá resultados idénticos a los contrastes independientes. Sin embargo, si la suposición de movimiento browniano es incorrecta, entonces estos valores se apartarán de 1.0 y la incorporación de los parámetros de escalamiento en el análisis dará lugar a estimaciones de parámetros diferentes que los contrastes indepemdientes para parámetros como la pendiente que relaciona dos variables.

Para empezar, instalen y carguen el paquete de R ‘caper’ y descarguen y lean los datos de las especies y el árbol filogenético:

# install.packages('caper')
library("caper")
primatedata <-read.table("./0Data/PGLS_example/Primatedata.txt", sep = "\t", header = TRUE)
primatetree <-read.nexus("./0Data/PGLS_example/consensusTree_10kTrees_Version2.nex")

Ahora podemos visualizar el árbol filogenético:

plot(primatetree)

Figura 9. Árbol filogenético de 226 especies de primates

Como en el primer ejemplo, verificamos que los nombres de las especies del árbol y el data frame coincidan (veremos que, al contrario del primer ejemplo, existen varias diferencias):

name.check(primatetree, primatedata)
#> $tree_not_data
#>   [1] "Allenopithecus_nigroviridis"                  
#>   [2] "Allocebus_trichotis"                          
#>   [3] "Alouatta_caraya"                              
#>   [4] "Alouatta_palliata"                            
#>   [5] "Alouatta_pigra"                               
#>   [6] "Alouatta_sara"                                
#>   [7] "Alouatta_seniculus"                           
#>   [8] "Aotus_azarae"                                 
#>   [9] "Aotus_azarae_infulatus"                       
#>  [10] "Aotus_lemurinus_griseimembra"                 
#>  [11] "Aotus_nancymaae"                              
#>  [12] "Aotus_trivirgatus"                            
#>  [13] "Arctocebus_aureus"                            
#>  [14] "Arctocebus_calabarensis"                      
#>  [15] "Ateles_belzebuth"                             
#>  [16] "Ateles_fusciceps"                             
#>  [17] "Ateles_geoffroyi"                             
#>  [18] "Ateles_geoffroyi_ornatus"                     
#>  [19] "Ateles_geoffroyi_vellerosus"                  
#>  [20] "Ateles_paniscus"                              
#>  [21] "Avahi_laniger"                                
#>  [22] "Avahi_occidentalis"                           
#>  [23] "Brachyteles_arachnoides"                      
#>  [24] "Bunopithecus_hoolock"                         
#>  [25] "Callicebus_donacophilus"                      
#>  [26] "Callicebus_moloch"                            
#>  [27] "Callimico_goeldii"                            
#>  [28] "Callithrix_(Mico)_emiliae"                    
#>  [29] "Callithrix_argentata"                         
#>  [30] "Callithrix_aurita"                            
#>  [31] "Callithrix_geoffroyi"                         
#>  [32] "Callithrix_humeralifera"                      
#>  [33] "Callithrix_jacchus"                           
#>  [34] "Callithrix_kuhli"                             
#>  [35] "Callithrix_penicillata"                       
#>  [36] "Callithrix_pygmaea"                           
#>  [37] "Cebus_albifrons"                              
#>  [38] "Cebus_apella"                                 
#>  [39] "Cebus_capucinus"                              
#>  [40] "Cercocebus_agilis"                            
#>  [41] "Cercocebus_atys"                              
#>  [42] "Cercocebus_galeritus"                         
#>  [43] "Cercocebus_torquatus"                         
#>  [44] "Cercopithecus_ascanius"                       
#>  [45] "Cercopithecus_cephus"                         
#>  [46] "Cercopithecus_cephus_cephus"                  
#>  [47] "Cercopithecus_cephus_ngottoensis"             
#>  [48] "Cercopithecus_diana"                          
#>  [49] "Cercopithecus_erythrogaster_erythrogaster"    
#>  [50] "Cercopithecus_erythrotis"                     
#>  [51] "Cercopithecus_hamlyni"                        
#>  [52] "Cercopithecus_lhoesti"                        
#>  [53] "Cercopithecus_lowei"                          
#>  [54] "Cercopithecus_mitis"                          
#>  [55] "Cercopithecus_mona"                           
#>  [56] "Cercopithecus_neglectus"                      
#>  [57] "Cercopithecus_nictitans"                      
#>  [58] "Cercopithecus_petaurista"                     
#>  [59] "Cercopithecus_preussi"                        
#>  [60] "Cercopithecus_solatus"                        
#>  [61] "Cercopithecus_wolfi"                          
#>  [62] "Cheirogaleus_crossleyi"                       
#>  [63] "Cheirogaleus_major"                           
#>  [64] "Cheirogaleus_medius"                          
#>  [65] "Chlorocebus_aethiops"                         
#>  [66] "Chlorocebus_pygerythrus"                      
#>  [67] "Chlorocebus_sabaeus"                          
#>  [68] "Chlorocebus_tantalus"                         
#>  [69] "Colobus_angolensis"                           
#>  [70] "Colobus_guereza"                              
#>  [71] "Colobus_polykomos"                            
#>  [72] "Daubentonia_madagascariensis"                 
#>  [73] "Erythrocebus_patas"                           
#>  [74] "Eulemur_albifrons"                            
#>  [75] "Eulemur_albocollaris"                         
#>  [76] "Eulemur_collaris"                             
#>  [77] "Eulemur_coronatus"                            
#>  [78] "Eulemur_fulvus"                               
#>  [79] "Eulemur_macaco_flavifrons"                    
#>  [80] "Eulemur_macaco_macaco"                        
#>  [81] "Eulemur_mongoz"                               
#>  [82] "Eulemur_rubriventer"                          
#>  [83] "Eulemur_rufus"                                
#>  [84] "Eulemur_sanfordi"                             
#>  [85] "Euoticus_elegantulus"                         
#>  [86] "Galago_alleni"                                
#>  [87] "Galago_demidoff"                              
#>  [88] "Galago_gallarum"                              
#>  [89] "Galago_moholi"                                
#>  [90] "Galago_senegalensis"                          
#>  [91] "Galago_zanzibaricus"                          
#>  [92] "Gorilla_gorilla_gorilla"                      
#>  [93] "Hapalemur_alaotrensis"                        
#>  [94] "Hapalemur_aureus"                             
#>  [95] "Hapalemur_griseus_griseus"                    
#>  [96] "Hapalemur_griseus_meridionalis"               
#>  [97] "Hapalemur_occidentalis"                       
#>  [98] "Homo_sapiens"                                 
#>  [99] "Hylobates_agilis"                             
#> [100] "Hylobates_klossii"                            
#> [101] "Hylobates_lar"                                
#> [102] "Hylobates_moloch"                             
#> [103] "Hylobates_muelleri"                           
#> [104] "Hylobates_pileatus"                           
#> [105] "Indri_indri"                                  
#> [106] "Lagothrix_lagotricha"                         
#> [107] "Lemur_catta"                                  
#> [108] "Leontopithecus_chrysomelas"                   
#> [109] "Leontopithecus_chrysopygus"                   
#> [110] "Leontopithecus_rosalia"                       
#> [111] "Lepilemur_aeeclis"                            
#> [112] "Lepilemur_ankaranensis"                       
#> [113] "Lepilemur_dorsalis"                           
#> [114] "Lepilemur_edwardsi"                           
#> [115] "Lepilemur_leucopus"                           
#> [116] "Lepilemur_microdon"                           
#> [117] "Lepilemur_mitsinjoensis"                      
#> [118] "Lepilemur_mustelinus"                         
#> [119] "Lepilemur_randrianasoli"                      
#> [120] "Lepilemur_ruficaudatus"                       
#> [121] "Lepilemur_sahamalazensis"                     
#> [122] "Lepilemur_seali"                              
#> [123] "Lepilemur_septentrionalis"                    
#> [124] "Lophocebus_albigena"                          
#> [125] "Lophocebus_aterrimus"                         
#> [126] "Loris_lydekkerianus_grandis"                  
#> [127] "Loris_lydekkerianus_malabaricus"              
#> [128] "Loris_tardigradus"                            
#> [129] "Macaca_arctoides"                             
#> [130] "Macaca_assamensis"                            
#> [131] "Macaca_cyclopis"                              
#> [132] "Macaca_fascicularis"                          
#> [133] "Macaca_fuscata"                               
#> [134] "Macaca_hecki"                                 
#> [135] "Macaca_leonina"                               
#> [136] "Macaca_maura"                                 
#> [137] "Macaca_mulatta"                               
#> [138] "Macaca_nemestrina"                            
#> [139] "Macaca_nigra"                                 
#> [140] "Macaca_nigrescens"                            
#> [141] "Macaca_ochreata"                              
#> [142] "Macaca_ochreata_brunnescens"                  
#> [143] "Macaca_pagensis"                              
#> [144] "Macaca_radiata"                               
#> [145] "Macaca_siberu"                                
#> [146] "Macaca_silenus"                               
#> [147] "Macaca_sinica"                                
#> [148] "Macaca_sylvanus"                              
#> [149] "Macaca_thibetana"                             
#> [150] "Macaca_tonkeana"                              
#> [151] "Mandrillus_leucophaeus"                       
#> [152] "Mandrillus_sphinx"                            
#> [153] "Microcebus_berthae"                           
#> [154] "Microcebus_bongolavensis"                     
#> [155] "Microcebus_danfossi"                          
#> [156] "Microcebus_griseorufus"                       
#> [157] "Microcebus_jollyae"                           
#> [158] "Microcebus_lehilahytsara"                     
#> [159] "Microcebus_lokobensis"                        
#> [160] "Microcebus_mittermeieri"                      
#> [161] "Microcebus_murinus"                           
#> [162] "Microcebus_myoxinus"                          
#> [163] "Microcebus_ravelobensis"                      
#> [164] "Microcebus_rufus"                             
#> [165] "Microcebus_sambiranensis"                     
#> [166] "Microcebus_simmonsi"                          
#> [167] "Microcebus_tavaratra"                         
#> [168] "Miopithecus_talapoin"                         
#> [169] "Mirza_coquereli"                              
#> [170] "Nasalis_larvatus"                             
#> [171] "Nomascus_concolor"                            
#> [172] "Nomascus_gabriellae"                          
#> [173] "Nomascus_leucogenys"                          
#> [174] "Nycticebus_coucang"                           
#> [175] "Nycticebus_pygmaeus"                          
#> [176] "Otolemur_crassicaudatus"                      
#> [177] "Otolemur_garnettii"                           
#> [178] "Pan_paniscus"                                 
#> [179] "Pan_troglodytes_schweinfurthii"               
#> [180] "Pan_troglodytes_troglodytes"                  
#> [181] "Pan_troglodytes_verus"                        
#> [182] "Papio_anubis"                                 
#> [183] "Papio_cynocephalus"                           
#> [184] "Papio_hamadryas"                              
#> [185] "Papio_papio"                                  
#> [186] "Papio_ursinus"                                
#> [187] "Perodicticus_potto"                           
#> [188] "Piliocolobus_badius"                          
#> [189] "Pongo_abelii"                                 
#> [190] "Pongo_pygmaeus_pygmaeus"                      
#> [191] "Presbytis_melalophos"                         
#> [192] "Prolemur_simus"                               
#> [193] "Propithecus_coquereli"                        
#> [194] "Propithecus_diadema"                          
#> [195] "Propithecus_edwardsi"                         
#> [196] "Propithecus_tattersalli"                      
#> [197] "Propithecus_verreauxi"                        
#> [198] "Pygathrix_nemaeus"                            
#> [199] "Rhinopithecus_avunculus"                      
#> [200] "Rhinopithecus_bieti"                          
#> [201] "Rhinopithecus_brelichi"                       
#> [202] "Rhinopithecus_roxellana"                      
#> [203] "Rungwecebus_kipunji"                          
#> [204] "Saguinus_fuscicollis"                         
#> [205] "Saguinus_geoffroyi"                           
#> [206] "Saguinus_imperator"                           
#> [207] "Saguinus_midas"                               
#> [208] "Saguinus_oedipus"                             
#> [209] "Saimiri_boliviensis_boliviensis"              
#> [210] "Saimiri_oerstedii"                            
#> [211] "Saimiri_sciureus"                             
#> [212] "Semnopithecus_entellus"                       
#> [213] "Symphalangus_syndactylus"                     
#> [214] "Tarsius_bancanus"                             
#> [215] "Tarsius_syrichta"                             
#> [216] "Theropithecus_gelada"                         
#> [217] "Trachypithecus_(Trachypithecus)_auratus"      
#> [218] "Trachypithecus_(Trachypithecus)_poliocephalus"
#> [219] "Trachypithecus_cristatus"                     
#> [220] "Trachypithecus_francoisi"                     
#> [221] "Trachypithecus_johnii"                        
#> [222] "Trachypithecus_obscurus"                      
#> [223] "Trachypithecus_phayrei"                       
#> [224] "Trachypithecus_pileatus"                      
#> [225] "Varecia_rubra"                                
#> [226] "Varecia_variegata_variegata"                  
#> 
#> $data_not_tree
#>  [1] "1"  "10" "11" "12" "13" "14" "15" "16" "17" "18" "19" "2"  "20" "21" "22"
#> [16] "23" "24" "25" "26" "27" "28" "29" "3"  "30" "31" "32" "33" "34" "35" "36"
#> [31] "37" "38" "39" "4"  "40" "41" "42" "43" "44" "45" "46" "47" "48" "49" "5" 
#> [46] "50" "51" "52" "53" "54" "55" "56" "57" "58" "59" "6"  "60" "61" "62" "63"
#> [61] "64" "65" "66" "67" "68" "69" "7"  "70" "71" "72" "73" "74" "75" "76" "77"
#> [76] "78" "79" "8"  "80" "81" "82" "83" "84" "85" "86" "87" "88" "89" "9"  "90"
#> [91] "91"

Observen que las filogenias en R no pueden tener espacios en las etiquetas de las puntas, en su lugar R utiliza guiones bajos para separar Género y especie (Genus_species). Los nombres de las especies en el árbol deben coincidir con los de los datos, por lo que en este conjunto de datos los espacios en los nombres de las especies se han sustituido por guiones bajos.

primatedata$Binomial<-gsub(" ", "_", primatedata$Binomial)

caper requiere que su filogenia y tus datos estén en un tipo especial de objeto R llamado objeto de datos comparativo (“comparative data object”). Esto hará coincidir los nombres de las especies en el árbol con los de sus datos de forma automática por lo que no hay necesidad de preocuparse de que los datos estén ordenados correctamente, etc. “phy” es tu árbol, “data” tu conjunto de datos, y “names.col” es el nombre de la columna que contiene los nombres de tus especies.

Cuadro 6: sobre los argumentos de esta función, “vcv = TRUE” almacena una matriz de varianza-covarianza de su árbol, que necesitará para pgls. “na.omit = FALSE” impide que la función elimine especies sin datos para ciertas variables. “warn.dropped = TRUE” te dirá si hay especies que no están en el árbol ni en el conjunto de datos y que, por tanto, se eliminan del objeto de datos.

Si eliminan especies, pueden ver la lista de especies eliminadas escribiendo “primate$dropped”. Asegúrense de comprobar que esta lista es la que esperaba, puede revelar errores tipográficos en los nombres de las especies. Si desean desactivar esta función, utilicen “warn.dropped = FALSE”. Tengan en cuenta que la función “comparative.data” hace todo el emparejamiento de datos y árbol por nosotros/as. Consideren además que esperamos descartar algunas especies aquí porque no todas las especies en “primatetree” están en “primatedata” y viceversa.

primate <- comparative.data(phy = primatetree,
                            data = primatedata,
                            names.col = Binomial,
                            vcv = TRUE,
                            na.omit = FALSE,
                            warn.dropped = TRUE)

primate$dropped
#> $tips
#>   [1] "Allenopithecus_nigroviridis"                  
#>   [2] "Cercopithecus_cephus_cephus"                  
#>   [3] "Cercopithecus_cephus_ngottoensis"             
#>   [4] "Cercopithecus_diana"                          
#>   [5] "Cercopithecus_erythrogaster_erythrogaster"    
#>   [6] "Cercopithecus_erythrotis"                     
#>   [7] "Cercopithecus_hamlyni"                        
#>   [8] "Cercopithecus_lhoesti"                        
#>   [9] "Cercopithecus_lowei"                          
#>  [10] "Cercopithecus_mona"                           
#>  [11] "Cercopithecus_petaurista"                     
#>  [12] "Cercopithecus_preussi"                        
#>  [13] "Cercopithecus_solatus"                        
#>  [14] "Cercopithecus_wolfi"                          
#>  [15] "Chlorocebus_aethiops"                         
#>  [16] "Chlorocebus_pygerythrus"                      
#>  [17] "Chlorocebus_sabaeus"                          
#>  [18] "Chlorocebus_tantalus"                         
#>  [19] "Allocebus_trichotis"                          
#>  [20] "Avahi_laniger"                                
#>  [21] "Avahi_occidentalis"                           
#>  [22] "Cheirogaleus_crossleyi"                       
#>  [23] "Eulemur_albifrons"                            
#>  [24] "Eulemur_albocollaris"                         
#>  [25] "Eulemur_collaris"                             
#>  [26] "Eulemur_macaco_flavifrons"                    
#>  [27] "Eulemur_macaco_macaco"                        
#>  [28] "Eulemur_rubriventer"                          
#>  [29] "Eulemur_rufus"                                
#>  [30] "Eulemur_sanfordi"                             
#>  [31] "Hapalemur_alaotrensis"                        
#>  [32] "Hapalemur_aureus"                             
#>  [33] "Hapalemur_griseus_griseus"                    
#>  [34] "Hapalemur_griseus_meridionalis"               
#>  [35] "Hapalemur_occidentalis"                       
#>  [36] "Indri_indri"                                  
#>  [37] "Lepilemur_aeeclis"                            
#>  [38] "Lepilemur_ankaranensis"                       
#>  [39] "Lepilemur_dorsalis"                           
#>  [40] "Lepilemur_edwardsi"                           
#>  [41] "Lepilemur_microdon"                           
#>  [42] "Lepilemur_mitsinjoensis"                      
#>  [43] "Lepilemur_randrianasoli"                      
#>  [44] "Lepilemur_ruficaudatus"                       
#>  [45] "Lepilemur_sahamalazensis"                     
#>  [46] "Lepilemur_seali"                              
#>  [47] "Lepilemur_septentrionalis"                    
#>  [48] "Microcebus_berthae"                           
#>  [49] "Microcebus_bongolavensis"                     
#>  [50] "Microcebus_danfossi"                          
#>  [51] "Microcebus_griseorufus"                       
#>  [52] "Microcebus_jollyae"                           
#>  [53] "Microcebus_lehilahytsara"                     
#>  [54] "Microcebus_lokobensis"                        
#>  [55] "Microcebus_mittermeieri"                      
#>  [56] "Microcebus_myoxinus"                          
#>  [57] "Microcebus_ravelobensis"                      
#>  [58] "Microcebus_sambiranensis"                     
#>  [59] "Microcebus_simmonsi"                          
#>  [60] "Microcebus_tavaratra"                         
#>  [61] "Propithecus_coquereli"                        
#>  [62] "Propithecus_diadema"                          
#>  [63] "Propithecus_edwardsi"                         
#>  [64] "Propithecus_tattersalli"                      
#>  [65] "Varecia_rubra"                                
#>  [66] "Varecia_variegata_variegata"                  
#>  [67] "Alouatta_caraya"                              
#>  [68] "Alouatta_sara"                                
#>  [69] "Ateles_fusciceps"                             
#>  [70] "Ateles_geoffroyi_ornatus"                     
#>  [71] "Ateles_geoffroyi_vellerosus"                  
#>  [72] "Brachyteles_arachnoides"                      
#>  [73] "Aotus_azarae"                                 
#>  [74] "Aotus_azarae_infulatus"                       
#>  [75] "Aotus_lemurinus_griseimembra"                 
#>  [76] "Aotus_nancymaae"                              
#>  [77] "Callithrix_(Mico)_emiliae"                    
#>  [78] "Callithrix_argentata"                         
#>  [79] "Callithrix_aurita"                            
#>  [80] "Callithrix_geoffroyi"                         
#>  [81] "Callithrix_humeralifera"                      
#>  [82] "Callithrix_kuhli"                             
#>  [83] "Callithrix_penicillata"                       
#>  [84] "Leontopithecus_chrysomelas"                   
#>  [85] "Leontopithecus_chrysopygus"                   
#>  [86] "Saguinus_geoffroyi"                           
#>  [87] "Saguinus_imperator"                           
#>  [88] "Saimiri_boliviensis_boliviensis"              
#>  [89] "Saimiri_oerstedii"                            
#>  [90] "Arctocebus_aureus"                            
#>  [91] "Arctocebus_calabarensis"                      
#>  [92] "Loris_lydekkerianus_grandis"                  
#>  [93] "Loris_lydekkerianus_malabaricus"              
#>  [94] "Nycticebus_coucang"                           
#>  [95] "Nycticebus_pygmaeus"                          
#>  [96] "Bunopithecus_hoolock"                         
#>  [97] "Gorilla_gorilla_gorilla"                      
#>  [98] "Homo_sapiens"                                 
#>  [99] "Hylobates_agilis"                             
#> [100] "Hylobates_klossii"                            
#> [101] "Hylobates_moloch"                             
#> [102] "Hylobates_muelleri"                           
#> [103] "Nomascus_gabriellae"                          
#> [104] "Nomascus_leucogenys"                          
#> [105] "Pan_troglodytes_schweinfurthii"               
#> [106] "Pan_troglodytes_troglodytes"                  
#> [107] "Pan_troglodytes_verus"                        
#> [108] "Pongo_abelii"                                 
#> [109] "Pongo_pygmaeus_pygmaeus"                      
#> [110] "Callicebus_donacophilus"                      
#> [111] "Cercocebus_agilis"                            
#> [112] "Cercocebus_atys"                              
#> [113] "Cercocebus_torquatus"                         
#> [114] "Lophocebus_aterrimus"                         
#> [115] "Macaca_arctoides"                             
#> [116] "Macaca_assamensis"                            
#> [117] "Macaca_cyclopis"                              
#> [118] "Macaca_hecki"                                 
#> [119] "Macaca_leonina"                               
#> [120] "Macaca_maura"                                 
#> [121] "Macaca_nigra"                                 
#> [122] "Macaca_nigrescens"                            
#> [123] "Macaca_ochreata"                              
#> [124] "Macaca_ochreata_brunnescens"                  
#> [125] "Macaca_pagensis"                              
#> [126] "Macaca_siberu"                                
#> [127] "Macaca_thibetana"                             
#> [128] "Macaca_tonkeana"                              
#> [129] "Mandrillus_leucophaeus"                       
#> [130] "Papio_papio"                                  
#> [131] "Rungwecebus_kipunji"                          
#> [132] "Colobus_angolensis"                           
#> [133] "Piliocolobus_badius"                          
#> [134] "Presbytis_melalophos"                         
#> [135] "Pygathrix_nemaeus"                            
#> [136] "Rhinopithecus_avunculus"                      
#> [137] "Rhinopithecus_bieti"                          
#> [138] "Rhinopithecus_brelichi"                       
#> [139] "Rhinopithecus_roxellana"                      
#> [140] "Trachypithecus_(Trachypithecus)_auratus"      
#> [141] "Trachypithecus_(Trachypithecus)_poliocephalus"
#> [142] "Trachypithecus_cristatus"                     
#> [143] "Trachypithecus_francoisi"                     
#> [144] "Trachypithecus_johnii"                        
#> [145] "Trachypithecus_phayrei"                       
#> [146] "Trachypithecus_pileatus"                      
#> [147] "Euoticus_elegantulus"                         
#> [148] "Galago_gallarum"                              
#> [149] "Galago_zanzibaricus"                          
#> 
#> $unmatched.rows
#>  [1] "Cercopithecus_campbelli" "Cercopithecus_pogonias" 
#>  [3] "Eulemur_macaco"          "Chiropotes_albinasus"   
#>  [5] "Chiropotes_satanas"      "Hapalemur_griseus"      
#>  [7] "Gorilla_gorilla"         "Pan_troglodytes"        
#>  [9] "Phaner_furcifer"         "Pithecia_pithecia"      
#> [11] "Pongo_pygmaeus"          "Tarsius_tarsier"        
#> [13] "Trachypithecus_vetulus"  "Varecia_variegata"

La función para análisis PGLS en caper es “pgls”. Para ajustar un modelo del logaritmo de la longitud de gestación contra el logaritmo de la masa corporal que utilice la estimación de máxima verosimilitud de \(\lambda\) utilizamos el siguiente código:

model.pgls<-pgls(log(GestationLen_d) ~ log(AdultBodyMass_g),
                 data = primate, lambda ='ML')

summary(model.pgls)
#> 
#> Call:
#> pgls(formula = log(GestationLen_d) ~ log(AdultBodyMass_g), data = primate, 
#>     lambda = "ML")
#> 
#> Residuals:
#>       Min        1Q    Median        3Q       Max 
#> -0.098899 -0.011661  0.003082  0.017758  0.075133 
#> 
#> Branch length transformations:
#> 
#> kappa  [Fix]  : 1.000
#> lambda [ ML]  : 0.892
#>    lower bound : 0.000, p = 1.1435e-14
#>    upper bound : 1.000, p = 0.00046393
#>    95.0% CI   : (0.753, 0.967)
#> delta  [Fix]  : 1.000
#> 
#> Coefficients:
#>                      Estimate Std. Error t value  Pr(>|t|)    
#> (Intercept)          4.290229   0.160355 26.7546 < 2.2e-16 ***
#> log(AdultBodyMass_g) 0.104864   0.019628  5.3426 9.479e-07 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> Residual standard error: 0.0261 on 75 degrees of freedom
#> Multiple R-squared: 0.2757,  Adjusted R-squared: 0.266 
#> F-statistic: 28.54 on 1 and 75 DF,  p-value: 9.479e-07

Si obtienen un error relacionado con la optimización (ver nota del error debajo), puede que necesiten ajustar los límites en la búsqueda del espacio de máxima verosimilitud. A cada uno de los parámetros que se pueden ajustar se les asignan límites, los valores por defecto para los límites son \(\lambda\): 1e-6 a 1, \(\kappa\): 1e-6 a 3 y \(\delta\): 1e-6 a 3. Para cambiar los límites, utilice el argumento “bounds” dentro de su modelo pgls, por ejemplo:

# model.pgls<-pgls(log(GestationLen_d)~ log(AdultBodyMass_g),
#                  data = primate,
#                  lambda ="ML",
#                  bounds = list(lambda=c(0.001,1),
#                                kappa=c(1e-6,3),
#                                delta=c(1e-6,3)))

Sección de código comentada dió error:
Error in pgls(log(GestationLen_d) ~ log(AdultBodyMass_g), data = primate, : Problem with optim:52ERROR: ABNORMAL_TERMINATION_IN_LNSRCH

Tengan en cuenta que hemos cambiado el límite inferior de lambda de 1e-6 a 0.001. Si desean cambiar los límites de uno de los parámetros, deben declarar también los límites de cada uno de los demás parámetros. Por lo tanto, arriba indicamos los límites para \(\delta\) y \(\kappa\) aunque no estamos ajustando estos parámetros con ML.

La salida debería ser la siguiente:

Call:
pgls(formula = log(GestationLen_d) ~ log(AdultBodyMass_g), data = primate, 
    lambda = "ML")
Residuals:
      Min        1Q    Median        3Q       Max 
-0.098899 -0.011661  0.003082  0.017758  0.075133 
Branch length transformations:
kappa  [Fix]  : 1.000
lambda [ ML]  : 0.892
   lower bound : 0.000, p = 1.1435e-14
   upper bound : 1.000, p = 0.00046393
   95.0% CI   : (0.753, 0.967)
delta  [Fix]  : 1.000
Coefficients:
                     Estimate Std. Error t value  Pr(>|t|)    
(Intercept)          4.290229   0.160355 26.7546 < 2.2e-16 ***
log(AdultBodyMass_g) 0.104864   0.019628  5.3426 9.479e-07 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 
Residual standard error: 0.0261 on 75 degrees of freedom
Multiple R-squared: 0.2757,     Adjusted R-squared: 0.266 
F-statistic: 28.54 on 2 and 75 DF,  p-value: 6.062e-10 

Además de los resultados de regresión estándar, la salida incluye el valor ML estimado de \(\lambda\) (0.892) y los valores \(p\) de las pruebas de razón de verosimilitud que muestran si la \(\lambda\) ML es significativamente diferente de 0 o 1 (véase 5.3: Estimación de la señal filogenética - \(\lambda\) de Pagel). En este caso, la estimación de máxima verosimilitud de \(\lambda\) es significativamente diferente tanto de cero como de 1; en otras palabras, hay evidencia de señal filogenética, pero el rasgo no ha evolucionado como se esperaba según un modelo de evolución de movimiento browniano (dadas las longitudes de rama utilizadas).

También podemos representar gráficamente los resultados, incluido un gráfico de la superficie de probabilidad de \(\lambda\) (que merece la pena investigar: a veces se encuentran patrones inesperados, ¡con picos en 0 o 1!):

plot(log(GestationLen_d)~ log(AdultBodyMass_g), data = primatedata)
abline(model.pgls)
profile_lambda=pgls.profile(model.pgls, which="lambda") # vary lambda
plot(profile_lambda)

Figura 10. Izquierda: Relación lineal entre el tamaño al nacer en función del tamaño corporal del adulto. Derecha: Distribución de probabilidad de valores de lambda óptimos para el modelo PGLS ajsutado

Los parámetros \(\kappa\) y \(\delta\) son transformaciones adicionales de longitud de rama implementadas en caper que pueden mejorar el ajuste de los datos al árbol (Véase 7.3: Transformaciones comunes de longitud de rama para más detalles). Para optimizar \(\kappa\) o \(\delta\), utilicen los argumentos delta = “ML” o kappa = “ML” en lugar de lambda = “ML” en el código anterior. Consideren que no es aconsejable optimizar más de uno de estos parámetros al mismo tiempo, ya que sería casi imposible interpretar los resultados. Además consideren que para \(\kappa\), necesitará una estructura de varianza-covarianza diferente de la utilizada anteriormente, concretamente un 3D vcv.array. Para ello, utilice el argumento “vcv.dim” al construir su objeto de datos comparativo:

primate_3d <- comparative.data(phy = primatetree,
                               data = primatedata,
                               names.col = Binomial,
                               vcv.dim=3, vcv = TRUE,
                               na.omit = FALSE,
                               warn.dropped = TRUE)


model.pgls.kappa<-pgls(log(GestationLen_d)~ log(AdultBodyMass_g),
                       data = primate_3d, kappa = "ML")

Debería encontrar que la estimación de máxima verosimilitud de \(\kappa\) es 0.458, y que también es significativamente diferente de cero (como se predijo bajo un modelo “especiacional” de evolución) y diferente de uno (como se predijo bajo el movimiento browniano utilizando las longitudes de rama dadas). Veamos la superficie de verosimilitud:

profile_kappa=pgls.profile(model.pgls.kappa, which="kappa") # vary kappa
plot(profile_kappa)

Figura 11. Distribución de probabilidad de valores de kappa óptimos (de máxima verosimilitud) para el modelo PGLS ajsutado

Otras consideraciones: Valores atípicos (“outliers”) con caper

Los valores atípicos pueden afectar gravemente a las estimaciones de los parámetros de cualquier modelo de regresión, y los modelos PGLS no son una excepción. A veces recomiendan la eliminación de valores atípicos con residuos “estudentizados” \(>\pm3\). El paquete caper no tiene actualmente una opción automatizada de eliminación de valores atípicos (compárese con “crunch”, sección 7.2.2), pero puede que la tenga en un futuro próximo. Comprueben las actualizaciones en ?pgls.

Para identificar los valores atípicos en un modelo PGLS, primero hay que extraer los residuos filogenéticos del modelo. A continuación, hay que estandarizar estos residuos dividiéndolos por la raíz cuadrada de su varianza. Por último, los residuos pueden compararse con los nombres de las especies e identificar el valor atípico (en el ejemplo debajo, Colobus_polykomos). A continuación, puede eliminar estas especies de los datos y del árbol, y volver a realizar el análisis. Consideren que puede que tenga que seguir eliminando especies hasta que no haya más valores atípicos.

res<- residuals(model.pgls, phylo = TRUE)
res<- res/sqrt(var(res))[1]
rownames(res)<-rownames(model.pgls$residuals)
rownames(res)[(abs(res)>3)]
#> [1] "Colobus_polykomos"

primate_nooutliers<-primate[-which(abs(res)>3),]
model.pgls_nooutliers<-pgls(log(GestationLen_d)~ log(AdultBodyMass_g), data = primate_nooutliers, lambda='ML')
summary(model.pgls_nooutliers)
#> 
#> Call:
#> pgls(formula = log(GestationLen_d) ~ log(AdultBodyMass_g), data = primate_nooutliers, 
#>     lambda = "ML")
#> 
#> Residuals:
#>       Min        1Q    Median        3Q       Max 
#> -0.056389 -0.013711 -0.002223  0.014340  0.085142 
#> 
#> Branch length transformations:
#> 
#> kappa  [Fix]  : 1.000
#> lambda [ ML]  : 0.889
#>    lower bound : 0.000, p = 2.1427e-14
#>    upper bound : 1.000, p = 0.00038538
#>    95.0% CI   : (0.748, 0.965)
#> delta  [Fix]  : 1.000
#> 
#> Coefficients:
#>                      Estimate Std. Error t value  Pr(>|t|)    
#> (Intercept)          4.289232   0.160836 26.6684 < 2.2e-16 ***
#> log(AdultBodyMass_g) 0.104958   0.019702  5.3272 1.034e-06 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> Residual standard error: 0.02616 on 74 degrees of freedom
#> Multiple R-squared: 0.2772,  Adjusted R-squared: 0.2674 
#> F-statistic: 28.38 on 1 and 74 DF,  p-value: 1.034e-06
par(mfrow=c(2,2))
plot(model.pgls_nooutliers)
abline(model.pgls_nooutliers)

Figura 12. Diagnóstico gráfico del modelo PGLS ajsutado

NOTA: Estos son residuos filogenéticos, por lo que detectan “valores atípicos filogenéticos”. Por lo tanto, los valores atípicos pueden no parecerlo en un gráfico de los datos brutos.

Gráficos de diagnóstico del modelo

Por lo general, merece la pena comprobar los gráficos de diagnóstico del modelo siempre que ajuste un modelo en R para comprobar que sus datos cumplen los supuestos del modelado lineal. El método es el mismo para los modelos OLS, contrastes independientes y PGLS (aunque los gráficos son ligeramente diferentes):

par(mfrow=c(2,2)) # permite ver los 4 gráficos juntos
plot(model.pgls)

Figura 13. GFáficos de diagnóstico de los modelos PGLS

# Para establecer la ventana gráfica a un solo gráfico
# par(mfrow=c(1,1))
# plot(model.pgls, which=1)

Básicamente, lo que se busca en estos gráficos es lo siguiente

1. Ningún residuo del modelo “estudentizados” \(>\pm 3\). Debe eliminarse cualquier especie con residuos tan grandes. Estos valores atípicos pueden influir excesivamente en los resultados de la regresión.

2. Los puntos del gráfico Q-Q deben formar aproximadamente una línea recta (en lugar de una forma de plátano).

3 y 4) Deben mostrar una dispersión bastante aleatoria de los puntos. Hay que evitar los patrones claros.

Se necesita práctica para saber lo que es “bueno”, “malo” y “aceptable” con estos gráficos. Estos gráficos tienen buen aspecto, pero parece que hay varios valores atípicos; puede que merezca la pena investigar cómo afectan estos valores a los resultados.

PGLS en modelos generalizados

OBSERVACIÓN: Esta sección está en construcción. Si tienen urgencia en implementar este tipo de modelos, por favor referirse a wzmli/phyloglmm: First release of phylogenetic comparative analysis in lme4-verse, ver haciendo click y avisarme si lo requieren.

Los ejemplos anteriores mostraron cómo ajustar modelos de regresión con una matriz de covarianza filogenética para relaciones lineales entre los rasgos (o linearizaciones desp{es de transformar la/s variable/s con comportamientos no lineales}). Es más interesante tener la posibilidad de combinar PGLS con la funcionalidad de paquetes que ajusten modelos lineales mixtos generalizados (GLMMs). En este ejemplo veremos cómo implementar PGLS con glmmTMB (paquete de Ben Bolker). glmmTMB es un paquete de R para ajustar GLMMs con una amplia gama de distribuciones estadísticas (Gaussiana, Poisson, binomial, binomial negativa, Beta…), así como extensiones para controlar problemas como la inflación cero (“zero-inflation”), la heteroscedasticidad y la autocorrelación. El ejemplo, tomado de Michael Li and Ben Bolker, reproduce los casos del capítulo 11 de Garamszegi (2014) utilizando GLMM filogenéticos basados en lme4 y glmmTMB.

Para el ajuste necesitamos un efecto aleatorio en la fórmula que incluya un término (…|phylo), para construir la estructura básica de efectos aleatorios multiplicada por la matriz phyloZ (véase la viñeta para más detalles sobre la matriz phyloZ).

# install.packages('glmmTMB')
library(ape)
library(Matrix)
library(lme4)
library(MASS)
library(glmmTMB)
library(coda)
library(lattice)
library(broom)
library(dplyr)

source("lme4_phylo_setup.R")
source("glmmTMB_phylo_setup.R")

Obtenención de datos

Del chapter 11 of @garamszegi2014modern: los datos están acá

if (!file.exists("data/phylo.nex")) {
    dir.create("data")
    download.file("http://mpcm-evolution.com/OPM/Chapter11_OPM/data.zip",
                  dest="data/OPM_ch11_data.zip")
    setwd("data")
    untar("OPM_ch11_data.zip")
    setwd("..")
}
phylo <- ape::read.nexus("data/phylo.nex")

Calculen la matriz \(Z\) apropiada por adelantado, para medir la velocidad (igualmente la necesitaremos más abajo):

system.time(phyloZ <- phylo.to.Z(phylo))
#>    user  system elapsed 
#>    0.37    0.00    0.39

Comparación de resultados con el ejemplo gaussiano del capítulo 11 de Garamszegi (2014)

datG <- read.table("data/data_simple.txt",header=TRUE)
## create observation-level grouping variable
datG$obs <- factor(seq(nrow(datG)))
datG$sp <- datG$phylo
# phylo_lmm_fit <- phylo_lmm(phen~cofactor+(1|sp)
#   , data=datG
#   , phylonm = "sp"
#   , phylo = phylo
#   , phyloZ=phyloZ
#   , REML = TRUE
#   , control=lmerControl(check.nobs.vs.nlev="ignore",check.nobs.vs.nRE="ignore")
# )
# print(summary(phylo_lmm_fit))

Sección de código comentada dió error:
Error in match.fun(FUN) : object ‘safeDeparse’ not found

Del mismo modo, el ajuste mediante glmmTMB:

# glmmTMB_fit <- glmmTMBphylo(phen~cofactor+(1|sp)
#     , data=datG
#     , phyloZ=phyloZ
#     , phylonm = "sp"
#     , doFit=TRUE
#     , dispformula = ~1
#     , REML = FALSE
#   ) 
# print(summary(glmmTMB_fit))

Sección de código comentada dió error:
Error in model.matrix(glmmTMB:::drop.special2(fixedform), fr, contrasts) : object ‘drop.special2’ not found

Comparación de resultados con el ejemplo gaussiano de mediciones repetidas del capítulo 11

datR <- read.table("data/data_repeat.txt",header=TRUE)
datR$obs <- factor(seq(nrow(datR)))
datR <- (datR 
  %>% mutate(sp = species
        , animals = phylo
        )
)
datR$spec_mean_cf <- sapply(split(datR$cofactor,datR$phylo),mean)[datR$phylo]
datR$within_spec_cf <- datR$cofactor-datR$spec_mean_cf
# phylo_lmm_fit <- phylo_lmm(phen~spec_mean_cf+within_spec_cf+(1|sp) + (1|animals)
#   , data=datR
#   , phylonm = "sp"
#   , phylo = phylo
#   , phyloZ=phyloZ
#   , REML = FALSE
#   , control=lmerControl(check.nobs.vs.nlev="ignore",check.nobs.vs.nRE="ignore")
# )
             
# print(summary(phylo_lmm_fit))

# glmmTMB_fit2 <- glmmTMBphylo(phen~spec_mean_cf+within_spec_cf+(1|sp) + (1|animals)
#     , data=datR
#     , phyloZ=phyloZ
#     , phylonm = "sp"
#     , doFit=TRUE
#     , dispformula = ~1
#     , REML = FALSE
#   ) 
# print(summary(glmmTMB_fit2))

Sección de código comentada dió error:
Error in model.matrix(glmmTMB:::drop.special2(fixedform), fr, contrasts) : object ‘drop.special2’ not found

Comparación de resultados con el ejemplo no gaussiano del capítulo 11

dat <- read.table("data/data_pois.txt",header=TRUE)
dat$obs <- factor(seq(nrow(dat)))


dat <- dat %>% mutate(sp=phylo) 
# phylo_glmm_fit <- phylo_glmm(phen_pois~cofactor+(1|sp)+(1|obs)
#   , data=dat
#   , phylonm = "sp"
#   , family = poisson
#   , phylo = phylo
#   , phyloZ=phyloZ
#   , control=lmerControl(check.nobs.vs.nlev="ignore",check.nobs.vs.nRE="ignore")
# )
# 
# summary(phylo_glmm_fit)

# glmmTMB_fit3 <- glmmTMBphylo(phen_pois~cofactor+(1|sp)+(1|obs)
#     , data=dat
#     , family = poisson
#     , phyloZ=phyloZ
#     , phylonm = "sp"
#     , doFit=TRUE
#     , dispformula = ~1
#     , REML = FALSE
#   ) 
# print(summary(glmmTMB_fit3))

Sección de código comentada dió error:
Error in model.matrix(glmmTMB:::drop.special2(fixedform), fr, contrasts) : object ‘drop.special2’ not found