Crecimiento logístico (denso-dependiente) continuo. Aplicación en R (paquete deSolve y primer)

Crecimiento logístico (denso-dependiente) continuo

#### CRECIMIENTO DENSO-DEPENDIENTE CONTINUO
#la versión clásica es la de la ecuación de crecimiento logístico continuo: dN/dt?rN*((K-N)/K)=rN(1-alfa*N)
#cómo relajar los supuestos y generalizar la ecuación: 1) r la tasa de crecimiento per cápita DI es realmente la tasa neta r=b-d. El efecto DI positivo de b en la tasa de crecimiento es balanceado por el efecto DD negativo de b. A su vez d incluye un pequeño efecto DD positivo que aumenta la tasa de crecimiento porque reduce el efecto negativo de alfa^N.
#DD en la tasa de mortalidad: representamos la DD (1-alfa*N) como la función más simple F(N)=1-alfa*N, donde alfa es un efecto neto. Podemos generalizar: dN/N*dt=r*N^F(N), donde F(N)=B(N)+D(N). Aquí, podemos modelar D(N)=g*N^2 (g es una constante), donde la densidad no tiene efecto sobre la tasa de mortalidad cuando N es bajo porque no existe ningún factor limitante, pero cuando N es grande beneficia paradójicamente a la tasa de crecimiento porque la libera de los recursos limitantes.
#DD en la tasa de natalidad: efecto Allee, donde B(N)?-a*N^2+e*N-f, donde a, e y f son constantes.
#DD en las tasas de mortalidad y natalidad
 B.N <- expression(-a * N^2 + e * N - f); D.N <- expression(g * N^2)
 a <- 1/1600; e <- 80 * a; f <- 0.2; g <- 4e-05; N <- 0:100
 plot(N, eval(B.N), type = "l", ylab = "B(N), D(N)");     lines(N, eval(D.N), lty = 2);     abline(h = 0, lty = 3);     legend("bottom", c("Effect on Birth Rate", "Effect on Death Rate"), lty = 1:2, bty = "n")
 plot(N, eval(B.N) + eval(D.N), type = "l", ylim = c(-0.4, 1), ylab = "Net Density Dependence, F(N)");     abline(h = 0, lty = 3);     curve(1 - 0.01 * x, 0, 100, lty = 2, add = T);     legend("topright", c("Generalized", "Logistic"), lty = 1:2, bty = "n")
#la extrema flexibilidad implica que B(N) y D(N) puede tomar infinitas formas mientras que su suma sea 1-alfa*N. Sin información adicional podemos utilizar las funciones más simples para B y D: dN/N*dt=(b-mb*N)-(d-md*N), donde mb y md son constantes que dsecriben el efecto de un individuo adicional en las tasas de nacimiento y mortalidad DI per cápita. También podemos simplificar aún más utilizando mb+md=m, b-d=r, entonces: dN/N*dt=r-m*N, entonces el efecto neto de N es lineal, aún cuando los efectos de N en la natalidad y mortalidad no lo son. Podemos también crear otra constrante alfa=m/r donde regresamos a dN/N*Dt=r*(1-alfa*N).
#La integral del crecimiento logístico: si integramos la ecuación dN/N*dt=r*(1-alfa*N) obtenemos Nt=No*exp(rt)/(1+alfa*No*(exp(rt)-1))
#equilibrio en el modelo de crecimiento logístico continuo: los valores de N para los cuales dN/dt=0. Son N*=0, N*=1/alfa.
#dinámica alrededor del equilibrio (estabilidad):: es una medida de cuánto tiende a estar igual, dado disturbios externos o cambios en el estado del sistema.El análisis de estabilidad lineal calcula el grado de atracción o repulsión.  Realizamos un análisis analítico de la estabilidad lineal: consideramos el gráfico de la tasa de crecimiento dN/dt vs. N.
 pop.growth.rate <- expression(r * N * (1 - alpha * N)); r <- 1; alpha <- 0.01; N <- 0:120
 plot(N, eval(pop.growth.rate), type = "l", ylab = "Population Growth Rate (dN/dt)", xlab = "N"); abline(h = 0); legend("topright", "r=1", lty = 1)
 N <- c(0, 10, 50, 100, 115);  points(N, eval(pop.growth.rate), cex = 1.5); text(N, eval(pop.growth.rate), letters[1:5], adj = c(0.5, 2)); arrows(20, 2, 80, 2, length = 0.1, lwd = 3); arrows(122, -2, 109, -2, length = 0.1, lwd = 3)
 #para calcular la pendiente de la función necesitamos las derivadas parciales de la tasa de crecimiento DN/dt respecto a N: N`=r*N*(1-alfa*N) donde N`indica la derivada temporal. dado N`su derivada parcial respecto a N es deltaN`/deltaN=r-2*r*alfa*N. Estamos interesados en el equilibrio donde N=1/alfa, por lo que lo sustituimos deltaN`/dN=-r y así la pendiente negativa en el equilibrio demuestra la estabilidad del equilibrio. Cerca de N* la tasa de crecimiento es -r*x, entonces dx/dt=-r*x, integramos respecto al tiempo y obtenemos xt=xo*exp(-rt), donde xo es el tamaño del desplazamiento inicial en relación el equilibrio y xt es el tamaño del desplazamiento en el tiempo t en relación al equilibrio. Podemos obtener el tiempo de retorno característico (tiempo requerido para alcanzar xt) tomando xt?xo/e y definiéndolo como la cantidad de tiempo necesario para reducir la perturbación en un 63% (i.e. 1/e): entonces xo/e=xo*edxp(-rt) y así t=1/r.
 #diferenciación simbólica para obtener las derivadas y evaluar los dos equilibrios:
 dF.dN <- deriv(pop.growth.rate, "N"); N <- c(0, 1/alpha); eval(dF.dN) #The first value, 1, corresponds to the first value of N, which is 0. Because it is positive, this indicates that the perturbation will increase with time, meaning that N = 0 is a repellor. The second value, -1, is negative, and so indicates that the perturbation will decrease with time, meaning that N = 1/
is an attractor.
 #dinámica: ODE
 clogistic <- function(times, y, parms) { n <- y[1]; r <- parms[1]; alpha <- parms[2]; dN.dt <- r * n * (1 - alpha * n); return(list(c(dN.dt))) }; prms<- c(r = 1, alpha = 0.01); init.N <- c(1); t.s <- seq(0.1,10, by = 0.1)
 library(deSolve); out <- ode(y = init.N, times = t.s, clogistic, parms = prms); plot(out[, 1], out[, 2], type = "l", xlab = "Time", ylab = "N") #la salida es una matriz que incluye los tiempos y N
 #DD retrasada o crecimiento logístico con tiempo de retardo: dN/dt=r*N*(1-alfa*Nt-theta), donde theta es el grado de retardo. Ahora la dinámica se vuelve bien complicada, como en el modelo discreto, dado un lag y r suficiente.
 #graficando poblaciones aleatorias: ceramos 20 poblaciones con diferentes características
 outmat <- matrix(NA, nrow = length(t.s), ncol = 20); for (j in 1:20) outmat[, j] <- { y <- runif(n = 1, min = 0, max = 120); prms <- c(r = runif(1, 0.01, 2), alpha = 0.01); ode(y, times = t.s, clogistic, prms)[, 2] }; matplot(t.s, outmat, type = "l", col = 1, ylab = "All Populations")
 #otras formas de DD (no logística): modelo Richards, modelo Gompertz, modelo de Ricker y el modelo theta-logístico que agrega un parámetro para aumentar la flexibilidad y la generalidad: dN/dt=r*N*log(1-(alfa*N)^theta), donde theta es estrictamente positivo y theta=0 es que no crece. Aplicaremos la función Theta-logistic: y la DD theta-logística:
 thetalogistic <- function(times, y, parms) { n <- y[1]; with(as.list(parms), { dN.dt <- r * n * (1 - (alpha * n)^theta); return(list(c(dN.dt))) }) }
 r <- 0.75; alpha <- 0.01; theta <- c(0.5, 1, 2); N <- 0:110; theta.out <- sapply(theta, function(th) { 1 - (alpha * N)^th }); matplot(N, theta.out, type = "l", col = 1); abline(h = 0); legend("topright", legend = paste("theta =", c(2, 1, 0.5)), lty = 3:1) #DD theta-logístico
 thetaGR.out <- sapply(theta, function(th) { r * N * (1 - (alpha * N)^th) }); matplot(N, thetaGR.out, type = "l", col = 1); abline(h = 0); legend("bottomleft", legend = paste("theta =", c(2, 1, 0.5)), lty = 3:1) #tasa de crecimiento theta-logístico
 prms <- c(r <- 0.75, alpha <- 0.01, theta = 1); thetaN <- sapply(theta, function(th) { prms["theta"] <- th; ode(y = 1, t.s, thetalogistic, prms)[, 2] }); matplot(t.s, thetaN, type = "l"); legend("topleft", legend = paste("theta =", c(2, 1, 0.5)), lty = 3:1, bty = "n") #dinámica theta-logístico
 #ajustando los modelos a los datos: e.g. examinar los efectos de los nutrientes en la interacción de un alga dentro de una red trófica (estudiamos el rol de los recursos en la alteración de las interacciones poblacionales dentro de una red trófica simple)
    #agregar recursos influye en la fuerza de la DD per cápita alfa
    #la variedad de respuestas observadas a la disponibilidad de recursos es mediada pro la configuración de a red trófica en al cual la población se encuentra.
    #estudiamos los efectos de la tasa de suplemento de recursos en los parámetros de DD lineal (crecimiento logístico) r y alfa de Clostrium acerosum, dado interacciones competitivas con una bacteria y predación por Colpidium striatum.
    #1) importamos y examinamos los datos, 2) miramos la dinámica poblacional para cada microcosmo replicadoo, 3)ver outliers debudos al muestreo, 4)estimamos parámetros para una de las poblaciones, 5)ajustamos modelos logísticos para cada pobalción, 6)ajustamos un modelo para todas las poblacionas y probamos los efectos de los nutrientes en los parámetros de crecimiento logístico
    library(nlme); library(lattice); data(ClostExp); summary(ClostExp)
    xyplot(No.per.ml ~ Day | Nutrients, ClostExp, groups = rep, type = "b", scales = list(relation = "free"), auto.key = list(columns = 4, lines = T))
    subset(ClostExp, Nutrients == "high" & No.per.ml > 1000); subset(ClostExp, Nutrients == "low" & No.per.ml > 100)
    Hi.c <- subset(ClostExp, Nutrients == "high" & rep == "c") #para estimar los parámetros sin considerar la dinámica temporal calculamos la tasa de crecimiento per cápita desde los datos y ajustamos una regresión lineal vs. N (dN/N*dt=r-r*alfa*N), el intercepto con y es la estimación de r y la pendiente una estimación de r*alfa. Estimamos la tasa de crecimiento per cápita mediante pgr=log(Nt+1/Nt)/i, donde i es el tiempo de intervalo entre dos densidades poblacionales.
    n <- nrow(Hi.c);  N.change <- Hi.c$No.per.ml[-1]/Hi.c$No.per.ml[-n];    interval <- diff(Hi.c$Day);       pgr <- log(N.change)/interval
    Nt <- Hi.c$No.per.ml[-n];       plot(pgr ~ Nt);       mod1 <- lm(pgr ~ Nt);       abline(mod1); summary(mod1)
    EachPop <- lapply(split(ClostExp, list(ClostExp$Nutrients, ClostExp$rep)), function(X) { n <- nrow(X); N.change <- (X$No.per.ml[-1]/X$No.per.ml[-n]); interval <- diff(X$Day); data.frame(Nutrients = as.factor(X$Nutrients[-n]), rep = as.factor(X$rep[-n]), pgr = log(N.change)/interval, Nt = X$No.per.ml[-n]) }) #estimamos los parámetros con todas las replicaciones
    AllPops <- NULL;
 for (i in 1:length(EachPop)) { AllPops <- rbind(AllPops, EachPop[[i]])};
 xyplot(pgr ~ Nt | rep * Nutrients, AllPops, layout = c(4, 2, 1), scales = list(x = list(rot = 90)), panel = function(x, y) { panel.grid(h = -1, v = -4); panel.xyplot(x, y, type = c("p", "r")) }) #ajustamos con todos los individuos en uno
    AllPops$ID <- with(AllPops, Nutrients:rep); modSlope <- lme(pgr ~ Nt + Nutrients + Nt:Nutrients, data = AllPops, random = ~1 | ID) #modelo estadístico para el crecimiento per cápita. Un modelo mixto incluye efectos aleatorios (e.g. bloques o subjects = microcosmo) y efectos fijos (e.g. tratamiento). Poner etiqueda ID para cada microcosmos.
    xyplot(resid(modSlope) ~ fitted(modSlope)); qqmath(~resid(modSlope) | ID, data = AllPops, layout = c(8, 1, 1)); qqmath(~ranef(modSlope)) #chequeamos los supuestos del modelo: que el ruido (residuales) y los efectos aleatorios tienen distribución normal. Observamos grandes residuos para grandes valores ajustados (esto no es bueno). el gráfico cuantil-cuantil de los res y efectos aleat se ajustan bastante a la linea diagonal. Los residuales para cada microcosmo parecen normales.
    modSlope2 <- update(modSlope, weights = varExp()); anova(modSlope, modSlope2); anova(modSlope2) #having more observations at larger values, but is also a typical and common problem — variables with small absolute values simple cannot vary as much as variables with large absolute value. Therefore, we specify that that model takes that into account by weighting residuals observations differently, depending upon the expected (i.e. fitted or mean) value. Observamos que el modelo más complicado (7df) es significativamente mejor (likelihood ratio=7.1, P<0.05) y más parsimonioso (menor AIC). Evaluamos los efectos fijos (nutrientes y tamaño poblacional) con el análisis de varianza final. Vemos que el efecto de Nt depende del nivel de nutrientes (Nt:Nutrients, P<0.0001). Dado esto debemos no realizar conclusiones bados en los efectos principales delos nutrientes.
    summary(modSlope2)$tTable; cfs <- fixef(modSlope2); -(cfs[2] + cfs[4])/(cfs[1] + cfs[3]) #estimamos interceptos y pendientes asociadas a distintos tratamientos. Este constraste por default (parámetros del modelo lineal) indica que el intercepto (estimación de r) para el grupo de mayores nutrientes es ~.13 y la pendiente para el microcosmo con altos nutrientes es -.000286; mientras para bajos niveles de nutrientes tenemos .13-.07 de intercepto (no sig. respecto al otro microcosmo) y pendiente -.000278-.00573 (sí sig. respecto al otro microcosmo). Los últimos calculso nos dicen que los efectos per cápita son mucho menores en el tratamiento de altos nutrientes.
    ilogistic <- function(t, alpha, N0, r) { N0 * exp(r * t)/(1 + alpha * N0 * (exp(r * t) - 1)) };
 x11(); plot(No.per.ml ~ Day, ClostExp, subset = Nutrients == "high"); curve(ilogistic(x, alpha = -cfs[2]/cfs[1], N0 = 6, r = cfs[1]), 1, 60, add = T, lty = 2)  ; Cmod.list <- nlsList(No.per.ml ~ ilogistic(Day, alpha, N0, r) | ID, data = ClostExp, start = c(alpha = 1/100, N0 = 5, r = 0.15), control = list(maxiter = 1000));
 x11();plot(coef(Cmod.list), layout = c(3, 1, 1));
 x11(); plot(Cmod.list)  #aproximación tiempo exlpícita: 1) creamos y evaluamos el modelo logístico no-lineal, 2) ajustamos el modelo a cada microcosmo y 3) ajustamos un modelo mizto no-lineal que pruebe qué parámetro de crecimiento logístico varía entre tratamientos. We note that models cannot be fit for three microcosms. Probamos los supuestos
    Cmod.all <- nlme(Cmod.list, random = pdDiag(form = alpha + N0 + r ~ 1), weights = varPower()) #creamos un modelo mixto no-lineal simple. utilizamos pdDiag para especificar que cada microcomo tiene su efecto único propio en cada parámetro, fuera de las concentraciones de nutrientes. Con weights especificamos la varianza de los residuales como una función de potencia (a^2t) de la media y que ajusta un parámetro extra t.
    cfs2 <- fixef(Cmod.all) #vemos que los nutrientes causaron grandes diferencias en las pendientes de pgrvs.N y estimamos el alpha para cada nivel de nutrientes.
    Cmod.all2 <- update(Cmod.all, fixed = list(alpha ~ Nutrients, N0 + r ~ 1), start = c(cfs2[1], 0, cfs2[2], cfs2[3])); summary(Cmod.all2) #actualizamos el modelo, especificando que alpha depende del nivel de nutrientes (estimamos cada alpha para cada nivel de nutrientes). Especificamos que los demás parámetros son los mismos para ambos tratamientos. Graficamos la diagnosis nuevamente. Observamos que la varianza de alpha o r son muy pequeños, por lo que los removemos del análisis y compararemos los modelos.
    Cmod.all3 <- update(Cmod.all2, random = N0 ~ 1); anova(Cmod.all2, Cmod.all3)  #vemos que el modelo menos complejo no es signif. diferente (P>0.05) y es más parsimonioso (menor AIC).
    intervals(Cmod.all3)$fixed  #examinamos los IC para tener las estimaciones de los coef de competemncia intraespecífica aii
    1/intervals(Cmod.all3)$fixed[1:2, c(1, 3)] #calculamos los límites de los K basados en nuestros IC inf y sup en alpha.
    X11(); plot(augPred(Cmod.all3, level = 0:1), layout = c(4, 2, 1), scales = list(y = list(relation = "free"))) #graficamos los efectos fijos y la variazión agregada debido a la variación aleatoria en n0 y r debido al ID del microcosmos.
    #para evaluar la disparidad de los modelos con los datos debemos considerar: error de muestreo, diseño experimental, dinámica.

Created by Pretty R at inside-R.org


Para descargar: https://www.dropbox.com/sh/qqnr01l8s3wglqm/ZJQBuCVJmk