Using devRate package to fit development rate models to an empirical dataset

This vignette exemplifies the use of the devRate package to fit a development rate model to an empirical dataset (laboratory experiments at air constant temperatures) and to build a simple phenological model.

The preliminary step is to perform an experiment where the arthropod study model is reared at different constant air temperatures, and to monitor the time at which individuals change of life stage (e.g., from eggs to larva). The monitoring of life stages is commonly expressed in development rate, corresponding to the inverse of the development time needed to reach the next life stage. Development time is usually expressed in days. For example, assuming that an individual needs 10 days to develop from eggs to larva at 15 degrees Celsius, its development rate would be 1/10 = 0.1. In this vignette we illustrate the different functions using a dataset from the literature.

The dataset for T. solanivora was retrieved from Crespo-Perez et al. 2011¹, using Web Plot Digitizer². The dataset is also included in the package in the exTropicalMoth object. This object is composed of the laboratory data and results of the modeling (see below). In this vignette we present how to organize your own dataset from scratch.

Organizing the dataset

In this example various individuals of T. solanivora were reared at different temperatures and the development rates were monitored for three life stages (eggs, larvae, pupae):

the eggs were reared at 10, 13, 15, 15.5, 16, 17, 20, 25, 30, and 35 degrees
the larva were reared at 10, 13, 15, 15.5, 17, 20, 25, 30, and 35 degrees
the pupa were reared at 10, 13, 15, 15.5, 16, 17, 20, 25, 30, and 35 degrees.

The dataset resulting from the rearing experiments looks like the following table (for one life stage) and represents the only input needed by the package devRate.

Temperature	Development Rate
10.0	0.031
10.0	0.039
13.0	0.072
…	…

If the dataset is stored into a file, it can be read directly from that file (e.g., using read.table()). In this example, we copy-pasted the experimental results to create vectors and data frames that are the required input types for the package.

### Experimental temperatures and development rate for the eggs
expeTempEggs <- c(10.0, 10.0, 13.0, 15.0, 15.0, 15.5, 16.0, 16.0, 17.0, 20.0, 20.0, 
              25.0, 25.0, 30.0, 30.0, 35.0)
expeDevEggs <- c(0.031, 0.039, 0.072, 0.047, 0.059, 0.066, 0.083, 0.1, 0.1, 0.1, 0.143, 
             0.171, 0.2, 0.2, 0.18, 0.001)
dfDevEggs <- data.frame(expeTempEggs, expeDevEggs)

### Experimental temperatures and development rate for the larva
expeTempLarva <- c(10.0, 10.0, 10.0, 13.0, 15.0, 15.5, 15.5, 15.5, 17.0, 20.0, 25.0, 
                   25.0, 30.0, 35.0)
expeDevLarva <- c(0.01, 0.014, 0.019, 0.034, 0.024, 0.029, 0.034, 0.039, 0.067, 0.05, 
                  0.076, 0.056, 0.0003, 0.0002)
dfDevLarva <- data.frame(expeTempLarva, expeDevLarva)

### Experimental temperatures and development rate for the pupa
expeTempPupa <- c(10.0, 10.0, 10.0, 13.0, 15.0, 15.0, 15.5, 15.5, 16.0, 16.0, 17.0, 
                  20.0, 20.0, 25.0, 25.0, 30.0, 35.0)
expeDevPupa <- c(0.001, 0.008, 0.012, 0.044, 0.017, 0.044, 0.039, 0.037, 0.034, 0.051, 0.051, 0.08, 0.092, 0.102, 0.073, 0.005, 0.0002)
dfDevPupa <- data.frame(expeTempPupa, expeDevPupa)

### Same dataset included in the package in the form of matrices
library("devRate")
data(exTropicalMoth)
str(exTropicalMoth[[1]])

## List of 3
##  $ eggs : num [1:16, 1:2] 10 10 13 15 15 15.5 16 16 17 20 ...
##  $ larva: num [1:14, 1:2] 10 10 10 13 15 15.5 15.5 15.5 17 20 ...
##  $ pupa : num [1:17, 1:2] 10 10 10 13 15 15 15.5 15.5 16 16 ...

Finding models in the literature

Before attempting to fit any model to the empirical data, the devRate function “devRateFind” search the database for previous articles fitting models to the organism, either by Order, Family, or species. The most used models are presented at the top of the data.frame.

devRateFind(orderSP = "Lepidoptera")

##              equation occu paramNumb
## 1         campbell_74  414         2
## 2          lactin2_95   32         4
## 3           logan6_76   27         4
## 4          briere1_99   23         3
## 5           taylor_81   22         3
## 6          briere2_99   18         4
## 7         davidson_44   15         3
## 8      harcourtYee_82   13         4
## 9          logan10_76   12         5
## 10    hilbertLogan_83   10         5
## 11            wang_82    9         6
## 12 sharpeDeMichele_77    8         6
## 13         lactin1_95    8         3
## 14        analytis_77    7         5
## 15              poly4    5         5
## 16         stinner_74    4         4
## 17              poly2    4         3
## 18            lamb_92    4         4
## 19      kontodimas_04    4         3
## 20           damos_08    4         3
## 21 schoolfieldHigh_81    3         4
## 22     schoolfield_81    1         6
## 23          rootsq_82    1         2
## 24       wangengel_98    1         3
## 25        regniere_12    1         6

Here we can see that at the time of this vignette, campbell_74 model was used 108 times to characterize the relationship between development rate and temperature for the Lepidoptera Order, and the taylor_81 model 22 times.

devRateFind(familySP = "Gelechiidae")

##          equation occu paramNumb
## 1     campbell_74   33         2
## 2      lactin2_95    7         4
## 3       logan6_76    4         4
## 4      lactin1_95    4         3
## 5      briere1_99    4         3
## 6        damos_08    4         3
## 7     analytis_77    3         5
## 8       taylor_81    2         3
## 9      logan10_76    1         5
## 10 harcourtYee_82    1         4
## 11          poly4    1         5
## 12      rootsq_82    1         2

If we focus on the Family (Gelechiidae), campbell_74 has been used 12 times, lactin2_95 7 times, logan6_76 4 times, lactin1_95 4 times, briere1_99 4 times, damos_08 4 times, analytis_77 3 times, taylor_81 2 times, and so on (this may change as the database is growing).

devRateFind(species = "Tecia solanivora")

##      equation occu paramNumb
## 1 campbell_74    1         2

Unfortunately, the species Tecia solanivora is not in the package database at this time, so that we have to find a closely related species. A rapid search in the literature indicates that T. solanivora is often associated with two tuber moths: Symmetrischema tangolias and Phthorimaea operculella, both being of the Gelechiidae family. The devRateFind function can be used to browse the database for these two Gelechiidae species.

devRateFind(species = "Symmetrischema tangolias")

##    equation occu paramNumb
## 1 taylor_81    2         3
## 2 rootsq_82    1         2

devRateFind(species = "Phthorimaea operculella")

##      equation occu paramNumb
## 1 campbell_74   10         2
## 2   logan6_76    3         4
## 3 analytis_77    3         5
## 4  lactin1_95    3         3
## 5  lactin2_95    3         4
## 6  briere1_99    3         3

The taylor_81 model was used for S. tangolias and among others, the lactin1_95 model was used for P. operculella.

The “devRateInfo” function provides additional information on these models and on parameter estimations.

devRateInfo(eq = taylor_81)

## ----------------------------------------
## Gauss 
## ----------------------------------------
## Taylor, F. (1981) Ecology and evolution of physiological time in insects.
## American Naturalist, 1-23. \nLamb, RJ. (1992) Developmental rate of
## Acyrthosiphon pisum (Homoptera: Aphididae) at low temperatures: implications
## for estimating rate parameters for insects. Environmental Entomology 21(1):
## 10-19.
## 
## rT ~ Rm * exp(-1/2 * ((T - Tm)/To)^2)
## <environment: namespace:devRate>
## 
## Parameter estimates from the literature (eq$startVal): 
## 
##        ordersp      familysp        genussp         species
## 1    Hemiptera     Lygaeidae       Geocoris       articolor
## 2    Hemiptera     Lygaeidae       Geocoris         pallens
## 3    Hemiptera     Lygaeidae       Geocoris       punctipes
## 4    Hemiptera       Miridae          Lygus       desetinus
## 5    Hemiptera       Miridae          Lygus        hesperus
## 6    Hemiptera     Aphididae  Acyrthosiphon           pisum
## 7    Hemiptera     Aphididae  Acyrthosiphon           pisum
## 8    Hemiptera     Aphididae    Brevicoryne       brassicae
## 9    Hemiptera     Aphididae    Brevicoryne       brassicae
## 10   Hemiptera     Aphididae      Hyadaphis pseudobrassicae
## 11   Hemiptera     Aphididae    Macrosiphum      euphorbiae
## 12   Hemiptera     Aphididae          Myzus        persicae
## 13   Hemiptera     Aphididae          Myzus        persicae
## 14   Hemiptera  Cicadellidae     Circulifer         tenelus
## 15   Hemiptera  Cicadellidae       Empoasca           fabae
## 16  Coleoptera     Bruchidae Callosobruchus     rhodesianus
## 17  Coleoptera Chrysomelidae      Crioceris        asparagi
## 18  Coleoptera Chrysomelidae         Oulema       melanopus
## 19  Coleoptera     Cucujidae   Cryptolestes     ferrugineus
## 20  Coleoptera Curculionidae     Anthonomus         grandis
## 21  Coleoptera Curculionidae         Hypera   brunneipennis
## 22  Coleoptera Curculionidae         Hypera         postica
## 23  Coleoptera   Dermestidae      Dermestes        frischii
## 24  Coleoptera   Dermestidae      Dermestes        frischii
## 25  Coleoptera   Dermestidae      Dermestes        frischii
## 26  Coleoptera   Dermestidae      Dermestes        frischii
## 27  Coleoptera Tenebrionidae      Tribolium       castaneum
## 28  Coleoptera Tenebrionidae      Tribolium        confusum
## 29 Lepidoptera     Arctiidae     Hyphantria           cunea
## 30 Lepidoptera     Noctuidae       Agrostis         segetum
## 31 Lepidoptera     Noctuidae       Amanthes        c-nigrum
## 32 Lepidoptera     Noctuidae       Mamestra     configurata
## 33 Lepidoptera     Noctuidae    Pseudaletia     unipunctata
## 34 Lepidoptera     Noctuidae         Simyra         henrici
## 35 Lepidoptera     Noctuidae     Spodoptera      frugiperda
## 36 Lepidoptera     Noctuidae      Triphaena         pronuba
## 37 Lepidoptera     Pyralidae       Anagasta      kuehniella
## 38 Lepidoptera     Pyralidae       Ostrinia       nubilalis
## 39 Lepidoptera   Tortricidae       Epiphyas     postvittana
## 40     Diptera   Chloropidae     Hippelates        bishoppi
## 41     Diptera   Chloropidae     Hippelates        pallipes
## 42     Diptera   Chloropidae     Hippelates           pusio
## 43     Diptera     Culicidae          Aedes      flavescens
## 44     Diptera     Culicidae          Aedes          vexans
## 45     Diptera     Culicidae      Anopheles quadrimaculatus
## 46     Diptera     Culicidae Toxorhynchites     brevipalpis
## 47     Diptera Drosophilidae     Drosophila    melanogaster
## 48     Diptera Drosophilidae     Drosophila    melanogaster
## 49     Diptera      Muscidae     Haematobia       stimulans
## 50     Diptera      Muscidae      Lyperosia        irritans
## 51     Diptera      Muscidae          Musca       domestica
## 52     Diptera      Muscidae       Stomoxys       calcitans
## 53 Hymenoptera    Braconidae      Apanteles     operculella
## 54 Hymenoptera    Braconidae      Apanteles     scutellaris
## 55 Hymenoptera    Braconidae      Apanteles      subandinus
## 56 Hymenoptera    Braconidae      Aphelinus      semiflavus
## 57 Hymenoptera    Braconidae       Aphidius           rapae
## 58 Hymenoptera    Braconidae         Bracon        mellitor
## 59 Hymenoptera    Braconidae         Bracon        mellitor
## 60 Hymenoptera    Braconidae          Praon        palitans
## 61 Hymenoptera    Braconidae        Trioxys          utilus
## 62   Hemiptera     Aphididae  Acyrthosiphon           pisum
## 63   Hemiptera     Aphididae  Acyrthosiphon           pisum
## 64 Lepidoptera   Gelechiidae Symmetrischema       tangolias
## 65 Lepidoptera   Gelechiidae Symmetrischema       tangolias
## 66  Coleoptera Coccinellidae         Nephus       includens
## 67  Coleoptera Coccinellidae         Nephus      bisignatus
## 68 Lepidoptera   Tortricidae          Cydia       pomonella
## 69 Lepidoptera   Tortricidae          Cydia       pomonella
## 70 Lepidoptera   Tortricidae          Cydia       pomonella
## 71 Lepidoptera   Tortricidae          Cydia       pomonella
## 72 Lepidoptera Yponomeutidae       Plutella      xylostella
## 73 Lepidoptera Yponomeutidae       Plutella      xylostella
## 74 Lepidoptera Yponomeutidae       Plutella      xylostella
## 75 Lepidoptera Yponomeutidae       Plutella      xylostella
## 76 Lepidoptera   Tortricidae  Choristoneura    occidentalis
## 77  Coleoptera Curculionidae   Dendroctonus      ponderosae
## 78  Coleoptera Chrysomelidae   Entomoscelis       americana
##                         genSp  stage param.Rm param.Tm param.To
## 1          Geocoris articolor    all  0.05500  37.2000   8.8000
## 2            Geocoris pallens    all  0.06300  37.0000   8.6000
## 3          Geocoris punctipes    all  0.04400  33.8000   7.7000
## 4             Lygus desetinus    all  0.05600  32.1000   9.8000
## 5              Lygus hesperus    all  0.06200  36.2000  12.8000
## 6         Acyrthosiphon pisum    all  0.16500  26.2000   9.0000
## 7         Acyrthosiphon pisum    all  0.15100  27.5000  11.0000
## 8       Brevicoryne brassicae    all  0.09700  26.4000  10.5000
## 9       Brevicoryne brassicae    all  0.10000  26.6000   8.4000
## 10  Hyadaphis pseudobrassicae    all  0.14400  26.0000   9.1000
## 11     Macrosiphum euphorbiae    all  0.14000  30.6000  14.6000
## 12             Myzus persicae    all  0.12900  24.7000   9.3000
## 13             Myzus persicae    all  0.15500  26.3000  10.4000
## 14         Circulifer tenelus    all  0.05200  35.2000   9.2000
## 15             Empoasca fabae    all  0.07100  30.7000   8.9000
## 16 Callosobruchus rhodesianus    all  0.03900  31.2000   8.0000
## 17         Crioceris asparagi    all  0.07100  37.6000  13.3000
## 18           Oulema melanopus    all  0.05100  33.7000  11.5000
## 19   Cryptolestes ferrugineus    all  0.04700  36.5000   9.4000
## 20         Anthonomus grandis    all  0.07500  36.7000   9.4000
## 21       Hypera brunneipennis    all  0.05800  36.8000  12.3000
## 22             Hypera postica    all  0.07400  37.8000  11.9000
## 23         Dermestes frischii 45% RH  0.01900  31.6000   6.2000
## 24         Dermestes frischii 60% RH  0.02500  33.2000   7.1000
## 25         Dermestes frischii 75% RH  0.03200  34.0000   7.7000
## 26         Dermestes frischii 90% RH  0.03900  34.3000   8.4000
## 27        Tribolium castaneum    all  0.04900  34.6000   7.1000
## 28         Tribolium confusum    all  0.03800  32.8000   7.2000
## 29           Hyphantria cunea    all  0.03100  32.9000   9.5000
## 30           Agrostis segetum    all  0.02600  33.9000  12.1000
## 31          Amanthes c-nigrum    all  0.02600  28.3000   9.8000
## 32       Mamestra configurata    all  0.02700  32.2000  12.2000
## 33    Pseudaletia unipunctata    all  0.03700  32.5000  10.8000
## 34             Simyra henrici    all  0.02900  37.7000  13.6000
## 35      Spodoptera frugiperda    all  0.05600  37.4000  11.3000
## 36          Triphaena pronuba    all  0.01800  28.8000  10.2000
## 37        Anagasta kuehniella    all  0.02600  29.8000   9.0000
## 38         Ostrinia nubilalis    all  0.04300  32.5000   9.7000
## 39       Epiphyas postvittana    all  0.03200  27.2000   9.0000
## 40        Hippelates bishoppi    all  0.07500  35.0000   9.8000
## 41        Hippelates pallipes    all  0.08400  32.1000   7.3000
## 42           Hippelates pusio    all  0.07900  32.5000   7.5000
## 43           Aedes flavescens    all  0.05100  22.2000   6.5000
## 44               Aedes vexans    all  0.14200  26.3000   7.7000
## 45  Anopheles quadrimaculatus    all  0.11400  32.8000   9.7000
## 46 Toxorhynchites brevipalpis    all  0.06200  28.4000   6.1000
## 47    Drosophila melanogaster    all  0.13100  30.2000   8.6000
## 48    Drosophila melanogaster    all  0.12200  29.2000   8.4000
## 49       Haematobia stimulans    all  0.08100  29.5000   8.8000
## 50         Lyperosia irritans    all  0.12100  34.1000  10.1000
## 51            Musca domestica    all  0.12500  33.6000   9.7000
## 52         Stomoxys calcitans    all  0.08600  32.2000   9.6000
## 53      Apanteles operculella    all  0.07400  38.1000  11.4000
## 54      Apanteles scutellaris    all  0.10000  35.2000  10.1000
## 55       Apanteles subandinus    all  0.08900  36.1000  11.5000
## 56       Aphelinus semiflavus    all  0.10000  31.8000   9.8000
## 57             Aphidius rapae    all  0.09800  27.0000   9.2000
## 58            Bracon mellitor female  0.11300  42.5000  15.3000
## 59            Bracon mellitor   male  0.10800  36.4000  11.9000
## 60             Praon palitans    all  0.08200  27.6000   7.7000
## 61             Trioxys utilus    all  0.10500  28.7000   8.8000
## 62        Acyrthosiphon pisum    all  0.19900  26.1000   9.9000
## 63        Acyrthosiphon pisum    all  0.19000  26.3000  10.3000
## 64   Symmetrischema tangolias  larva  0.04760  28.5800  10.5000
## 65   Symmetrischema tangolias   pupa  0.09900  31.8000  11.1000
## 66           Nephus includens    all  0.04410  34.9999  11.0220
## 67          Nephus bisignatus    all  0.03380  32.4999  10.7097
## 68            Cydia pomonella    egg  0.24680  30.9122   9.0574
## 69            Cydia pomonella  larva  0.06370  29.6918   8.7660
## 70            Cydia pomonella   pupa  0.08360  30.6414   8.8934
## 71            Cydia pomonella    all  0.03160  29.9478   8.5990
## 72        Plutella xylostella   eggs  0.03900  29.5200   9.3800
## 73        Plutella xylostella  larva  0.14100  26.6800   6.4500
## 74        Plutella xylostella   pupa  0.28100  31.2800   1.7300
## 75        Plutella xylostella    all  0.07800  26.7600   6.1400
## 76 Choristoneura occidentalis   eggs  0.19000  29.7000   9.5000
## 77    Dendroctonus ponderosae   eggs  0.19400  24.5000   7.6000
## 78     Entomoscelis americana   eggs  0.05063  30.5200   8.9620
##                            ref
## 1                  Taylor 1981
## 2                  Taylor 1981
## 3                  Taylor 1981
## 4                  Taylor 1981
## 5                  Taylor 1981
## 6                  Taylor 1981
## 7                  Taylor 1981
## 8                  Taylor 1981
## 9                  Taylor 1981
## 10                 Taylor 1981
## 11                 Taylor 1981
## 12                 Taylor 1981
## 13                 Taylor 1981
## 14                 Taylor 1981
## 15                 Taylor 1981
## 16                 Taylor 1981
## 17                 Taylor 1981
## 18                 Taylor 1981
## 19                 Taylor 1981
## 20                 Taylor 1981
## 21                 Taylor 1981
## 22                 Taylor 1981
## 23                 Taylor 1981
## 24                 Taylor 1981
## 25                 Taylor 1981
## 26                 Taylor 1981
## 27                 Taylor 1981
## 28                 Taylor 1981
## 29                 Taylor 1981
## 30                 Taylor 1981
## 31                 Taylor 1981
## 32                 Taylor 1981
## 33                 Taylor 1981
## 34                 Taylor 1981
## 35                 Taylor 1981
## 36                 Taylor 1981
## 37                 Taylor 1981
## 38                 Taylor 1981
## 39                 Taylor 1981
## 40                 Taylor 1981
## 41                 Taylor 1981
## 42                 Taylor 1981
## 43                 Taylor 1981
## 44                 Taylor 1981
## 45                 Taylor 1981
## 46                 Taylor 1981
## 47                 Taylor 1981
## 48                 Taylor 1981
## 49                 Taylor 1981
## 50                 Taylor 1981
## 51                 Taylor 1981
## 52                 Taylor 1981
## 53                 Taylor 1981
## 54                 Taylor 1981
## 55                 Taylor 1981
## 56                 Taylor 1981
## 57                 Taylor 1981
## 58                 Taylor 1981
## 59                 Taylor 1981
## 60                 Taylor 1981
## 61                 Taylor 1981
## 62                   Lamb 1992
## 63                   Lamb 1992
## 64       Sporleder et al. 2016
## 65       Sporleder et al. 2016
## 66      Kontodimas et al. 2004
## 67      Kontodimas et al. 2004
## 68          Aghdam et al. 2009
## 69          Aghdam et al. 2009
## 70          Aghdam et al. 2009
## 71          Aghdam et al. 2009
## 72 Marchioro and Foerster 2011
## 73 Marchioro and Foerster 2011
## 74 Marchioro and Foerster 2011
## 75 Marchioro and Foerster 2011
## 76        Regniere et al. 2012
## 77        Regniere et al. 2012
## 78            Lamb et al. 1984
## 
## Comments: 
## [...] The curve must be truncated to the right of Tm because of lethal effects
## of short exposures to high temperatures. The rate at which development rate
## falls away from Tm is measured by To. Taylor 1981

devRateInfo(eq = lactin1_95)

## ----------------------------------------
## Lactin-1 
## ----------------------------------------
## Lactin, Derek J, NJ Holliday, DL Johnson, y R Craigen (1995) Improved rate
## model of temperature-dependent development by arthropods. Environmental
## Entomology 24(1): 68-75.
## 
## rT ~ exp(aa * T) - exp(aa * Tmax - (Tmax - T)/deltaT)
## <environment: namespace:devRate>
## 
## Parameter estimates from the literature (eq$startVal): 
## 
##        ordersp      familysp       genussp      species
## 1   Coleoptera Chrysomelidae  Leptinotarsa decemlineata
## 2   Coleoptera Chrysomelidae  Leptinotarsa decemlineata
## 3   Coleoptera Chrysomelidae  Leptinotarsa decemlineata
## 4   Coleoptera Chrysomelidae  Leptinotarsa decemlineata
## 5   Coleoptera Chrysomelidae  Leptinotarsa decemlineata
## 6   Coleoptera Chrysomelidae  Entomoscelis    americana
## 7   Coleoptera Chrysomelidae  Entomoscelis    americana
## 8   Orthoptera     Acrididae    Melanoplus  sanguinipes
## 9    Homoptera     Aphididae Acyrthosiphon        pisum
## 10     Diptera    Simuliidae      Simulium     arcticum
## 11     Diptera      Muscidae         Musca    domestica
## 12 Lepidoptera   Gelechiidae   Phthorimaea  operculella
## 13 Lepidoptera   Gelechiidae   Phthorimaea  operculella
## 14 Lepidoptera   Gelechiidae   Phthorimaea  operculella
## 15 Hymenoptera Ichneumonidae      Diadegma       anurum
## 16 Hymenoptera Ichneumonidae      Diadegma       anurum
## 17 Lepidoptera Yponomeutidae      Plutella   xylostella
## 18 Lepidoptera Yponomeutidae      Plutella   xylostella
## 19 Lepidoptera Yponomeutidae      Plutella   xylostella
## 20 Lepidoptera Yponomeutidae      Plutella   xylostella
## 21 Lepidoptera   Gelechiidae          Tuta     absoluta
##                        genSp        stage param.aa param.Tmax param.deltaT
## 1  Leptinotarsa decemlineata         eggs 0.155430   38.04873     6.421234
## 2  Leptinotarsa decemlineata           L1 0.154034   38.95336     6.467896
## 3  Leptinotarsa decemlineata           L2 0.154035   37.52610     6.460183
## 4  Leptinotarsa decemlineata           L3 0.169451   36.39726     5.883764
## 5  Leptinotarsa decemlineata           L4 0.166364   35.91467     5.997446
## 6     Entomoscelis americana         eggs 0.147574   39.23976     6.704902
## 7     Entomoscelis americana        larva 0.159722   36.35009     6.257617
## 8     Melanoplus sanguinipes          all 0.121649   49.73228     8.217372
## 9        Acyrthosiphon pisum          all 0.163142   31.41951     6.110100
## 10         Simulium arcticum         eggs 0.216847   23.40644     4.603325
## 11           Musca domestica eggs + larva 0.146438   42.88196     6.821329
## 12   Phthorimaea operculella         eggs 0.177000   36.58600     5.631000
## 13   Phthorimaea operculella        larva 0.169000   37.91400     5.912000
## 14   Phthorimaea operculella         pupa 0.193000   36.29100     5.180000
## 15           Diadegma anurum          all 0.165000   35.00100     6.048000
## 16           Diadegma anurum          all 0.166000   35.00700     5.984000
## 17       Plutella xylostella         eggs 0.140000   38.02000     6.950000
## 18       Plutella xylostella        larva 0.170000   35.15000     5.880000
## 19       Plutella xylostella         pupa 0.160000   36.49000     6.220000
## 20       Plutella xylostella          all 0.170000   35.10000     5.700000
## 21             Tuta absoluta          all       NA         NA           NA
##                            ref
## 1           Lactin et al. 1995
## 2           Lactin et al. 1995
## 3           Lactin et al. 1995
## 4           Lactin et al. 1995
## 5           Lactin et al. 1995
## 6           Lactin et al. 1995
## 7           Lactin et al. 1995
## 8           Lactin et al. 1995
## 9           Lactin et al. 1995
## 10          Lactin et al. 1995
## 11          Lactin et al. 1995
## 12       Golizadeh et al. 2012
## 13       Golizadeh et al. 2012
## 14       Golizadeh et al. 2012
## 15       Golizadeh et al. 2008
## 16       Golizadeh et al. 2008
## 17 Marchioro and Foerster 2011
## 18 Marchioro and Foerster 2011
## 19 Marchioro and Foerster 2011
## 20 Marchioro and Foerster 2011
## 21         Ozgokce et al. 2016
## 
## Comments: 
## None

Here we can see for example that S. tangolias study by Sporleder et al. 2016 has used the taylor_81 model, and that P. operculella study by Golizadeh et al. 2012 has used the lactin1_95 model.

To return only the information on the closely related species, a specific query can be performed on the model.

taylor_81$startVal[taylor_81$startVal["genSp"] == "Symmetrischema tangolias",]

##        ordersp    familysp        genussp   species                    genSp
## 64 Lepidoptera Gelechiidae Symmetrischema tangolias Symmetrischema tangolias
## 65 Lepidoptera Gelechiidae Symmetrischema tangolias Symmetrischema tangolias
##    stage param.Rm param.Tm param.To                   ref
## 64 larva   0.0476    28.58     10.5 Sporleder et al. 2016
## 65  pupa   0.0990    31.80     11.1 Sporleder et al. 2016

lactin1_95$startVal[lactin1_95$startVal["genSp"] == "Phthorimaea operculella",]

##        ordersp    familysp     genussp     species                   genSp
## 12 Lepidoptera Gelechiidae Phthorimaea operculella Phthorimaea operculella
## 13 Lepidoptera Gelechiidae Phthorimaea operculella Phthorimaea operculella
## 14 Lepidoptera Gelechiidae Phthorimaea operculella Phthorimaea operculella
##    stage param.aa param.Tmax param.deltaT                   ref
## 12  eggs    0.177     36.586        5.631 Golizadeh et al. 2012
## 13 larva    0.169     37.914        5.912 Golizadeh et al. 2012
## 14  pupa    0.193     36.291        5.180 Golizadeh et al. 2012

Information from the database can be plotted using the “devRatePlotInfo” function if we want to have a first look on the taylor_81 or lactin1_95 models.

devRatePlotInfo (eq = taylor_81, sortBy = "ordersp",
  ylim = c(0, 0.20), xlim = c(0, 50))

devRatePlotInfo (eq = lactin1_95, sortBy = "ordersp",
  ylim = c(0, 1.00), xlim = c(0, 50))

If there is no a priori model selection (e.g., guided by a specific research question), the taylor_81 and/or lactin1_95 models can be used as candidate models.

Fitting models to empirical datasets

Now that we have candidate models and starting parameter estimates from closely related species, we can start the fitting process with our empirical data (NLS fit). The empirical data can be fitted to any model in the database with the “devRateModel” function, which returns an object of class “nls”.

### using the vectors from section "Organizing the dataset"
############################################################

### for the taylor_81 model
mEggs01 <- devRateModel(eq = taylor_81, 
  temp = expeTempEggs, 
  devRate = expeDevEggs,
  startValues = list(Rm = 0.05, Tm = 30, To = 5))

mLarva01 <- devRateModel(eq = taylor_81, 
  temp = expeTempLarva, 
  devRate = expeDevLarva,
  startValues = list(Rm = 0.05, Tm = 25, To = 10))
      
mPupa01 <- devRateModel(eq = taylor_81, 
  temp = expeTempPupa, 
  devRate = expeDevPupa,
  startValues = list(Rm = 0.1, Tm = 30, To = 10))

### for the lactin1_95 model
mEggs01b <- devRateModel(eq = lactin1_95, 
  temp = expeTempEggs, 
  devRate = expeDevEggs,
  startValues = list(aa = 0.177, Tmax = 36.586, deltaT = 5.631))

# mLarva01b <- devRateModel(eq = lactin1_95, 
#   temp = expeTempLarva, 
#   devRate = expeDevLarva,
#   startValues = list(aa = 0.169, Tmax = 37.914, deltaT = 5.912))
### The algorithm has not found a solution after 50 iterations
### One possibility is to increase the maximum number of iterations
### using the "control" argument (see ?nls() for more details).

mLarva01b <- devRateModel(eq = lactin1_95, 
  temp = expeTempLarva, 
  devRate = expeDevLarva,
  startValues = list(aa = 0.169, Tmax = 37.914, deltaT = 5.912), 
  control = list(maxiter = 500))
      
mPupa01b <- devRateModel(eq = lactin1_95, 
  temp = expeTempPupa, 
  devRate = expeDevPupa,
  startValues = list(aa = 0.193, Tmax = 36.291, deltaT = 5.18), 
  control = list(maxiter = 500))

### using the data frames from section "Organizing the dataset"
############################################################

mEggs02 <- devRateModel(eq = taylor_81, 
  dfData = dfDevEggs,
  startValues = list(Rm = 0.05, Tm = 30, To = 5))

mLarva02 <- devRateModel(eq = taylor_81, 
  dfData = dfDevLarva,
  startValues = list(Rm = 0.05, Tm = 25, To = 10))
      
mPupa02 <- devRateModel(eq = taylor_81, 
  dfData = dfDevPupa,
  startValues = list(Rm = 0.1, Tm = 30, To = 10))

### using the dataset included in the package (only for taylor_81 model)
############################################################

mEggs <- devRateModel(eq = taylor_81, 
  temp = exTropicalMoth$raw$eggs[,1], 
  devRate = exTropicalMoth$raw$eggs[,2],
  startValues = list(Rm = 0.05, Tm = 30, To = 5))
      
mLarva <- devRateModel(eq = taylor_81, 
  temp = exTropicalMoth$raw$larva[,1], 
  devRate = exTropicalMoth$raw$larva[,2],
  startValues = list(Rm = 0.05, Tm = 25, To = 10))
      
mPupa <- devRateModel(eq = taylor_81, 
  temp = exTropicalMoth$raw$pupa[,1], 
  devRate = exTropicalMoth$raw$pupa[,2],
  startValues = list(Rm = 0.1, Tm = 30, To = 10))

The summary of the “devRateModel” can be obtained with the “devRatePrint” function.

resultNLS <- devRatePrint(myNLS = mLarva)

## ##################################################
## ### Parameter estimates and overall model fit
## ##################################################
## 
## Formula: rT ~ Rm * exp(-1/2 * ((T - Tm)/To)^2)
## 
## Parameters:
##     Estimate Std. Error t value Pr(>|t|)    
## Rm  0.068383   0.009473   7.219 1.71e-05 ***
## Tm 21.408732   0.729062  29.365 8.41e-12 ***
## To  5.544986   0.865264   6.408 5.02e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.0133 on 11 degrees of freedom
## 
## Number of iterations to convergence: 13 
## Achieved convergence tolerance: 6.657e-06
## 
## ##################################################
## ### Confidence intervals for parameters
## ##################################################
##          2.5 %      97.5 %
## Rm  0.04981547  0.08694955
## Tm 19.97979623 22.83766819
## To  3.84909961  7.24087243
## 
## ##################################################
## ### Residuals distribution and independence
## ##################################################
## ### Normality of the residual distribution
## 
##  Shapiro-Wilk normality test
## 
## data:  stats::residuals(myNLS)
## W = 0.97915, p-value = 0.9696
## 
## ### Regression of the residuals against a lagged version of themselves
## ### and testing if the slope of the resulting relationship is significantly
## ### different from 0:
## 
## Call:
## stats::lm(formula = stats::residuals(myNLS)[-N] ~ stats::residuals(myNLS)[-1])
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.021388 -0.011116  0.002149  0.009943  0.020078 
## 
## Coefficients:
##                               Estimate Std. Error t value Pr(>|t|)
## (Intercept)                  0.0005691  0.0036255   0.157    0.878
## stats::residuals(myNLS)[-1] -0.1654787  0.2965227  -0.558    0.588
## 
## Residual standard error: 0.01307 on 11 degrees of freedom
## Multiple R-squared:  0.02753,    Adjusted R-squared:  -0.06087 
## F-statistic: 0.3114 on 1 and 11 DF,  p-value: 0.588
## 
## ##################################################
## ### Comparing models
## ##################################################
## ### Using AIC and BIC
## Akaike Information Criterion (AIC): -76.600135294506
## Bayesian Information Criterion (BIC): -74.043905976045

resultNLSb <- devRatePrint(myNLS = mLarva01b)

## ##################################################
## ### Parameter estimates and overall model fit
## ##################################################
## 
## Formula: rT ~ exp(aa * T) - exp(aa * Tmax - (Tmax - T)/deltaT)
## 
## Parameters:
##        Estimate Std. Error t value Pr(>|t|)    
## aa      0.10600    0.02346   4.519 0.000874 ***
## Tmax   34.36276    1.09320  31.433 4.01e-12 ***
## deltaT  9.40120    2.05900   4.566 0.000809 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.01864 on 11 degrees of freedom
## 
## Number of iterations to convergence: 101 
## Achieved convergence tolerance: 9.417e-06
## 
## ##################################################
## ### Confidence intervals for parameters
## ##################################################
##              2.5 %    97.5 %
## aa      0.06002062  0.151977
## Tmax   32.22011918 36.505397
## deltaT  5.36564562 13.436764
## 
## ##################################################
## ### Residuals distribution and independence
## ##################################################
## ### Normality of the residual distribution
## 
##  Shapiro-Wilk normality test
## 
## data:  stats::residuals(myNLS)
## W = 0.95297, p-value = 0.6078
## 
## ### Regression of the residuals against a lagged version of themselves
## ### and testing if the slope of the resulting relationship is significantly
## ### different from 0:
## 
## Call:
## stats::lm(formula = stats::residuals(myNLS)[-N] ~ stats::residuals(myNLS)[-1])
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.037702 -0.008981 -0.000671  0.005363  0.029264 
## 
## Coefficients:
##                              Estimate Std. Error t value Pr(>|t|)
## (Intercept)                 -0.001559   0.005071  -0.307    0.764
## stats::residuals(myNLS)[-1]  0.073343   0.306144   0.240    0.815
## 
## Residual standard error: 0.01828 on 11 degrees of freedom
## Multiple R-squared:  0.005191,   Adjusted R-squared:  -0.08525 
## F-statistic: 0.05739 on 1 and 11 DF,  p-value: 0.8151
## 
## ##################################################
## ### Comparing models
## ##################################################
## ### Using AIC and BIC
## Akaike Information Criterion (AIC): -67.1607406792436
## Bayesian Information Criterion (BIC): -64.6045113607826

Empirical data can be plotted against the model using the “devRatePlot” function.

par(mfrow = c(1, 2), mar = c(4, 4, 0, 0))
devRatePlot(eq = taylor_81, 
  nlsDR = mEggs, 
  pch = 16, ylim = c(0, 0.2))

devRatePlot(eq = lactin1_95, 
  nlsDR = mEggs01b, 
  pch = 16, ylim = c(0, 0.2))

devRatePlot(eq = taylor_81, 
  nlsDR = mLarva, 
  pch = 16, ylim = c(0, 0.1))

devRatePlot(eq = lactin1_95, 
  nlsDR = mLarva01b, 
  pch = 16, ylim = c(0, 0.1))

devRatePlot(eq = taylor_81, 
  nlsDR = mPupa, 
  pch = 16, ylim = c(0, 0.15))

devRatePlot(eq = lactin1_95, 
  nlsDR = mPupa01b, 
  pch = 16, ylim = c(0, 0.15))

Models comparison

Models can be compared using the AIC or any function compatible with an “nls” object, such as BIC or logLik (see the R manual for the use and interpretation of these functions, outside of the scope of this vignette).

### Models for the larva life stage
c(AIC(mLarva), AIC(mLarva01b))

## [1] -76.60014 -67.16074

c(BIC(mLarva), BIC(mLarva01b))

## [1] -74.04391 -64.60451

c(logLik(mLarva), logLik(mLarva01b))

## [1] 42.30007 37.58037

Forecasting phenologies from empirical temperature

From the “nls” object obtained above, we can build a simple phenology model using temperature time series (e.g., to forecast the organism outbreaks). In this example the temperature dataset is built from a Normal distribution of mean 15 and a standard deviation of 1, with a time step of one day. The development models used are those previously fitted with the Taylor model for the three stages. We assumed that the average time for female adults to lay eggs was of 1 day. We simulated 50 individuals, with a stochasticity in development rate centered on the development rate, with a standard deviation of 0.015 (Normal distribution).

forecastTsolanivora <- devRateIBM(
  tempTS = rnorm(n = 100, mean = 15, sd = 1),
  timeStepTS = 1,
  models = list(mEggs, mLarva, mPupa),
  numInd = 50,
  stocha = 0.015,
  timeLayEggs = 1)
print(forecastTsolanivora)

## [[1]]
##       g1s1 g1s2 g1s3 g2s1
##  [1,]   18   49   74   90
##  [2,]   17   48   76   91
##  [3,]   17   45   74   90
##  [4,]   18   46   73   91
##  [5,]   17   47   72   89
##  [6,]   18   44   73   91
##  [7,]   17   47   74   91
##  [8,]   20   49   78   94
##  [9,]   18   50   81   97
## [10,]   17   43   69   85
## [11,]   18   48   74   90
## [12,]   17   46   72   89
## [13,]   18   42   71   88
## [14,]   18   48   81   97
## [15,]   19   48   76   94
## [16,]   17   47   74   90
## [17,]   18   44   71   88
## [18,]   18   47   75   91
## [19,]   18   49   76   94
## [20,]   17   50   81   97
## [21,]   16   47   75   91
## [22,]   17   44   76   92
## [23,]   17   44   73   89
## [24,]   19   45   72   89
## [25,]   18   45   72   89
## [26,]   17   52   80   96
## [27,]   18   48   74   90
## [28,]   18   47   73   89
## [29,]   18   47   74   89
## [30,]   18   47   75   92
## [31,]   17   49   82   97
## [32,]   18   47   75   92
## [33,]   16   45   72   88
## [34,]   18   48   79   95
## [35,]   19   46   78   94
## [36,]   18   46   71   88
## [37,]   18   46   74   91
## [38,]   17   47   73   89
## [39,]   20   51   81   96
## [40,]   18   50   79   96
## [41,]   21   48   77   93
## [42,]   18   44   73   89
## [43,]   17   44   76   91
## [44,]   18   42   73   90
## [45,]   15   45   72   88
## [46,]   20   49   76   91
## [47,]   18   46   76   92
## [48,]   18   49   77   91
## [49,]   16   46   77   94
## [50,]   18   47   78   95
## 
## [[2]]
## [[2]][[1]]
## Nonlinear regression model
##   model: rT ~ Rm * exp(-1/2 * ((T - Tm)/To)^2)
##    data: data.frame(rT = devRate, T = temp)
##      Rm      Tm      To 
##  0.1934 25.3418 -6.8939 
##  residual sum-of-squares: 0.01303
## 
## Number of iterations to convergence: 13 
## Achieved convergence tolerance: 8.255e-06
## 
## [[2]][[2]]
## Nonlinear regression model
##   model: rT ~ Rm * exp(-1/2 * ((T - Tm)/To)^2)
##    data: data.frame(rT = devRate, T = temp)
##       Rm       Tm       To 
##  0.06838 21.40873  5.54499 
##  residual sum-of-squares: 0.001947
## 
## Number of iterations to convergence: 13 
## Achieved convergence tolerance: 6.657e-06
## 
## [[2]][[3]]
## Nonlinear regression model
##   model: rT ~ Rm * exp(-1/2 * ((T - Tm)/To)^2)
##    data: data.frame(rT = devRate, T = temp)
##      Rm      Tm      To 
##  0.0978 22.0114  4.7515 
##  residual sum-of-squares: 0.002394
## 
## Number of iterations to convergence: 14 
## Achieved convergence tolerance: 3.906e-06
## 
## 
## [[3]]
##   [1] 13.70278 13.74307 14.09737 14.69949 14.43374 14.15131 14.52781 16.32328
##   [9] 13.86916 13.99791 14.08161 14.81894 15.56448 15.66115 14.48841 15.22075
##  [17] 15.26903 15.77744 13.74930 14.19039 15.75059 15.22916 14.06040 13.89372
##  [25] 14.22037 13.24711 14.32201 14.98959 13.36837 15.44089 16.40496 15.65374
##  [33] 14.66768 15.29338 14.19792 15.23356 14.82910 13.68258 14.90226 13.82633
##  [41] 15.58481 17.12122 14.71389 15.00275 14.53283 14.83400 16.81026 14.27259
##  [49] 14.22460 16.11355 17.34444 15.71858 15.83886 14.85734 14.86458 14.64696
##  [57] 14.86827 14.40930 14.20924 15.39221 14.95297 14.03661 14.85702 14.67295
##  [65] 17.03639 14.77553 15.67051 16.30841 14.43236 14.95410 14.61064 16.03883
##  [73] 14.93001 15.62551 14.01097 15.29350 13.92375 15.95290 15.06230 15.14135
##  [81] 15.69510 14.26808 15.06204 15.38073 14.56196 13.45738 16.40630 16.09072
##  [89] 15.74254 14.74044 15.45186 14.48353 15.18935 15.62541 15.95085 15.15426
##  [97] 15.76322 14.46496 14.28494 15.66517

Results can be plotted using the “devRateIBMPlot” function. From this simple model we can observe that if eggs of the first generation are laid at time t = 0, adults should be seen at t = [65:85], if the temperature time series correspond to the temperatures experienced by the organism and in the absence of other limiting factors.

devRateIBMPlot(ibm = forecastTsolanivora, typeG = "density")

## Threshold for visualization = 10% of individuals

Crespo-Pérez, V., Rebaudo, F., Silvain, J.-F. and Dangles, O. (2011) Modeling invasive species spread in complex landscapes: the case of potato moth in Ecuador. Landscape Ecology, 26, 1447–1461.↩︎
Rohatgi, A. (2015) WebPlotDigitalizer: HTML5 Based Online Tool to Extract Numerical Data from Plot Images.↩︎