search

kingaa / pomp

R package for statistical inference using partially observed Markov processes
https://kingaa.github.io/pomp
GNU General Public License v3.0
111 stars 27 forks source link

Inconsistency between R-Desolve Model and Pomp Csnippet! #78

Closed pakdamie closed 5 years ago

pakdamie commented 5 years ago

I'm utilizing distributed delay differential equations to model some insect population dynamics. The image below represents my overall model here where I have individuals moving in and out of life-stages. I introduced the distributed delay by using the linear-chain trick. The external variables that would be driving this model would mainly be the daily temperatures which I have as a vector.

image

It's a simple compartmental model except when it comes to Larval Instar 4, Larval Instar 5, Diapausing Larval Instar 4, and Diapausing Larval Instar 5. Individuals in larval instar 4 can either get out of the subcompartment by continuing development, dying through background mortality, or by going into the first subcompartment of the stage Diapausing Instar 4. The rate of going into this life-stage changes through the season but remains the same for each subcompartment. Individuals in the Diapausing Instar 4 will have to spend time in this stage and after some period, the individuals will go back into the first sub-compartment of the Larval Instar 4. This is the same thing that will be seen in the fifth instar/diapausing fifth instar larvae.

Because I'm using a series of ODEs and running it for a large time-step, I was trying to use pomp to make it faster and to later use it for some parameter estimation. I combined some of the documentations and guide to finagle out a pomp code. The problem I'm having right now is that my original R desolve model and my Pomp C-snippet, even if they have the exact same parameters and external variables, seem to be giving me very, very different results.

Here is my original deSolve code here: (I ran the output for 5 years so 1825 days)

library(deSolve)
library(geosphere)
library(here)

###You need to load the temperature function
###(Alternate way)###
###load(here("Data","sub.RData"))

load(url("https://github.com/pakdamie/modeling_diapause_CM/blob/master/Data/sub.RData?raw=true"))

PHOTO_PERIOD_GEOSPHERE <- geosphere::daylength(40, 1:366)
DIFF <- diff(PHOTO_PERIOD_GEOSPHERE)
###Change this based on how many years you watnt to run out for
PP_DIFF= c(rep(DIFF,length.out=365*5))

MODEL_DDE_CM_DIA = function(t, y, parms) {

  ###This is for controlling the erlang distribution in each stage.

  nE=8#The number of subcompartments in the egg-stage
  nL1 =8 #The number of subcompartments in larval instar 1
  nL2 =8 #The number of subcompartments in larval instar 2
  nL3 = 8#The number of subcompartments in larval instar 3
  nL4 =8#The number of subcomppartmetns in larval instar 4
  nL5=8#The number of subcompartments in larval instar 5
  nDL4 =8  #The number of subcompartments in diapausing larval instar 4
  nDL5 =8 #The number of subcomparments in diapausing larval instar 5
  nP =8 # The number of subcompartments in pupae
  nA.r =8 #The number of subcompartments in reproductive adults

  ###The environmental variable- this is the daily mean temperature 
  ###from Biglerville, PA. It contains 33 years of data 
  ### (January 1st, 1984- December 31, 2016)
  TT =(SUBSET_TEMP)$avg[t]
  DIFF = (PP_DIFF)[t]

  ###The parameters that I am interested in (###IGNORE HERE: parameter values
### were directly put into the code and not as a separate list ###

  ###Development
  ###The parameter associated with the "diapause termination rate"
  ###for individuals in the diapausing instar 4 and diapausing instar 5
  alphaDL_a = parms[1]; alphaDL_b = parms[2];  alphaDL_c = parms[3];

  ###Mortality
  ###The parameters associated with the mortality of the Adults
  deltaAa= parms[4];deltaAb= parms[5];deltaAc= parms[6];
  ###The parameters associated with the mortality of the Diapausing Larvae
  deltaDLa= parms[7]; deltaDLb= parms[8]; deltaDLc = parms[9];

  ###The parameter associated with diapause induction:
  ###Diapause induction rate 1 is for larvae 4
  ###Diapause induction rate 2 is for larvae 5
   DIA_INDUC1_A= parms[10];DIA_INDUC1_B = parms[11];
   DIA_INDUC2_A=parms[12];DIA_INDUC2_B = parms[13]

  #########
  #STAGES##
  #########
  ###The stages
  E = y[(1:nE)]
  L1 = y[(nE+1):(nE + nL1)]
  L2=  y[(nE+nL1+1):(nE+nL1 + nL2 )]
  L3 = y[(nE+nL1+nL2+1):(nE+nL1+nL2+nL3)]
  L4 =  y[(nE+nL1+nL2+nL3+ 1):(nE+nL1+nL2+nL3+nL4)]
  L5 =  y[(nE+nL1+nL2+nL3+nL4 +1):(nE+nL1+nL2+nL3+nL4+nL5)]
  DL4 = y[(nE+nL1+nL2+nL3+nL4+nL5+1):(nE+nL1+nL2+nL3+nL4+nL5+nDL4)]
  DL5 = y[(nE+nL1+nL2+nL3+nL4+nL5+nDL4+ 1):(nE+nL1+nL2+nL3+nL4+nL5+nDL4+nDL5)]
  P =  y[(nE+nL1+nL2+nL3+nL4+nL5+nDL4+nDL5+1):(nE+nL1+nL2+nL3+nL4+nL5+nDL4+nDL5+nP)]
  A.r= y[(nE+nL1+nL2+nL3+nL4+nL5+nDL4+nDL5+nP+1):(nE+nL1+nL2+nL3+nL4+nL5+nDL4+nDL5+nP+nA.r) ]
  A.s = y[(nE+nL1+nL2+nL3+nL4+nL5+nDL4+nDL5+nP+nA.r+1) ]

  #Birth Rate
  B =10.52* exp(-1 / 2 * (TT - 24.629) ^ 2 / ( 3.427  ^ 2))

  #Development rate (including the diapause termination rate)
  alphaE <- 0.1927/(1+ exp(-0.3039*(TT - 18.5929)))
  alphaL1 <- 0.30/(1+ exp(- 0.327*(TT - 17.60  )))
  alphaL2 <- 0.457/(1+ exp(- 0.2301*(TT - 21.23)))
  alphaL3 <- 0.338/(1+ exp(- 0.350*(TT - 18.45)))
  alphaL4 <- 0.305/(1+ exp(- 0.429*(TT - 21.54)))
  alphaL5 <- 0.33/(1+ exp(- 0.2144*(TT - 20.94)))
  alphaP <- 0.09287/(1+ exp(-0.28966 * (TT - 18.48736 )))
  alphaA <-   0.1707 /(1+ exp(-0.1830*(TT - 20.2206)))
  alphaDL <- 0.7543/(1+ exp(0.2392 *(TT + 7.1729 )))

  ###THE MORTALITY FUNCTIONS
  deltaE <- 0.00056*(TT)^2 -0.038 *TT + 0.65
  deltaL <- 0.00067*(TT)^2 -0.017*TT+0.04
  deltaP <- -1.598e-05*(TT)^2 + (6.787e-05  *TT) + 1.606e-02
  deltaA <- 3.328e-04*(TT)^2-  1.179e-03  *TT + 2.452e-02 
  deltaDL <-  3.328e-05*(TT)^2-  1.179e-03  *TT + 2.452e-02 

  ###DIAPAUSE INDUCTION RATES- I used the change in photoperiod for right now 

  DIA_INDUC1 <-0.2130/(1+exp(80 *DIFF))
  DIA_INDUC2 <-0.5130/(1+exp(80 *DIFF))

  ######
  #Egg##
  ######
  dE= rep(0,nE)
  ###For the first substage 
  dE[1] = B * sum(A.r) - (nE*alphaE + deltaE)*E[1] 
  for (i in seq(2, nE)){
    dE[i] = nE*alphaE*(E[i-1]- E[i]) - deltaE*E[i]
  } 
  ###LARVAL INSTAR 1
  dL1= rep(0,nL1)
  dL1[1] = nE*alphaE*E[nE] - (nL1 *alphaL1 +deltaL) * L1[1]

  for (i in seq(2, nL1)){
    dL1[i] = nL1 * alphaL1  * (L1[i - 1] - L1[i]) - deltaL * L1[i]   
  } 

  ###LARVAL INSTAR 2
  dL2= rep(0,nL2)
  dL2[1] = nL1*alphaL1*L1[nL1] - (nL2 * alphaL2 + deltaL) * L2[1] 

  for (i in seq(2, nL2)){
    dL2[i] = nL2 *   alphaL2  * (L2[i-1]-L2[i]) - deltaL*L2[i] 
  } 
  ###LARVAE 3
  dL3= rep(0,nL3)
  dL3[1] = (nL2*alphaL2*L2[nL2]) - (nL3 * alphaL3 * L3[1]) - deltaL * L3[1]  

  for (i in seq(2, nL3)){
    dL3[i] = nL3 *   alphaL3  * (L3[i - 1] - L3[i]) -  
      deltaL * L3[i] #Background mortality 
  } 

  ###LARVAE 4
  dL4= rep(0,nL4)

  dL4[1] = nL3 * alphaL3 * L3[nL3]  + 
    nDL4 * alphaDL* DL4[nDL4] - 
    nL4 * alphaL4 * L4[1]- 
    DIA_INDUC1*L4[1]-                       
    (deltaDL* L4[1])          

  for (i in seq(2, nL4)){
    dL4[i] =  nL4 *   alphaL4  * (L4[i - 1]- L4[i])-   
      (DIA_INDUC1*L4[i])- #Recruitment into the diapausing stage 
      (deltaL * L4[i]) #Background mortality 
  } 
  ###LARVAE 5

  dL5= rep(0,nL5)

  dL5[1] = nL4*alphaL4*L4[nL4] + 
    nDL5 * alphaDL* DL5[nDL5] -          
    nL5 * alphaL5  * L5[1] -
    DIA_INDUC2*L5[1] -
    (deltaL*L5[1])                    

  for (i in seq(2, nL5)){
    dL5[i] = nL5 *   alphaL5  * (L5[i - 1]-L5[i])-  #Recruitment from the previous stage
      #Recruitment intot he next substage      
      (  DIA_INDUC2 * L5[i])-                       #Recruitment into the diapausing stage
      (deltaL * L5[i])                             #Background mortality 
  } 

  ###The diapausing stage - 4th instar
  dDL4 = rep(0,nDL4)
  dDL4[1] =DIA_INDUC1 * sum(L4)-  #Recruitment from the fourth instar larvae
    alphaDL *nDL4*DL4[1] - #Recruitment out into the next substage
    deltaDL*DL4[1] #Backgroudn mortality 

  for (i in seq(2, nDL4))
  {
    dDL4[i] = nDL4 *   alphaDL  * (DL4[i - 1]-DL4[i])-  
      deltaDL * DL4[i] #Background mortality 
  } 
  ###The diapausing stage - 5th instar
  dDL5 = rep(0,nDL5)
  dDL5[1] =   DIA_INDUC2  * sum(L5) - #Diapause in
    alphaDL * nDL5 * DL5[1]- #Development out 
    deltaDL*DL5[1] # Mortality

  for (i in seq(2, nDL5)){
    dDL5[i] =nDL5 *   alphaDL  * (DL5[i - 1]-DL5[i]) - 
      deltaDL * DL5[i]  }

  #######
  #Pupae#
  #######
  dP = rep(0,nP)
  dP[1] =  nL5*alphaL5* L5[nL5] - #Recruitment in from previous stage
    (nP * alphaP * P[1]) - deltaP * P[1] #
  for (i in seq(2, nP))
  {
    dP[i] = nP * alphaP * (P[i - 1] - P[i])-  #Recruitment in from prevoius substage
      deltaP* P[i] #Backgroudn mortality 
  }

  #Reproductive Adult 
  dA.r = rep(0,nA.r)
  dA.r[1] = nP * alphaP* P[nP] - #Recruitment in from previous stage
    (nA.r* alphaA + deltaA ) * A.r[1]  #Recruitment out
  for (i in seq(2, nA.r))
  {
    dA.r[i] = nA.r* alphaA * (A.r[i - 1] - A.r[i]) -deltaA* A.r[i]
  }

  #senscent adult
  dA.s = nA.r*alphaA * A.r[nA.r]  - deltaA*A.s

  res = c(dE,dL1,dL2,dL3,dL4,dL5,dDL4,dDL5,dP,dA.r,dA.s)
  list(res)

}
tstep = seq(1,365*5,1)

nE=8 #The number of subcompartments in the egg-stage
nL1 =8 #The number of subcompartments in larval instar 1
nL2 =8 #The number of subcompartments in larval instar 2
nL3 = 8 #The number of subcompartments in larval instar 3
nL4 =8#The number of subcomppartmetns in larval instar 4
nL5=8 #The number of subcompartments in larval instar 5
nDL4 =8  #The number of subcompartments in diapausing larval instar 4
nDL5 =8 #The number of subcomparments in diapausing larval instar 5
nP =8 # The number of subcompartments in pupae
nA.r =8

yinit = c(
  E = rep(0, nE),
  L1 = rep(0,nL1),
  L2 = rep(0,nL2),
  L3 = rep(0,nL3),
  L4 = rep(0,nL4),
  L5 = rep(0,nL5),
  DL4 =rep(2,nDL4),
  DL5 =rep(0,nDL5),
  P = rep(0,nP),
  A.r = rep(0,nA.r),
  A.s = 0
)

###These doesn't really matter, cause i'm inserting the value in directly to 
### the code for debugging purposes
p <- c( 
  alphaDL_a =0.7543 , alphaDL_b =0.2392,  alphaDL_c =7.1729,
  deltaAa=  3.328e-05,deltaAb= 1.179e-03,deltaAc=2.452e-02 ,
  deltaDLa=  3.328e-05, deltaDLb= 1.179e-03 , deltaDLc = 2.452e-02 ,
  DIA_INDUC1_A=0.2130,DIA_INDUC1_B = 80.0004,
  DIA_INDUC2_A= 0.5130,DIA_INDUC2_B = 80.0004)

yout<- ode(
  y = yinit ,
  times = tstep,
  func = MODEL_DDE_CM_DIA,
  parms=p)

E<-(rowSums(yout[,(2:(nE+1))]))
plot(E,type='l',main='E')

L1<- rowSums(yout[,(nE+2):(nE+nL1+1)])
plot(L1,type='l',main='L1')

L2 <- rowSums(yout[,(nE+nL1+2):(nE+nL1+nL2 +1)])
plot(L2, type='l',main ='L2')

L3 <- rowSums(yout[,(nE+nL1+nL2+2):(nE+nL1+nL2 + nL3+1)])
plot(L3, type='l',main ='L3')

L4 <- (rowSums(yout[,(nE+nL1+nL2+nL3+2):(nE+nL1+nL2 + nL3+nL4+1)]))
plot(L4, type='l',main ='L4')

L5 <- rowSums(yout[,(nE+nL1+nL2+nL3+nL4+2):(nE+nL1+nL2 + nL3+nL4+nL5+1)])
plot(L5, type='l',main ='L5')

DL4 <-as.double(rowSums(yout[,(nE+nL1+nL2+nL3+nL4+nL5+2):(nE+nL1+nL2 + nL3+nL4+nL5+nDL4 + 1)]))
plot(DL4,main ='DL4',type='l')

DL5<- rowSums(yout[,(nE+nL1+nL2 + nL3+nL4+nL5+nDL4 + 2):(nE+nL1+nL2 + nL3+nL4+nL5+nDL4 +nDL5+ 1)])
plot(DL5, type='l',main ='DL5')

P<- rowSums(yout[,(nE+nL1+nL2 + nL3+nL4+nL5+nDL4 + nDL5+2):(nE+nL1+nL2 + nL3+nL4+nL5+nDL4 +nDL5+nP+ 1)])
plot(P, type='l',main ='P')

A.r<- rowSums(yout[,(nE+nL1+nL2 + nL3+nL4+nL5+nDL4 + nDL5+nP+2):
                     (nE+nL1+nL2 + nL3+nL4+nL5+nDL4 +nDL5+nP+nA.r +1)])
plot(scale(A.r),type='l',main='A.R')

Here is the output of the reproductive adults: image

Here is my pomp code here where I only use the deterministic skeleton to use the trajectory function.

library(here)
library(pomp)

###This loads the covariate table which has the temperature and photoperiod data 
###alternate###
###load(here("Data","covar.RData"))
load(url("https://github.com/pakdamie/modeling_diapause_CM/blob/master/Data/covar.RData?raw=true"))

###This sets the number of subcompartments, or the number within each array 
n = 8                          
nE <- as.integer(n)
nL1 <-  as.integer(n)
nL2 <-  as.integer(n)
nL3 <-  as.integer(n)
nL4 <-  as.integer(n)
nL5 <-  as.integer(n)
nP <-  as.integer(n)
nAr <-  as.integer(n)
nDL4 <-  as.integer(n)
nDL5  <-  as.integer(n)

globs <- c(paste0("static int nE = ",nE,";"),
           paste0("static int nL1 = ",nL1,";"),
           paste0("static int nL2 = ",nL2,";"),
           paste0("static int nL3 = ",nL3,";"),
           paste0("static int nL4 = ",nL4,";"),
           paste0("static int nL5 = ",nL5,";"),
           paste0("static int nDL4 = ",nDL4,";"),
           paste0("static int nDL5 = ",nDL5,";"),
           paste0("static int nP = ",nP,";"),
           paste0("static int nAr = ",nAr,";"))

DIA_SKEL <- Csnippet("
                      //I need this function for the birth rate and diapause rate
                      // example- the number of eggs is dependent on all the reproductive adults 
                     // in the subcompartments. So total eggs = Birth Rate * All Adults

                      const double  sum(const double arr[], int n) 
                      { 
                   double  sum = 0.0;  
                         for (int i = 0; i < n; i++) 
                        sum += arr[i]; 

                           return sum; 
                      } 

                  // Here are the states that will require pointers to form 
                 // an array of vectors 

                    double  *E = &E1;
                    double *DE = &DE1;

                    double *L1 = &L11;
                    double *DL1 = &DL11;

                    double *L2 = &L21;
                    double *DL2 = &DL21;

                    double *L3 = &L31;
                    double *DL3 = &DL31;

                    double  *L4= &L41;
                    double *DL4 = &DL41;

                    double *L5= &L51;
                    double *DL5 = &DL51;

                    double *DIAL4= &DIAL41;
                    double *DDIAL4 = &DDIAL41;

                    double *DIAL5= &DIAL51;
                    double *DDIAL5 = &DDIAL51;

                    double *P= &P1;
                    double *DP = &DP1;

                    double *Ar= &Ar1;
                    double *DAr = &DAr1;

                  int i; 

                  //Birth Rate
                  long  double B;
                   B =  10.52 * exp(-0.50* pow ((TT- 24.629), 2.0) / pow (3.427 , 2.0));

                   //Development Rate
                   double alphaE, alphaL1, alphaL2, alphaL3, 
                   alphaL4, alphaL5, alphaP, alphaA, alphaDL;

                     alphaE  =  0.1927/ (1+ exp(-0.3039*(TT - 18.5929)));
                     alphaL1  =  0.30/(1+ exp(-0.327*(TT - 17.60  )));
                     alphaL2 =  0.457/(1+ exp(- 0.2301 * (TT - 21.23)));
                     alphaL3  =  0.338/(1+ exp(- 0.350*(TT - 18.45)));
                     alphaL4 = 0.305/(1+ exp(- 0.429*(TT - 21.54)));
                     alphaL5 =  0.33/(1+ exp(- 0.2144*(TT - 20.94)));
                     alphaP  = 0.09287/(1+ exp(-0.28966*(TT - 18.48736 )));
                     alphaA =    0.1707/(1+ exp(-0.1830*(TT - 20.2206 )));
                     alphaDL =  0.7543 / (1+ exp(0.2392 *(TT + 7.1729)));

                      //Mortality Rate
                    double muE, muL, muP, muA, muDL;

                     muE = 0.00056 * pow(TT, 2)  - (0.038* TT) + 0.65;

                     muL =  0.00067 * pow(TT, 2)  - (0.017 * TT) + 0.04;

                     muP =  -1.598e-05 * pow(TT, 2)  - (6.787e-05  * TT) + 1.606e-02;

                     muA =3.328e-04 * pow(TT, 2)  - ( 1.179e-03 * TT) +  2.452e-02 ;

                     muDL =3.328e-05 * pow(TT, 2)  - ( 1.179e-03 * TT) + 2.452e-02;

                     //Diapause induction

                      double DIA_INDUC_L4;
                      double DIA_INDUC_L5;

                     DIA_INDUC_L4=0.2130/( 1+exp(80 *PP_DIFF));
                     DIA_INDUC_L5 = 0.5130/( 1+exp(80 *PP_DIFF));

                     //THE DIFFERENTIAL EQUATIONS                      
                     //EGG
                     DE[0] = B*sum(Ar,nAr)  - (nE * alphaE* E[0]) - (muE* E[0]);
                     for ( i = 1; i  < nE;  i ++){
                     DE[i] = nE * alphaE * (E[i-1]- E[i]) - muE * E[i];}

                     // LARVAE 1
                     DL1[0] = nE * alphaE * E[nE-1] - (nL1 * alphaL1  * L1[0])-  (muL*L1[0]);
                     for ( i = 1;  i < nL1; i++){
                     DL1[i] = nL1 * alphaL1 * (L1[i-1]- L1[i])- muL * L1[i];}

                     //LARVAE 2
                     DL2[0] =  nL1 * alphaL1 * L1[nL1-1]  - (nL2 * alphaL2 * L2[0])- (muL*L2[0]);
                     for (i = 1;  i < nL2 ;  i++){
                     DL2[i] = nL2 * alphaL2 * (L2[i-1]- L2[i])- muL*L2[i];}

                     //LARVAE 3
                     DL3[0] =  nL2 * alphaL2 * L2[nL2-1] -  (nL3 * alphaL3 * L3[0])- (muL*L3[0]);
                     for (i = 1; i < nL3; i++){
                     DL3[i]= nL3 * alphaL3 *(L3[i-1]- L3[i])- muL* L3[i];
                     }

                     //LARVAE 4
                     DL4[0] =   (nL3 * alphaL3 * L3[nL3-1]) +
                     (nDL4 * alphaDL * DIAL4[nDL4-1]) - 
                     (nL4 * alphaL4 * L4[0]) -
                     (DIA_INDUC_L4 * L4[0])-
                     (muL * L4[0]);

                     for ( i = 1;  i <nL4; i++){
                     DL4[i]=  nL4 * alphaL4 * (L4[i-1] - L4[i]) - (DIA_INDUC_L4*L4[i]) - 
                     (muL* L4[i]);
                     }

                     //LARVAE 5
                     DL5[0] = (nL4 * alphaL4 * L4[nL4-1]) +
                     (nDL5 * alphaDL * DIAL5[nDL5-1]) -           
                     (nL5 * alphaL5 * L5[0])  -
                      (DIA_INDUC_L5 * L5[0]) - 
                     (muL*L5[0]);

                     for ( i = 1;  i < nL5; i++){
                     DL5[i]=  nL5 * alphaL5 * (L5[i-1] - L5[i]) - 
                    (DIA_INDUC_L5*L5[i]) - 
                    (muL*L5[i]);
                     }

                     //DIAPAUSING LARVAL STAGE 4 

                     DDIAL4[0] =  DIA_INDUC_L4 * sum(L4,nL4) -
                     (nDL4 * alphaDL* DIAL4[0])-(muDL *DIAL4[0]);

                      for ( i = 1;  i <nDL4; i++){
                     DDIAL4[i] =  nDL4 * alphaDL *(DIAL4[i-1] -DIAL4[i]) - muDL * DIAL4[i];
                     }

                     //DIAPAUSING LARVAL STAGE 5 
                     DDIAL5[0] = DIA_INDUC_L5*sum(L5,nL5)  - 
                     (nDL5 * alphaDL*DIAL5[0])-  (muDL * DIAL5[0]);

                     for ( i = 1;  i <  nDL5; i++){
                     DDIAL5[i]=  nDL5* alphaDL*(DIAL5[i-1] - DIAL5[i]) - muDL * DIAL5[i];
                     }

                     //PUPAE 

                     DP[0] = nL5 * alphaL5 * L5[nL5-1] - nP * alphaP *  P[0] - muP * P[0] ;
                     for ( i = 1;  i <nP; i++){
                     DP[i] =  nP * alphaP * (P[i-1] - P[i]) - muP*P[i];
                     }

                     //REPRODUCTIVE ADULT

                     DAr[0] = nP * alphaP * P[nP-1] -nAr * alphaA*  Ar[0] - muA * Ar[0] ;

                      for ( i = 1;  i < nAr; i++){
                     DAr[i] =  nAr * alphaA * (Ar[i-1] - Ar[i]) - muA*Ar[i];
                      }

                     //SENESCENT ADULT
                     DAs =  nAr * alphaA * Ar[nAr-1] - muA * As;
                     ")

###This initializes the function. I wrote this is as an R code instead of a Csnippet
init <-function(params,t0,...)
  {
  c(E = rep(0, nE),
  L1 = rep(0,nL1),
  L2 = rep(0,nL2),
  L3 = rep(0,nL3),
  L4 = rep(0,nL4),
  L5 = rep(0,nL5),
  DIAL4 =rep(2,nDL4),
  DIAL5 =rep(0,nDL5),
  P = rep(0,nP),
  Ar = rep(0,nAr),
  As = 0)}

POMP_CM<- pomp(covartable[,c(1,4)],
         times="time",
         t0=1,
         globals=globs,
         initializer=init,
         covar = covartable ,
         tcovar = "time",
         skeleton = vectorfield(DIA_SKEL),
         statenames=c(sprintf("E%1d",seq_len(nE)),
                      sprintf("L1%1d",seq_len(nL1)),
                      sprintf("L2%1d",seq_len(nL2)),
                      sprintf("L3%1d",seq_len(nL3)),
                      sprintf("L4%1d",seq_len(nL4)),
                      sprintf("L5%1d",seq_len(nL5)),
                      sprintf("DIAL4%1d",seq_len(nDL4)),
                      sprintf("DIAL5%1d",seq_len(nDL5)),
                      sprintf("P%1d",seq_len(nP)),
                      sprintf("Ar%1d",seq_len(nAr)),
                      "As"))

###Make a dummy parameter, because I put all the parameter 
###values unto the extra code for debugging easier-ness
params = c('alphaTmp'=1.000)

###This runs the desolver and we integrate
TRAJ_MODEL <- trajectory(POMP_CM,params,times=seq(1,365*5),as.data.frame=TRUE)

###To get the total number of individuals within a life-stage, 
###We have to sum up the whole subcompartments 

EGG<-rowSums(TRAJ_MODEL [,1:(nE)])

LARVAE1<- rowSums(TRAJ_MODEL [,(nE+1): (nE+nL1)])

LARVAE2<- rowSums(TRAJ_MODEL[,(nE+nL1+1): (nE+nL1+nL2)])
LARVAE3<- rowSums(TRAJ_MODEL[,(nE+nL1+nL2+1): (nE+nL1+nL2+nL3)])
LARVAE4<- rowSums(TRAJ_MODEL[,(nE+nL1+nL2+nL3+1): (nE+nL1+nL2+nL3+nL4)])
LARVAE5<- rowSums(TRAJ_MODEL[,(nE+nL1+nL2+nL3+nL4+1): (nE+nL1+nL2+nL3+nL4+nL5)])

DIAL4<-  rowSums(TRAJ_MODEL[,(nE+nL1+nL2+nL3+nL4+nL5+1): (nE+nL1+nL2+nL3+nL4+nL5+nDL4)])

DIAL5<- rowSums(TRAJ_MODEL[,(nE+nL1+nL2+nL3+nL4+nL5+nDL4+1): (nE+nL1+nL2+nL3+nL4+nL5+nDL4+nDL5)])

Pupae<- rowSums(TRAJ_MODEL[,(nE+nL1+nL2+nL3+nL4+nL5+nDL4+nDL5+1): (nE+nL1+nL2+nL3+nL4+nL5+nDL4+nDL5+nP)])

R.adult<- rowSums(TRAJ_MODEL[,(nE+nL1+nL2+nL3+nL4+nL5+nDL4+nDL5+nP+1): (nE+nL1+nL2+nL3+nL4+nL5+nDL4+nDL5+nP+nAr)])

plot(EGG,col='red',type='l')
plot(LARVAE1,col='red')
plot(LARVAE2,col='red')
plot(LARVAE3,col='red')
plot(LARVAE4,col='red')
plot(LARVAE5,col='red')
plot(DIAL4,col='red')
plot(DIAL5,col='red')
plot(Pupae,col='red')
plot(R.adult,type='l',col='red')

Here is the output of the reproductive adults and you can see that the numbers are quite different than from the R desolve model. image

I know it's bad coding, but I put the values of the parameter directly into the model as I was trying to figure out what was happening. All models should be relying on the same photoperiod/temperature and should have the same initial starting values with the same number of subcompartments (8 across all stages).

I did some experimenting to figure out where the problem may be coming from. I decided to only run for 365 days and I set all the mortality rates to be 0, the two outputs were very similar (Figure below with black as original desolve and red as pomp). I set the mortality rate to a constant of 0.05 and it was the same case.

image

However, when I add the effect of temperature as (0.05 * TT) for all life stages in both the desolve and pomp code, that's when there's a huge discrepancy. The patterns seem to match up but there seem to be a big difference in the abundance. image

I have no idea why this is happening: but my hunch is that it's something to do with the mortality functions and with the covariates?

Here is the github-link- just in case more is needed. (https://github.com/pakdamie/modeling_diapause_CM)

Apologies for the somewhat messy code Thank you!

kingaa commented 5 years ago

@pakdamie, I won't have the time to read through your very thoroughly documented and carefully constructed, but extremely long and detailed, help request. Sorry for that.