Gráficas de líneas y puntos. Modelos

Libreria

ggplot2, dplyr

Datos

obdung.csv, son los datos observados de emergencia de escarabajos estercoleros de la especie Euniticellus intermedius (Cruz Rosales M. et al. 2012. Effect of ivermectin on the survival and fecundity of Euniticellus intermedius (Coleotera: Scarabaeidae). Rev. Biol. Trop. 60(1): 333-345)

La estructura de los datos es como sigue:

Se hicieron lecturas de la emergencia de adultos a traves del tiempo (DIA) sujetos a tres tratamientos. TEST, TEST-A e IV-B (éste contenia ivermectina); SEXO es el sexo de los individuos. PROP es la proporción de emergencia de los adultos emergidos por día en relación al total emergido.

exdungb.csv, archivo con los valores esperados (=predichos) por el modelo de emergencia de escarabajos.

obELP.csv, son los datos observados de las tasas de desarrollo de la palomilla gitana Cydia pomonella (Aghdam et al. 2009. Temperature-dependent development and temperature thresholds of codling moth. Environ. Entomol. 38(3): 885-895).

La estructura de los datos es como sigue:

Se hicieron lecturas del tiempo desarrollo de los diferentes estados de la palomilla (Stage), aqui solamente se analyzan larvas y huevo-larva, a diferentes temperaturas constantes (Temp), la respuesta fue la tasa promedio de desarrollo (dev, =1/dias).

predELP.csv

Son los valores predichos por diferentes modelos para los datos de C. pomonella.

Scripts

emergModel1.R

Genera valores predichos y construye tabla de datos.

devmodels.R

Contiene las funciones que representan los modelos de las tasas de desarrollo. Esta lista no es exhaustiva, el lector se le recomienda seleccionar su modelo de acuerdo a sus requerimientos, de ser necesario.

Los archivos se encuentran aquí:

https://sites.google.com/site/digitcognem/workfiles

Desarrollo

En esta sección vamos a cubrir el proceso de elaborar gráficas de puntos (observados) y el traslape con gráficas de linea que representan valores predichos (= esperados) por algún modelo matemático.

Los pasos generales para hacer un análisis de este tipo y las gráficas son los siguientes:

1. Seleccionar el modelo o modelos que pueden representar nuestros datos

2. Estimar los parámetros de estos modelos

3. Si el modelo es Ok, calcular los valores predichos

4. Arreglar los valores predichos en la misma estructura que los valores observados.

5. Realizar las gráficas.

1. Seleccionar el modelo o modelos que pueden representar nuestros datos

En el paso 1 se seleccionó el modelo de Kostal & Havelka (2001. Cryobiology. 42: 112-120):

PROP = 1/(1 + 10^(s0*(e50 - DIA)))

Donde PROP es la proporción de adultos emergidos, en el DIA establecido. Los parámetros del modelo son s0 y e50, el primero es la pendiente de la curva y el segundo es el tiempo mediano a la emergencia.

2. Estimación de parámetros.

En este caso particular, puesto que es un modelo no lineal, se empleó la función nls() de R. El archivo ecum.csv se filtró (ver Manejo de Datos) para cada uno de los tratamientos y sexo. Con esto, se empleó la función previa con valores iniciales obtenidos de la observación de los datos y los publicados por Kostal & Havelka.

Como ejemplo, se estiman los parámetros del modelo para datos del testigo A y machos:

A) Selección del conjunto de datos (tiene que activar la libreria dplyr):

tam<-filter(obdungb,TRAT=="TEST-A",SEXO=="MACHO")

Los registros de la tabla anterior son los siguientes:

DIA,TRAT,SEXO,INDIVIDUOS,PROP

34 TEST-A MACHO 0

36 TEST-A MACHO 0.04

38 TEST-A MACHO 0.07

41 TEST-A MACHO 0.22

43 TEST-A MACHO 0.43

45 TEST-A MACHO 0.78

48 TEST-A MACHO 0.93

50 TEST-A MACHO 0.95

52 TEST-A MACHO 0.97

55 TEST-A MACHO 1

57 TEST-A MACHO 1

Ahora se procede a aplicar la función nls() y estimar los parámetros del modelo, como ejemplo se utiliza el conjunto de datos anteriores (tam):

modtam<- nls(PROP ~ 1/(1 + 10^(s0*(e50 -DIA))),data= tam, start= list(s0=0.5, e50= 43), trace= TRUE)

0.05210755 : 0.5 43.0

0.01366298 : 0.3228733 43.0997889

0.006146556 : 0.2525288 43.2287525

0.006142391 : 0.2512933 43.2273048

0.00614224 : 0.2510583 43.2270763

0.006142234 : 0.2510131 43.2270331

0.006142234 : 0.2510044 43.2270248

0.006142234 : 0.2510027 43.2270232

Las instrucciones de nls son como sigue;

PROP ~ 1/(1 + 10^(s0*(e50 -DIA)))

Es la especificación del modelo en notación de R y las variables de la tabla de datos, en este caso data= tam que es el segundo argumento de la función.

start=list(s0=0.5,e50= 43)

Aquí se declaran los valores iniciales de los parámetros. El lector debe consultar trabajos similares para asignar valores iniciales, lo que facilita que la función encuentre los parámetros. En este caso, s0 se obtuvo observando el trabajo de Kostal & Havelka (op. cit.) y la mediana a la emergencia observando los datos.

trace=TRUE

hace que R imprima en la consola la secuencia iterativa de la búsqueda de los parámetros. La lineal final son los parámetros del modelo.

Si escribimos modtam en la consola, R imprime el modelo, sus parámetros y la suma de cuadrados de los residuales, observe que los parámetros son los mismos que los impresos anteriormente con el trace:

modtam

Nonlinear regression model

model: PROP ~ 1/(1 + 10^(s0 * (e50 - DIA)))

data: tam

s0 e50

0.251 43.227

residual sum-of-squares: 0.006142

Number of iterations to convergence: 7

Achieved convergence tolerance: 6.253e-06

Si escribimos summary(modtam), se presenta la estimación de los parámetros y el error estándar de cada uno de ellos, junto con una prueba de t, para la hipótesis nula par=0, es decir, que el valor del parámetro respectivo es cero, así como el valor de probabilidad de la prueba:

summary(modtam)

Formula: PROP ~ 1/(1 + 10^(s0 * (e50 - DIA)))

Parameters:

Estimate Std. Error t value Pr(>|t|)

s0 0.25100 0.01744 14.39 1.62e-07 ***

e50 43.22702 0.12340 350.30 < 2e-16 ***

---

Signif. codes:

0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.02612 on 9 degrees of freedom

Number of iterations to convergence: 7

Achieved convergence tolerance: 6.253e-06

De esta manera se calculan los parámetros de cada tratamiento y sexo, y se pueden registrar en el siguiente cuadro:

La última columna es el número de insectos original.

3. Si el modelo es Ok, calcular los valores predichos

En este paso, usualmente se calcular los valores predichos para el rango de valores de los observados, en este caso, son los días.

Para agilizar los cálculos, se puede escribir un script en R y escribir la secuencia de los instrucciones (script "emergModel1.R"):

# SCRIPT PARA CALCULAR VALORES PREDICHOS, EJEMPLO DE TEST-A MACHOS:

# escribir ecuación en lenguaje R

# En este caso es el modelo de Kostal & Havelka 2001

func.emerge<-function(dia, s0, e50){

result<- 1/(1+10^(s0*(e50-dia)));

return(result)

}

#Generar los valores de x, en el rango de los valores observados

DIA<-seq(36,57, by= 1)

#Generar los valores esperados aplicando la ecuación

# y con los parámetros estimados previamente

PROP<-func.emerge(DIA,0.251, 43.227)

# generar las etiquetas de las columnas de tratamiento y sexo

# La idea es tener una tabla con la estructura igual a la de

# los valores observados

# length(DIA) nos da el tamaño de la lista que contiene

# los valores predichos

TRAT<-rep("TEST-A",length(DIA))

SEXO<-rep("MACHO",length(DIA))

# Contruir la tabla =data frame, en el mismo orden que los observados

expTAM<-data.frame(TRAT,SEXO,DIA,PROP)

# Escribir tabla a archivo csv:

write.table(expTAM,"expTAM.csv",sep=',',row.names=FALSE)

4. Arreglar los valores predichos en la misma estructura que los valores observados.

Realmente, el arreglo de los datos comienza en el paso previo, construyendo las tablas con la misma estructura que el archivo de los datos observados (obdungb.csv), aunque se hizo por partes para cada tratamiento y sexo.

Una vez calculados los valores predichos de cada uno de los casos (tratamientos y sexo), los archivos de los valores predichos se combinan en uno solo empleando LIBREOFFICE u EXCEL y se guardan en un solo archivo csv, con el nombre exdunbg.csv.

5. Realizar las gráficas

Usualmente, cuando se grafican valores predichos por un modelo y los observados, los primeros (predichos) se representan con una línea y los observados con puntos. Por lo anterior, se requieren dos archivos de datos con la misma estructura (variables = columnas) y que tengan el mismo nombre las variables en el mismo orden ambos.

figCE<-ggplot(obdungb,aes(x=DIA,y=PROP,shape=TRAT)) + geom_point(size= 6,alpha= 0.7)+ theme_bw() + geom_line(data= exdungb, aes(x= DIA,y= PROP), size=1) + facet_grid(.~SEXO)

las instrucciones anteriores se interpretan de la siguiente manera:

ggplot(obdungb,aes(x=DIA,y=PROP,shape=TRAT))

Aquí, se toma la tabla obdungb para generar los puntos, con los valores de X=DIA y los de Y=PROP, instruyendo a ggplot que asigne diferentes símbolos a los tratamientos = TRAT

geom_point(size=6,alpha=0.7)

Genera la gráfica de puntos, con tamaño y opacidad definidas

theme_bw()

Es el fondo en blanco y negro de la gráfica

geom_line(data=exdungb,aes(x=DIA,y=PROP),size=1)

Esta es la instrucción para graficas los valores predichos, aquí es la parte importante pues el conjunto de datos a graficar la línea proviene de la tabla de valores predichos (=exdungb), especificando que graficar en X=DIA y en y=PROP. El tamaño de línea es 1.

facet_grid(.~SEXO)

Despliega dos paneles por tipo de sexo.

Si queremos ver la figura pero separada por tratamientos para distinguir el efecto de los sexos:

figCE2<-ggplot(obdungb,aes(x=DIA,y=PROP, shape=SEXO)) + geom_point(size=6,alpha=0.7) + theme_bw()+geom_line(data=exdungb,aes(x= DIA, y=PROP), size=1) + facet_wrap(~TRAT)

En este caso las instrucciones son similares, pero la forma de los puntos es por SEXO y los paneles es por TRAT.

II. Multiples modelos para un conjunto de datos

En esta sección se construyen gráficas con dos o mas modelos para un solo conjuntos de datos. La secuencia es similar a la anterior. Sin embargo, en este caso es necesario tener una estructura de datos como la siguiente (ver archivo obELP.csv):

Stage,Temp,dev,Modelo

larva,14,0.017140898183065,Observado

larva,20,0.03088326127239,Observado

larva,25,0.053276505061268,Observado

...

1. Selección de modelos

para este conjunto de datos, se van a estimar los parámetros de tres modelos (Briere, Davidson y Lactin; ver por ejemplo Kontodimas et al. 2004. Environ. Entomol., 33(1): 1-11 y Davidson. 1944. J. Anim. Ecol., 13(1): 26-38). Estos modelos están implementados en r en el script devmodels.R.

2. Estimación de parámetros

Se realiza para cada estado de desarrollo (larva y huevo-pupa) y se estiman los parámetros de los tres modelos, en este caso, puesto que el script devmodels. R contiene las funciones de los modelos, es posible correr el script y llamar la funcion nls() con la función respectiva. Por ejemplo, para seleccionar larvas:

larvdat<- filter(obELP, Stage=='larva')

Ahora se puede llamar la función nls(), con el ajuste del modelo de Davidson:

fit.lvdav<- nls(dev~func.davidson(Temp,aa,bb,K),data= larvdat, start= list(aa= 4, bb= 0.2,K= 10), trace= TRUE)

Lo mismo con el modelo de Briere:

fit.lvBr<- nls(dev~func.briere1(Temp,t1,t2,c1),data= larvdat, start= list(t1=10, t2= 35,c1= 0.01), trace= TRUE)

El modelo de Lactin:

fit.lvlact<- nls(dev~func.lactin2(Temp,c1,c2,c3,c4),data= larvdat, start= list(c1= 0.003, c2= 41, c3= 2.5,c4= -1 ), trace= TRUE)

En este caso, el procedimiento no tuvo convergencia, asi que no se llevó a cabo el ajuste.

Para el caso de huevoPupa si se pudo realizar el ajuste de los tres modelos.

Los parámetros estimados fueron:

Modelo de Davidson:

Modelo de Briere:

Para conocer mas con relación al ajuste de modelos no lineales, puede consultar la referencia:

Ritz, C. & J.C. Streibig. 2008. Nonlinear regression with R. Springer. New York.