La PRCF ¿qué es y cómo se calcula?

Recientemente me he puesto a analizar unos datos de series temporales de dinámica de poblaciones de procesionaria (estimadas a partir de mediciones del grado de infestación de rodales) en Mora de Rubielos, Teruel, en el marco de una colaboración a tres bandas entre el Laboratorio de Sanidad Vegetal de Mora de Rubielos (Rodólfo Sánchez, Gerárdo Sánchez), la Universidad de Granada (José Antonio Hódar, Regino Zamora) y la Universidad Rey Juan Carlos (en mi humilde persona).

Para poder discernir entre factores endógenos y exógenos que podrían estar causando las oscilaciones de las series, una de las primeras cosas que hay que hacer es tratar de determinar si hay dependencia de la densidad y con cuantos retardos. Es decir, si el hecho de que, por ejemplo, uno, dos o tres años atrás, la población estuviera en un máximo, va a influir en que ahora, en tiempo t, la población haya sufrido un descenso en el número de individuos como consecuencia del agotamiento de los recursos u otros factores que podrían estar operando a nivel intraespecífico. Pues bien, para determinar con cuantos retardos la propia variable tiene un efecto sobre la observación de esta misma variable en tiempo t, Alan Berryman (1999) sugiere el uso de lo que llama la PRCF (partial rate correlation function), que no traduzco porque no se si lo haría bien. Para aquellos que estéis un poco familiarizados con la jerga de series temporales, no os será desconocido la función de autocorrelación (ACF) y la función de autocorrelación parcial (PACF). Estas son herramientas exploratorias básicas para identificar la estructura de una serie temporal. Pues bien, en el caso de dinámica de poblaciones, Berryman sugiere que es más apropiado utilizar la PRCF... ¿y qué es esto de la PRCF? Se trata de correlacionar la tasa de cambio poblacional en tiempo t, con la densidad poblacional en tiempo t-1, t-2 y así sucesivamente. La tasa de cambio de la población es lo que Berryman denomina como variable Rt, y que en realidad es el log(Nt/Nt-1). Es decir, una medida de cómo la población cambia entre un periodo y el siguiente.

En mi periplo analítico, me encontré con que no tenía ni idea de cómo calcular la PRCF en R, que es ya el único software con el que trabajo. Busqué en los foros y nada. Lo más que encontré fue la página de un software (PAS) que diseñó Alan Berryman ya hace tiempo y que explicaba con fórmulas cómo calcular la PACF y, a partir de ahí, la PRCF. Pero honestamente, no entendí mucho. Así que di en mis indagaciones con un artículo de Mauricio Lima publicado en PLoS ONE donde decía explícitamente que ellos habían calculado la PRCF con un código escrito en R. Escribí a Mauricio y este rápidamente me dirigió a Sergio Estay, quien con mucha diligencia me pasó el código, que con permiso suyo y agradeciéndole de antemano su enorme ayuda, pongo aquí a disposición de todos aquellos interesados (he modificado un par de cosas para generalizarlo, pero en esencia, todo el código es suyo).

Código para calcular la función PRCF en R 

Al final, y viendo el código de Sergio (y con alguna aclaración suya) entendí que era esto de la PRCF... y es bastante simple. Para más de un retardo la PRCF es equivalente a la PACF. Y para un retardo es simplemente la correlación entre Rt (la tasa de cambio) y Nt-1. La (última) pregunta que yo me hice fue ¿por qué para dos o más retardos la PRCF es igual a la PACF? Y dándole vueltas mientras volvía a casa en coche, llegué a una conclusión, que expongo a continuación y que no se si será más o menos acertada (si alguien conoce una respuesta más acertada por favor que me corrija).

Básicamente lo que se está diciendo es que para dos retardos (y generalizable a otros retardos por encima de uno):

r(Rt ~ Nt-2) habiendo quitado el efecto de Nt-1 sobre Rt 
r(Nt ~ Nt-2) habiendo quitado el efecto de Nt-1 sobre Nt 

¿Tiene esto algún sentido? Lo tiene si pensamos que Rt se calcula como el log(Nt/Nt-1). Por tanto si quitamos el efecto de Nt-1 sobre Rt, lo que nos queda es básicamente Nt. Por tanto, con dos o más retardos, la PRCF será equivalente a la PACF.

Y ya para terminar, incluyo aquí los resultados de la PRCF para las 11 parcelas con las que estamos trabajando en Mora de Rubielos, y en donde se ve que la dinámica de la procesionaria está regulada por efectos de densidad con un único retardo.