R/IPMpack-Hidden.r

Defines functions createFecObj createSurvObj createGrowthObj makeListFmatrix makeListPmatrix alteredFit deathDataAugment stochLifeExpect sensParamsDiscrete predictMuX fecRaw fecPreCensusInteger fecPreCensus fecPostCensusInteger offspringCum fecPostCensus convergeLifeExpectancy convergeR0 convergeLambda makeListIPMs makeCovDf plotResultsStochStruct dentifyPossibleYearsCarlina getListRegObjectsFec getListRegObjects getIPMoutputDirect getIPMoutput generateDataStoch generateDataDiscrete

# TODO: Add comment
# 
# Author: ctg
###############################################################################
# This file contains functions that are hidden because they are no longer supported or because we like hiding them.

# =============================================================================
# =============================================================================
## Generate a simple data-frame for continuous and discrete covariates
#  with a total of 1000 measurements in columns called
# size, sizeNext, surv, fec, stage, stageNext number
# Stage contains names including "continuous", and then a range
# of names for discrete stages, e.g., in this example,
#  "dormant" "seedAge1"   "seedOld" 
#
# 
.generateDataDiscrete <- function(nSamp=1000){
  size <- rnorm(nSamp,5,2)
  sizeNext <- 1+0.8*size+rnorm(nSamp,0,1)
  surv <- rbinom(nSamp,1,invLogit(-1+0.2*size))
  sizeNext[surv==0] <- NA
  fec <- rnorm(length(size),exp(-7+0.9*size),1)
  fec[size<quantile(size,0.20) | fec<0] <- 0
  stage <- rep("continuous",nSamp)
  stageNext <- rep("continuous",nSamp)
  sizeNext[surv==0] <- NA
  stageNext[surv==0] <- c("dead")
  number <- rep(1,nSamp)
  become.dormant <- which(rank(size)%in%sample(rank(size),50,prob=surv*fec))
  sizeNext[become.dormant] <- NA; stageNext[become.dormant] <- c("dormant")
  were.dormant <- which(rank(sizeNext)%in%sample(rank(sizeNext),50,prob=surv*fec))
  size[were.dormant] <- NA; stage[were.dormant] <- c("dormant")
  dataf <- rbind(data.frame(size=size,sizeNext=sizeNext,surv=surv,
          fec=fec,stage=stage,stageNext=stageNext,number=number),
      data.frame(size=NA,sizeNext=NA,surv=rep(c(1,0),2),fec=0,
          stage=rep(c("seedAge1","seedOld"),each=2),
          stageNext=rep(c("seedOld","dead"),2),
          number=c(202,220,115,121)),
      data.frame(size=NA,sizeNext=rnorm(113,3,2),surv=1,fec=0,
          stage=c(rep("seedAge1",33),rep("seedOld",30),rep(NA,50)),
          stageNext=c("continuous"),number=1))
  
  
  return(dataf)
}

# =============================================================================
# =============================================================================
## Generate a simple data-frame for continuous and discrete covariates
#  with a total of 1000 measurements in columns called
# size, sizeNext, surv, fec, stage, stageNext number
#
# 
.generateDataStoch <- function(nSamp=1000){
  covariate1 <- rnorm(nSamp)
  covariate2 <- rnorm(nSamp)
  covariate3 <- rnorm(nSamp)
  size <- rnorm(nSamp,5,2)
  sizeNext <- 1+0.9*size+3*covariate1+0.01*covariate2+0.2*covariate3+rnorm(nSamp,0,0.1)
  
  fec <- surv <- rep(NA, length(size))
  surv[!is.na(size)] <- rbinom(sum(!is.na(size)),1,invLogit(-1+0.2*size[!is.na(size)]))
  fec[!is.na(size)] <- rnorm(sum(!is.na(size)),exp(-7+0.9*size[!is.na(size)]),1)
  fec[size<quantile(size,0.20,na.rm=TRUE) | fec<0] <- 0
  fec <- 10*fec
  
  seedlings <- sample(1:nSamp,size=100,replace=TRUE)
  size[seedlings] <- NA; 
  sizeNext[seedlings] <- rnorm(100,-2,0.1)
  surv[seedlings] <- 1
  #set to flower when covariate1 is around 1.5
  pfec <- 1*(runif(length(size))<invLogit(size+covariate1)); #print(pfec)
  fec[pfec==0] <- 0
  #fill in stage
  stage <- stageNext <- rep("continuous",nSamp)
  stage[is.na(size)] <- NA
  stageNext[is.na(sizeNext)] <- "dead"
  
  dataf <- data.frame(size=size,sizeNext=sizeNext,surv=surv,
      covariate1=covariate1,covariate2=covariate2,covariate3=covariate3,
      fec=fec, stage=stage,stageNext=stageNext, number=rep(1,length(size)))
  
  dataf$sizeNext[dataf$surv==0] <- NA
  
  return(dataf)
}

# =============================================================================
# =============================================================================
# Function to extract IPM output from a list
# of P (survival + growth) and F (fecundity) matrices
# (usually from Bayes fit) 
#
# Parameters - PmatrixList
#            - targetSize - the size you want passage time estimated for.
#            - FmatrixList
#
# Returns - a list 

.getIPMoutput <- function(PmatrixList,targetSize=c(),FmatrixList=NULL){
  
  if (length(targetSize)==0)  { 
    print("no target size for passage time provided; taking meshpoint median")
    targetSize <- median(PmatrixList[[1]]@meshpoints)
  }
  nsamps <- length(PmatrixList)
  h1 <- PmatrixList[[1]]@meshpoints[2]-PmatrixList[[1]]@meshpoints[1]
  stableStage <- LE <- pTime <- matrix(NA,nsamps,length(PmatrixList[[1]]@.Data[1,]))
  lambda <- rep(NA,nsamps)
  for (k in 1:nsamps) {
    Pmatrix <- PmatrixList[[k]]
    LE[k,]<-meanLifeExpect(Pmatrix) 
    pTime[k,]<-passageTime(targetSize,Pmatrix) 
    
    if (class(FmatrixList)!="NULL") {
      IPM <- Pmatrix + FmatrixList[[k]]
      lambda[k] <- Re(eigen(IPM)$value[1])
      stableStage[k,] <- eigen(IPM)$vector[,1]
      #normalize stable size distribution
      stableStage[k,] <- stableStage[k,]/(h1*sum(stableStage[k,]))
    }
  }
  
  return(list(LE=LE,pTime=pTime,lambda=lambda,stableStage=stableStage))
  
}

# =============================================================================
# =============================================================================
# Function to extract IPM output from posteriors 
# (usually from Bayes fit)  - quicker to do all at once
# rather than build list of IPM P matrices, then list of IPM F matrices
#
# Parameters - survObjlist - list of survival objects
#            - growObjList - list of growth objects
#            - targetSize - the size you want passage time estimated for.
#            - nBigMatrix - the number of bins
#            - minSize - the minimum size
#            - maxSize - the maximum size
#            - cov - do you want to fit a discrete covariate
#            - fecObjList - list of fecundity objects
#            - envMat - matrix of env transitions (only if cov=TRUE)
#            - nsizeToAge - numeric describing how many size to age defined (0 - 100s)
# 
#
# Returns - a list 

.getIPMoutputDirect <- function(survObjList,growObjList,targetSize=c(),
    nBigMatrix,minSize,maxSize,discreteTrans = 1,
    cov=FALSE,fecObjList=NULL, envMat=NULL,
    nsizeToAge=0, sizeStart = 10,
    integrateType = "midpoint", correction = "none", storePar=TRUE,
    chosenCov = data.frame(covariate = 1),
    onlyLowerTriGrowth=FALSE){
  
  # adjust the sample lengths to they are all the same
  if (length(targetSize)==0)  targetSize <- 0.2*(minSize+maxSize)
  nsamp <- max(length(growObjList),length(survObjList),length(fecObjList))
  if (length(survObjList)<nsamp)  
    survObjList <- sample(survObjList,size=nsamp,replace=TRUE)
  if (length(growObjList)<nsamp)  
    growObjList <- sample(growObjList,size=nsamp,replace=TRUE)
  if (class(fecObjList)!="NULL") {
    if (length(fecObjList)<nsamp)  
      fecObjList <- sample(fecObjList,size=nsamp,replace=TRUE)
  }
  
  # store chosen parameters
  if (storePar){
    surv.par <- matrix(NA,nsamp,length(survObjList[[1]]@fit$coefficients))
    grow.par <- matrix(NA,nsamp,length(growObjList[[1]]@fit$coefficients)+1)
    for (k in 1:nsamp) {
      surv.par[k,] <- survObjList[[k]]@fit$coefficients
      grow.par[k,] <- c(growObjList[[k]]@fit$coefficients,growObjList[[k]]@sd)
    }} else { surv.par <- grow.par <- c()}
  
  #set up storage
  if (class(discreteTrans)=="discreteTrans") nDisc <- (ncol(discreteTrans@discreteTrans)-1) else nDisc <- 0
  
  if (class(envMat)!="NULL") nEnv <- envMat@nEnvClass else nEnv <- 1
  LE <- pTime <- matrix(NA,nsamp,(nBigMatrix+nDisc)*nEnv)
  if (class(fecObjList)=="NULL") {
    lambda <- stableStage <- c()
  } else {
    stableStage <- matrix(NA,nsamp,(nBigMatrix+nDisc)*nEnv)
    lambda <- rep(NA,nsamp)
  }
  if (nsizeToAge==0) { resAge <- resSize <- c() } else { resAge <- resSize <- matrix(NA,nsamp,nsizeToAge)} 
  if (length(sizeStart)==0) { if (minSize<0) sizeStart <- 0.5*minSize else sizeStart <- 2*minSize }
  
  #go!
  for (k in 1:nsamp) {
    
    if (!cov) {
      Pmatrix <- makeIPMPmatrix(nBigMatrix = nBigMatrix, minSize = minSize, 
          maxSize = maxSize,  growObj = growObjList[[k]],
          survObj = survObjList[[k]],discreteTrans=discreteTrans,
          integrateType=integrateType, correction=correction) 
      
    } else {
      Pmatrix <- makeCompoundPmatrix(nEnvClass = nEnv,
          nBigMatrix = nBigMatrix, minSize = minSize,
          maxSize = maxSize, envMatrix=envMat,growObj = growObjList[[k]],
          survObj = survObjList[[k]],discreteTrans=discreteTrans,
          integrateType=integrateType, correction=correction)    
      
    }
    
    if (onlyLowerTriGrowth & !cov) {
      Pmatrix@.Data <- Pmatrix@.Data*lower.tri(Pmatrix@.Data, diag = TRUE)
      nvals <- colSums(Pmatrix@.Data,na.rm=TRUE)
      Pmatrix@.Data <- t((t(Pmatrix@.Data)/nvals) *
              surv(size = Pmatrix@meshpoints, 
                  cov = chosenCov,
                  survObj = survObjList[[k]]))                                    
    }
    
    
    LE[k,] <- meanLifeExpect(Pmatrix) 
    pTime[k,] <- passageTime(targetSize,Pmatrix) 
    if (k==1) h1 <- diff(Pmatrix@meshpoints)[1]
    
    if (class(fecObjList)!="NULL") {
      if (!cov) { 
        Fmatrix <- makeIPMFmatrix(nBigMatrix = nBigMatrix, minSize = minSize, 
            maxSize = maxSize,  
            fecObj=fecObjList[[k]],
            integrateType=integrateType, correction=correction)
      } else {
        Fmatrix <- makeCompoundFmatrix(nEnvClass = nEnv,
            nBigMatrix = nBigMatrix, minSize = minSize, 
            maxSize = maxSize, envMatrix=envMat,
            fecObj=fecObjList[[k]],integrateType=integrateType, correction=correction)
      }
      
      
      
      IPM <- Pmatrix + Fmatrix
      lambda[k] <- Re(eigen(IPM)$value[1])
      stableStage[k,] <- eigen(IPM)$vector[,1]
      #normalize stable size distribution
      stableStage[k,] <- stableStage[k,]/(h1*sum(stableStage[k,]))
      
      #print("here2")
    }
    
    # get size to age results
    if (nsizeToAge>0) { 
      res2 <- sizeToAge(Pmatrix=Pmatrix,startingSize=minSize*1.3,
          targetSize=seq(sizeStart,maxSize*0.9,length=nsizeToAge))
      resAge[k,] <- res2$timeInYears
      resSize[k,] <- res2$targetSize
    }
    
  }
  
  return(list(LE=LE,pTime=pTime,lambda=lambda,stableStage=stableStage,
          meshpoints=Pmatrix@meshpoints,resAge=resAge,resSize=resSize,
          surv.par=surv.par,grow.par=grow.par))
  
}

# =============================================================================
# =============================================================================
## Function to take fit of these and output a list of growth objects
.getListRegObjects <- function(Obj,nsamp=1000) {
  
  require(mvtnorm)
  
  #generate new set parameters from mvn
  npar <- length(Obj@fit$coefficients)
  newpar <- rmvnorm(nsamp, mean = Obj@fit$coefficients, sigma = vcov(Obj@fit))
  
  objList <- list()
  
  for (j in 1:nsamp) {
    objList[[j]] <- Obj
    objList[[j]]@fit$coefficients <- newpar[j,]
  }
  
  return(objList)
}


# =============================================================================
# =============================================================================
## Function to take fit of these and output a list of growth objects
.getListRegObjectsFec <- function(Obj,nsamp=1000) {
  
  require(mvtnorm)
  
  objList <- list()
  
  #generate new set parameters from mvn
  for (j in 1:nsamp) {
    for (k in 1:length(Obj@fitFec)) { 
      npar <- length(Obj@fitFec[[k]]$coefficients)
      newpar <- rmvnorm(nsamp, mean = Obj@fitFec[[k]]$coefficients, 
          sigma = vcov(Obj@fitFec[[k]]))
      objList[[j]] <- Obj
      objList[[j]]@fitFec[[k]]$coefficients <- newpar[j,]
    }
  }	
  
  return(objList)
}

# =============================================================================
# =============================================================================

#Find years where can estimate all three stochastic vital rates(survival, growth and baby size)
.identifyPossibleYearsCarlina <- function(dataf){
  
  yr1 <- table(dataf$year[!is.na(dataf$size) & 
              !is.na(dataf$sizeNext) & is.na(dataf$offspringNext)])
  yr2 <- table(dataf$year[!is.na(dataf$size) & 
              !is.na(dataf$surv) & is.na(dataf$offspringNext)])
  yr3 <- table(dataf$year[!is.na(dataf$sizeNext) & 
              !is.na(dataf$offspringNext)])
  
  good.yrs <- intersect(as.numeric(as.character(names(yr1)[yr1>1])),
      as.numeric(as.character(names(yr2))[yr2>1]))
  good.yrs <- intersect(good.yrs,as.numeric(as.character(names(yr3)[yr3>1])))
  
  return(is.element(dataf$year,good.yrs))
}

# =============================================================================
# =============================================================================
# Function to plot the results of a stochastic simulation
# structure run
#
# Parameters - tVals - time points
#            - st - output of stochGrowthRateManyCov
#            - covTest - the key covariate for germination / flowering
#            - nRunIn - how many to leave off pics
# Returns - 
#
.plotResultsStochStruct <- function(tVals,meshpoints,st,covTest,nRunIn=15,log="y",...) { 
  
  par(mfrow=c(1,3),bty="l", pty="s")
  plot(tVals[nRunIn:length(tVals)],
      colSums(st[,nRunIn:length(tVals)]),
      xlab="Time", 
      ylab="Population size",type="l",...)
  abline(v=1:max(tVals),lty=3)
  
  for (j in 1:ncol(covTest)) {
    #normalized
    covTest[,j] <- (covTest[,j]-mean(covTest[,j]))/sd(covTest[,j])
  }
  
  matplot(tVals[nRunIn:length(tVals)],covTest[nRunIn:length(tVals),],
      type="l",lty=3,col=1:ncol(covTest),xlab="Time", ylab="Covariates")
  abline(v=1:max(tVals),lty=3)
  
  
  if (log=="y") st <- log(st)
  
  image(tVals[nRunIn:length(tVals)],
      meshpoints,
      t(st[,nRunIn:length(tVals)]),
      ylab="Size", xlab="Time")
}


# =============================================================================
# =============================================================================
# makeCovDf creates a dataframe of size variables for prediction
# TODO: make able to use 'covariates'
.makeCovDf <- function(size, explanatoryVariables) {
  sizeSorted <- sort(size)
  splitVars <- strsplit(explanatoryVariables, split = "\\+")
  expVar <- unlist(splitVars[grep("size", splitVars)])
  covDf <- data.frame(size = sizeSorted)
  for(i in 1:length(expVar)) {
    if(is.null(expVar[i])) next
    expVar[i] <- sub('[[:space:]]', '', expVar[i])
    if(expVar[i] == "size2") {
      covDf$size2 <- sizeSorted ^ 2
    }
    if(expVar[i] == "size3"){
      covDf$size3 <- sizeSorted ^ 3
    }
    if(expVar[i] == "expsize") {
      covDf$expsize <- exp(sizeSorted)
    }
    if(expVar[i] == "logsize") {
      covDf$logsize <- log(sizeSorted)
    }
    if(expVar[i] == "logsize2") {
      covDf$logsize2 <- log(sizeSorted ^ 2)
    }
  }
  return(covDf)
}

# =============================================================================
# =============================================================================
.makeListIPMs<- function(dataf, nBigMatrix=10, minSize=-2,maxSize=10, 
    integrateType="midpoint", correction="none",
    explSurv=surv~size+size2+covariates,
    explGrow=sizeNext~size+size2+covariates, 
    regType="constantVar",explFec=fec~size,Family="gaussian", 
    Transform="none",fecConstants=data.frame(NA)) {
  
  #convert to 1:n for indexing later
  dataf$covariates <- as.factor(dataf$covariates)
  levels(dataf$covariates) <- 1:length(unique(dataf$covariates))
  
  print(explSurv)
  sv1 <- makeSurvObj(dataf=dataf,
      Formula=explSurv)
  gr1 <- makeGrowthObj(dataf=dataf,
      Formula=explGrow,
      regType=regType)
  
  fv1 <- makeFecObj(dataf=dataf,Formula=explFec, Family=Family, Transform=Transform, 
      fecConstants=fecConstants) 
  
  covs <- unique(dataf$covariates)
  covs <- covs[!is.na(covs)]
  
  #print(covs)
  
  IPM.list <- list()
  for (k in 1:length(covs)) { 
    
    tpF <- makeIPMFmatrix(nBigMatrix = nBigMatrix, minSize = minSize,
        maxSize = maxSize, chosenCov = data.frame(covariates=as.factor(k)),
        fecObj = fv1,integrateType=integrateType, correction=correction)
    tpS <- makeIPMPmatrix(nBigMatrix = nBigMatrix, minSize = minSize,
        maxSize = maxSize, chosenCov = data.frame(covariates=as.factor(k)),
        growObj = gr1, survObj = sv1,
        integrateType=integrateType, correction=correction)
    IPM.list[[k]] <- tpF+tpS
  }
  return(IPM.list)
  
}

# =============================================================================
# =============================================================================
.convergeLambda<-function(growObj, survObj, fecObj, nBigMatrix, minSize, maxSize, 
    discreteTrans = 1, integrateType = "midpoint", correction = "none", preCensus = TRUE, tol=1e-4,binIncrease=5){
  
  lambda.new<-1000
  delta<-1000
  
  print(paste(c("delta: ","new lambda:", "New number of grid points:
                  "),collapse=""))
  
  
  while(delta>tol) {
    lambda.old <-lambda.new
    nBigMatrix <- nBigMatrix + binIncrease
    Pmatrix <- makeIPMPmatrix(nBigMatrix = nBigMatrix, minSize = minSize, 
        maxSize = maxSize, growObj = growObj, survObj = survObj, 
        discreteTrans = discreteTrans, integrateType = integrateType, 
        correction = correction)
    Fmatrix <- makeIPMFmatrix(nBigMatrix = nBigMatrix, minSize = minSize, 
        maxSize = maxSize, fecObj = fecObj, integrateType = integrateType, 
        correction = correction, preCensus = preCensus,survObj=survObj,growObj=growObj)
    
    IPM <- Pmatrix + Fmatrix
    lambda.new <- Re(eigen(IPM)$value[1])
    
    delta<-abs(lambda.new-lambda.old)
    print(paste(c( delta, lambda.new,nBigMatrix),collapse=" "))
    
  }
  
  print(c("Final lambda from iteration:",lambda.new))
  print(c("Number of bins:",nBigMatrix))
  
  output <- list(nBigMatrix = nBigMatrix, IPM = IPM, lambda = lambda.new)
  
  return(output)
}


# =============================================================================
# =============================================================================
.convergeR0<-function(growObj, survObj, fecObj, nBigMatrix, minSize, maxSize, 
    discreteTrans = 1, integrateType = "midpoint", correction = "none", preCensus = TRUE, tol=1e-4,binIncrease=5){
  
  R0.new<-1000
  delta<-1000
  
  
  print(paste(c("delta: ","new R0:", "New number of grid points:"),collapse=""))
  
  while(delta>tol) {
    R0.old <-R0.new
    nBigMatrix <- nBigMatrix + binIncrease
    Pmatrix <- makeIPMPmatrix(nBigMatrix = nBigMatrix, minSize = minSize, 
        maxSize = maxSize, growObj = growObj, survObj = survObj, 
        discreteTrans = discreteTrans, integrateType = integrateType, 
        correction = correction)
    Fmatrix <- makeIPMFmatrix(nBigMatrix = nBigMatrix, minSize = minSize, 
        maxSize = maxSize, fecObj = fecObj, integrateType = integrateType, 
        correction = correction, preCensus = preCensus,survObj=survObj,growObj=growObj)
    
    R0.new <- R0Calc(Pmatrix,Fmatrix)
    
    delta<-abs(R0.new-R0.old)
    print(paste(c( delta, R0.new,nBigMatrix),collapse=" "))
  }
  
  IPM <- Pmatrix + Fmatrix
  
  print(c("Final R0 from iteration:",R0.new))
  print(c("Number of bins:",nBigMatrix))
  
  output<-list(nBigMatrix = nBigMatrix, binIncrease=binIncrease,IPM=IPM,R0=R0.new)
  
  return(output)
}

# =============================================================================
# =============================================================================
.convergeLifeExpectancy <- function(growObj, survObj,nBigMatrix, minSize, maxSize, 
    discreteTrans = 1, integrateType = "midpoint", correction = "none", 
    tol=1e-1,binIncrease=5, chosenBin=1){
  
  LE.new<-1000
  delta<-1000
  
  print(paste(c("delta: ","new LE:", "New number of grid points:"),collapse=""))
  
  while(delta>tol) {
    LE.old <- LE.new
    nBigMatrix <- nBigMatrix + binIncrease
    Pmatrix <- makeIPMPmatrix(nBigMatrix = nBigMatrix, minSize = minSize, 
        maxSize = maxSize, growObj = growObj, survObj = survObj, 
        discreteTrans = discreteTrans, integrateType = integrateType, 
        correction = correction)
    
    LE.new <- meanLifeExpect(Pmatrix)    
    
    delta <- abs(LE.new[chosenBin]-LE.old[chosenBin])
    print(paste(c( delta, LE.new[chosenBin],nBigMatrix),collapse=" "))
  }
  
  print(c("Final life expectancy of chosen bin from iteration:",LE.new[1]))
  print(c("Number of bins:",nBigMatrix))
  
  output<-list(nBigMatrix = nBigMatrix,
      binIncrease=binIncrease,Pmatrix=Pmatrix,LE=LE.new)
  
  return(output)
}

# =============================================================================
# =============================================================================
# TO DO .fecPostCensus properly (i.e. include size changes between the census and reproduction event) the following steps have to be made:
# 1. use growSurv to determine what the distribution of size is at the reproduction event given initial size x
# 2. over this distribution of sizes x2 at the reproduction event, what are the expected number of offspring: per x2: .fecRaw(x=x2,...)
# 3. multiply the expected number of offspring per x2 with the probability that an offspring is of size y, using dnorm(y,predict(..., newdata=newd2, ...) where newd2 is calculated for all levels of x2
# all in all, Eelke wonders if the outer-solution is still useful in the .fecPostCensus case, since x2 needs to have a certain distribution, which is not passed down to .fecPostCensus... 

## REMOVE GROWTH FROM THIS - note that this means 
#### growth obj generally not needed down below.....
## A function that outer can use showing numbers from x to y via production, growth, survival and distribution offspring
.fecPostCensus <- function(x,y,cov=data.frame(covariate=1),
    fecObj, growObj,survObj,offspringObj=NULL) {
  newd <- data.frame(cbind(cov,size=x),
      stringsAsFactors = FALSE)
  
  newd$size2 <- x^2
  newd$size3 <- x^3
  
  if (is.null(offspringObj)){
    if (length(grep("expsize",fecObj@offspringRel$formula))>0 |
        length(grep("expsize",growObj@fit$formula))>0) { newd$expsize <- exp(x)}            	
    if (length(grep("logsize",fecObj@offspringRel$formula))>0 |
        length(grep("logsize",growObj@fit$formula))>0) { newd$logsize <- log(x)}            
    u <- .fecRaw(x=x,cov=cov,fecObj=fecObj)[[1]]*
        dnorm(y,predict(fecObj@offspringRel,newdata=newd, type="response"),fecObj@sdOffspringSize)*
        surv(size=x, cov=cov, survObj=survObj)
  } else {
    u <- .fecRaw(x=x,cov=cov,fecObj=fecObj)[[1]]*
        growth(y,y,newd,offspringObj)*surv(size=x, cov=cov, survObj=survObj)				
  }
  return(u)
}

# =============================================================================
# =============================================================================
## A function that outer can use giving pnorm for offspring reprod
.offspringCum <- function(x,y,cov=data.frame(covariate=1),fecObj,offspringObj=NULL) {
  newd <- data.frame(cbind(cov,size=x),
      stringsAsFactors = FALSE)
  
  if (is.null(offspringObj)){
    newd$size2 <- x^2
    newd$size3 <- x^3
    if (length(grep("expsize",fecObj@offspringRel$formula))>0) { newd$expsize <- exp(x)}            
    if (length(grep("logsize",fecObj@offspringRel$formula))>0) { newd$logsize <- log(x)}            
    u <- pnorm(y,predict(fecObj@offspringRel,newdata=newd, type="response"),fecObj@sdOffspringSize)
  } else { 
    u <- growthCum(y,y,newd, offspringObj)	
  }
  return(u)
}

# =============================================================================
# =============================================================================
#### growth obj generally not needed down below.....
## A function that outer can use showing numbers from x to y via production, growth, survival and distribution offspring
.fecPostCensusInteger <- function(x,y,cov=data.frame(covariate=1),fecObj, growObj,survObj, offspringObj=NULL) {
  newd <- data.frame(cbind(cov,size=x),
      stringsAsFactors = FALSE)
  
  u <- .fecRaw(x=x,cov=cov,fecObj=fecObj)[[1]]*
      surv(size=x, cov=cov, survObj=survObj)
  
  if (is.null(offspringObj)){
    newd$size2 <- x^2
    newd$size3 <- x^3
    if (length(grep("expsize",fecObj@offspringRel$formula))>0 |
        length(grep("expsize",growObj@fit$formula))>0) { newd$expsize <- exp(x)}            
    if (length(grep("logsize",fecObj@offspringRel$formula))>0 |
        length(grep("logsize",growObj@fit$formula))>0) { newd$logsize <- log(x)}            
    
    if (fecObj@distOffspring=="poisson")
      u <- u*dpois(y,predict(fecObj@offspringRel,newdata=newd, type="response"))
    if (fecObj@distOffspring=="negBin")
      u <- u*dnbinom(y,mu=predict(fecObj@offspringRel,newdata=newd, type="response"),
          size=fecObj@thetaOffspringSize)
  } else {
    u <- u*growth(y,y,newd, offspringObj)
    
  }
  
  
  return(u)
}
# =============================================================================
# =============================================================================

## A function that outer can use showing numbers from x to y via production and distribution offspring
.fecPreCensus <- function(x,y,cov=data.frame(covariate=1),fecObj,offspringObj=NULL) {
  newd <- data.frame(cbind(cov,size=x),
      stringsAsFactors = FALSE)	
  newd$size2 <- x^2
  newd$size3 <- x^3
  
  if (is.null(offspringObj)){
    if (length(grep("expsize",
            fecObj@offspringRel$formula))>0) { newd$expsize <- exp(x)}
    if (length(grep("logsize",
            fecObj@offspringRel$formula))>0) { newd$logsize <- log(x)}
    u <- .fecRaw(x=x,cov=cov,fecObj=fecObj)[[1]]*
        dnorm(y,predict(fecObj@offspringRel,newdata=newd, type="response"),
            fecObj@sdOffspringSize)
  } else {
    u <- .fecRaw(x=x,cov=cov,fecObj=fecObj)[[1]]*
        growth(y,y,newd,offspringObj)
    
  }
  #print(cbind(y,predict([email protected])))
  
  
  return(u)
}

# =============================================================================
# =============================================================================
## A function that outer can use showing numbers from x to y via production and distribution offspring
.fecPreCensusInteger <- function(x,y,cov=data.frame(covariate=1),fecObj,offspringObj=NULL) {
  newd <- data.frame(cbind(cov,size=x),
      stringsAsFactors = FALSE)	
  newd$size2 <- x^2
  newd$size3 <- x^3
  
  u <- .fecRaw(x=x,cov=cov,fecObj=fecObj)[[1]]
  
  if (is.null(offspringObj)){
    if (length(grep("expsize",
            fecObj@offspringRel$formula))>0) { newd$expsize <- exp(x)}
    if (length(grep("logsize",
            fecObj@offspringRel$formula))>0) { newd$logsize <- log(x)}
    if (fecObj@distOffspring=="poisson")
      u <- u*dpois(y,predict(fecObj@offspringRel,newdata=newd, type="response"))
    if (fecObj@distOffspring=="negBin")
      u <- u*dnbinom(y,mu=predict(fecObj@offspringRel,newdata=newd, type="response"),
          size=fecObj@thetaOffspringSize)
  } else {
    u <- u*growth(y,y,newd, offspringObj)
    
  }
  
  #print(cbind(y,predict([email protected])))
  
  return(u)
}

# =============================================================================
# =============================================================================
## Get raw numbers of offspring produced by every size class by multiplying up the constants,
## and doing all the "predict: values needed; and taking out only the babies that go to the continuous classes

.fecRaw <- function(x,cov=data.frame(covariate=1),fecObj) { 
  
  newd <- data.frame(cbind(cov,size=x),
      stringsAsFactors = FALSE)
  
  newd$size2 <- x^2
  newd$size3 <- x^3
  newd$expsize <- exp(x)
  
  #if (length([email protected])>1) {
  #	if (length(grep("logsize", [email protected]$formula))>0) { newd$logsize <- log(x)}
  #}
  
  fecObj@fecConstants[is.na(fecObj@fecConstants)] <- 1
  
  #fecundity rates
  fecValues <- matrix(c(rep(1,length(fecObj@fitFec)),unlist(fecObj@fecConstants)),
      ncol=length(x),nrow=length(fecObj@fitFec)+
          length(fecObj@fecConstants))
  #rownames(fecValues) <- c([email protected],names([email protected]))
  for (i in 1:length(fecObj@fitFec)) fecValues[i,] <- predict(fecObj@fitFec[[i]],newd,type="response")
  if (length(grep("log",fecObj@Transform))>0) for (i in grep("log",fecObj@Transform)) fecValues[i,]<-exp(fecValues[i,])
  if (length(grep("exp",fecObj@Transform))>0) for (i in grep("exp",fecObj@Transform)) fecValues[i,]<-log(fecValues[i,])
  if (length(grep("sqrt",fecObj@Transform))>0) for (i in grep("sqrt",fecObj@Transform)) fecValues[i,]<-(fecValues[i,])^2
  if (length(grep("-1",fecObj@Transform))>0) for (i in grep("-1",fecObj@Transform)) fecValues[i,]<-fecValues[i,]+1
  if (length(which(fecObj@vitalRatesPerOffspringType[,"continuous"]==1))>1) {
    prodFecValues <- apply(fecValues[which(fecObj@vitalRatesPerOffspringType[,"continuous"]==1),],2,prod)*unlist(fecObj@offspringSplitter["continuous"])
  } else {
    prodFecValues <- fecValues[which(fecObj@vitalRatesPerOffspringType[,"continuous"]==1),]*unlist(fecObj@offspringSplitter["continuous"])
  }
  return(list(prodFecValues,fecValues))
}

# =============================================================================
# =============================================================================

# function to predict growth object - NOTE THAT THIS will not work with factors / contrasts
.predictMuX <- function(grObj, newData, covPred = 1) {
  dataBase <- newData
  coefNames <- attr(grObj@fit$coefficients, "names")
  coefValues <- as.matrix(grObj@fit$coefficients)
  covNames <- coefNames[grep("covariate", coefNames)]
  covPos <- as.numeric(unlist(strsplit(covNames, "covariate")))
  covPos <- as.numeric(covPos[!is.na(covPos)])
  covDf <- as.data.frame(matrix(0, nrow = dim(newData)[1], ncol = length(covPos)))
  names(covDf) <- covNames
  # Turn "on" intended covariate and build newDataFrame
  if(covPred != 1) {
    covDf[, (covPred - 1)] <- 1
  }
  newData <- cbind(dataBase, covDf)
  newDataSubset <- as.matrix(cbind(1, newData[, (names(newData) %in% coefNames)]))
  predValues <- as.matrix(newDataSubset) %*% matrix(coefValues, length(coefValues), 1)
  return(as.numeric(predValues))
}

# =============================================================================
# =============================================================================
### sens params for discrete survival bit

.sensParamsDiscrete <-  function (growObj, survObj, fecObj, nBigMatrix, minSize, maxSize, 
    chosenCov = data.frame(covariate=1),		
    discreteTrans = 1, integrateType = "midpoint", correction = "none", preCensus = TRUE, delta=1e-4) {
  
  
  #all the transitions - note that this is in the order of the columns, 
  nmes <- paste("",c(outer(colnames(discreteTrans@discreteTrans),colnames(discreteTrans@discreteTrans),paste,sep="-")),sep="")
  
  # all survival out of discrete stages
  nmes <- c(nmes,paste("survival from",colnames(discreteTrans@discreteSurv),sep=""))
  
  # all means out of discrete stages
  nmes <- c(nmes,paste("mean from ",colnames(discreteTrans@meanToCont),sep=""))
  
  # all sds out of discrete stages
  nmes <- c(nmes,paste("sd from ",colnames(discreteTrans@sdToCont),sep=""))
  
  #  not sure makes sense to do distribToDiscrete???	
  
  #coefficients linking continuous survival into discrete
  nmes <- c(nmes, paste("survival from continuous ", names(discreteTrans@moveToDiscrete$coefficients),sep=""))
  
  
  elam <- rep(NA,length(nmes))
  names(elam) <- nmes
  slam <- elam
  Pmatrix <- makeIPMPmatrix(nBigMatrix = nBigMatrix, minSize = minSize, 
      maxSize = maxSize, growObj = growObj, survObj = survObj, 
      chosenCov = chosenCov,
      discreteTrans = discreteTrans, integrateType = integrateType, 
      correction = correction)
  Fmatrix <- makeIPMFmatrix(nBigMatrix = nBigMatrix, minSize = minSize, 
      chosenCov = chosenCov,
      maxSize = maxSize, fecObj = fecObj, integrateType = integrateType, 
      correction = correction,preCensus = preCensus,survObj=survObj,growObj=growObj)
  IPM <- Pmatrix + Fmatrix
  lambda1 <- Re(eigen(IPM)$value[1])
  
  #all the transitions
  param.test <- 0
  for (j in 1:nrow(discreteTrans@discreteTrans)) {
    for (k in 1:nrow(discreteTrans@discreteTrans)) {
      
      #print(c(rownames([email protected])[k],colnames([email protected])[j]))
      
      param.test <- param.test+1
      
      if (discreteTrans@discreteTrans[k,j]==0) next()
      
      #pick element of matrix
      adj <- rep(discreteTrans@discreteTrans[k,j]*delta,nrow(discreteTrans@discreteTrans)-1)
      discreteTrans@discreteTrans[k,j] <- discreteTrans@discreteTrans[k,j] * (1 + delta)
      #alter the other values in the columns so as continue to sum to one
      if (sum(discreteTrans@discreteTrans[k,-j]>0)>0) adj <- adj/sum(discreteTrans@discreteTrans[k,-j]>0)
      adj[discreteTrans@discreteTrans[-k,j]==0] <- 0
      discreteTrans@discreteTrans[-k,j] <- discreteTrans@discreteTrans[-k,j]-adj
      
      #print(colSums([email protected]))
      
      
      Pmatrix <- makeIPMPmatrix(nBigMatrix = nBigMatrix, 
          minSize = minSize, maxSize = maxSize, growObj = growObj, 
          chosenCov = chosenCov,
          survObj = survObj, discreteTrans = discreteTrans, 
          integrateType = integrateType, correction = correction)
      
      Fmatrix <- makeIPMFmatrix(nBigMatrix = nBigMatrix, 
          minSize = minSize, maxSize = maxSize, fecObj = fecObj, 
          chosenCov = chosenCov,
          integrateType = integrateType, correction = correction,
          preCensus = preCensus,survObj=survObj,growObj=growObj)
      
      IPM <- Pmatrix + Fmatrix
      lambda2 <- Re(eigen(IPM)$value[1])
      discreteTrans@discreteTrans[k,j] <-  discreteTrans@discreteTrans[k,j]/(1 + delta)	
      discreteTrans@discreteTrans[-k,j] <- discreteTrans@discreteTrans[-k,j]+adj
      
      slam[param.test] <- (lambda2 - lambda1)/(discreteTrans@discreteTrans[k,j] * delta)
      elam[param.test] <- (lambda2 - lambda1)/(lambda1*delta)
    }}
  
  
  # all survival out of discrete stages
  count <- param.test
  for (param.test in 1:length(discreteTrans@discreteSurv)) { 
    
    discreteTrans@discreteSurv[param.test] <- discreteTrans@discreteSurv[param.test]* (1 + delta)
    
    Pmatrix <- makeIPMPmatrix(nBigMatrix = nBigMatrix, 
        minSize = minSize, maxSize = maxSize, growObj = growObj, 
        survObj = survObj, discreteTrans = discreteTrans, 
        chosenCov = chosenCov,
        integrateType = integrateType, correction = correction)
    
    Fmatrix <- makeIPMFmatrix(nBigMatrix = nBigMatrix, 
        minSize = minSize, maxSize = maxSize, fecObj = fecObj, 
        chosenCov = chosenCov,
        integrateType = integrateType, correction = correction,
        preCensus = preCensus,survObj=survObj,growObj=growObj)
    
    IPM <- Pmatrix + Fmatrix
    lambda2 <- Re(eigen(IPM)$value[1])
    
    discreteTrans@discreteSurv[param.test] <- discreteTrans@discreteSurv[param.test]/ (1 + delta)
    slam[param.test+count] <- (lambda2 - lambda1)/(discreteTrans@discreteSurv[param.test] * delta)
    elam[param.test+count] <- (lambda2 - lambda1)/(lambda1*delta)
    
  }
  
  
  # all mean values coming out of discrete stages
  count <- param.test+count
  for (param.test in 1:length(discreteTrans@meanToCont)) { 
    
    discreteTrans@meanToCont[param.test] <- discreteTrans@meanToCont[param.test]* (1 + delta)
    
    Pmatrix <- makeIPMPmatrix(nBigMatrix = nBigMatrix, 
        minSize = minSize, maxSize = maxSize, growObj = growObj, 
        chosenCov = chosenCov,
        survObj = survObj, discreteTrans = discreteTrans, 
        integrateType = integrateType, correction = correction)
    
    Fmatrix <- makeIPMFmatrix(nBigMatrix = nBigMatrix, 
        minSize = minSize, maxSize = maxSize, fecObj = fecObj, 
        chosenCov = chosenCov,
        integrateType = integrateType, correction = correction,
        preCensus = preCensus,survObj=survObj,growObj=growObj)
    
    IPM <- Pmatrix + Fmatrix
    lambda2 <- Re(eigen(IPM)$value[1])
    
    discreteTrans@meanToCont[param.test] <- discreteTrans@meanToCont[param.test]/ (1 + delta)
    slam[param.test+count] <- (lambda2 - lambda1)/(discreteTrans@meanToCont[param.test] * delta)
    elam[param.test+count] <- (lambda2 - lambda1)/(lambda1*delta)
    
  }
  
  
  # all sd values coming out of discrete stages
  count <- param.test+count
  for (param.test in 1:length(discreteTrans@sdToCont)) { 
    
    discreteTrans@sdToCont[param.test] <- discreteTrans@sdToCont[param.test]* (1 + delta)
    
    Pmatrix <- makeIPMPmatrix(nBigMatrix = nBigMatrix, 
        minSize = minSize, maxSize = maxSize, growObj = growObj, 
        chosenCov = chosenCov,
        survObj = survObj, discreteTrans = discreteTrans, 
        integrateType = integrateType, correction = correction)
    
    Fmatrix <- makeIPMFmatrix(nBigMatrix = nBigMatrix, 
        minSize = minSize, maxSize = maxSize, fecObj = fecObj, 
        chosenCov = chosenCov,
        integrateType = integrateType, correction = correction,
        preCensus = preCensus,survObj=survObj,growObj=growObj)
    
    IPM <- Pmatrix + Fmatrix
    lambda2 <- Re(eigen(IPM)$value[1])
    
    discreteTrans@sdToCont[param.test] <- discreteTrans@sdToCont[param.test]/ (1 + delta)
    slam[param.test+count] <- (lambda2 - lambda1)/(discreteTrans@sdToCont[param.test][param.test] * delta)
    elam[param.test+count] <- (lambda2 - lambda1)/(lambda1*delta)
    
  }
  
  # parameters linking size to survival
  count <- param.test+count
  for (param.test in 1:length(discreteTrans@moveToDiscrete$coefficients)) { 
    
    discreteTrans@moveToDiscrete$coefficients[param.test] <- discreteTrans@moveToDiscrete$coefficients[param.test]* (1 + delta)
    
    Pmatrix <- makeIPMPmatrix(nBigMatrix = nBigMatrix, 
        minSize = minSize, maxSize = maxSize, growObj = growObj, 
        chosenCov = chosenCov,
        survObj = survObj, discreteTrans = discreteTrans, 
        integrateType = integrateType, correction = correction)
    
    Fmatrix <- makeIPMFmatrix(nBigMatrix = nBigMatrix, 
        minSize = minSize, maxSize = maxSize, fecObj = fecObj, 
        chosenCov = chosenCov,
        integrateType = integrateType, correction = correction,
        preCensus = preCensus,survObj=survObj,growObj=growObj)
    
    IPM <- Pmatrix + Fmatrix
    lambda2 <- Re(eigen(IPM)$value[1])
    
    discreteTrans@moveToDiscrete$coefficients[param.test] <- discreteTrans@moveToDiscrete$coefficients[param.test]/ (1 + delta)
    slam[param.test+count] <- (lambda2 - lambda1)/(discreteTrans@moveToDiscrete$coefficient[param.test] * delta)
    elam[param.test+count] <- (lambda2 - lambda1)/(lambda1*delta)
    
  }
  
  print("Did not calculate sensitivities and elasticities for:")
  print(c(names(slam[is.na(slam)])))
  print("Values of zero")
  
  slam <- slam[!is.na(slam)]
  elam <- elam[!is.na(elam)]
  
  return(list(slam = slam, elam = elam))
  
}

# =============================================================================
# =============================================================================
#Generic for life expectancy in a modeled stoch env  -NEEDS DOUBLE-CHECKING
#parameters - a compound IPM of dim nEnvClass*nBigMatrix
#           - an environmental matrix
# returns - the life expectancy for each of the sizes in the IPM (columns)
#           for each of the starting env states
#
#
.stochLifeExpect <- function(IPMmatrix,envMatrix){
  
  matrix.dim <- length(IPMmatrix[1,])
  nstages <- IPMmatrix@nBigMatrix
  nstates <- IPMmatrix@nEnvClass
  
  pis <- Re(eigen(envMatrix)$vector[,1])
  pis <- pis/(sum(pis))
  
  #print(pis)
  
  #ckron <- kronecker(envMatrix,diag(nstages))  #doesnt work??
  m <- IPMmatrix  #%*%ckron  #eq 29 in carols paper
  
  Itilda <- diag(matrix.dim)
  
  #not used in this one
  eatildas <- array(dim=c(nstates*nstages,nstages,nstates))
  eatildas[,,] <-0
  for (i in 1:nstates){
    eatildas[,,i] <- Itilda[,((i-1)*nstages+1):(nstages*i)]
  }
  
  qatildas <- array(dim=c(nstates*nstages,nstages,nstates));
  qatildas[,,]<-0
  for (i in 1:nstates) {
    indext <-((i-1)*nstages+1):(nstages*i) 
    qatildas[indext,,i] <- IPMmatrix[indext,indext]/envMatrix[i,i]
    #print( IPMmatrix[cbind(indext,indext)]/envMatrix[i,i] )
  }                            #array of qatildas, eqn 26
  #need to remove env effect since mega-matrix pre-built 
  
  #print(qatildas)
  
  etilda <- array(dim=c(nstates*nstages,nstages));
  etilda[,]<-0
  for (i in 1:nstates) { 
    etilda[((i-1)*nstages+1):(nstages*i),] <- diag(nstages);
  }                            #etilda, eqn 27
  
  I <- diag(nstages);                 #identity matrix of dimension K x K
  
  Ns.markov <- array(dim=c(nstages,nstages,nstates)); #set up for array of Ns
  Ns.markov[,,]<-0
  for (i in 1:nstates){ 
    Ns.markov[,,i] <- I + t(etilda)%*%(solve(Itilda-m))%*%qatildas[,,i];
    #eqn 28, conditional on initial state 1
  }
  
  #print(Ns.markov)
  
  #average over all initial states %%%%%
  Nbar.markov <- rep(0,nstages);
  for (i in 1:nstates) { 
    Nbar.markov <-  Nbar.markov + pis[i] * Ns.markov[,,i]; 
  }                         #eqn 29, weight each fundamntal matrix by expctd frequency 
  
  
  
  #lifeexp, column sums
  lifeexp.markov<-matrix(0,nstates,nstages); #set up array
  for (i in 1:nstates) { 
    lifeexp.markov[i,] <- apply(Ns.markov[,,i], 2, sum);
  }
  
  return(lifeexp.markov)
}

# =============================================================================
# =============================================================================
# =============================================================================
# =============================================================================
# Function to augment mortality data for large trees
# Parameters dataf - existing data frame
#            size.thresh - the size above which tree death is being augmented
#            prop.dead - the proportion of these expected to be dead

.deathDataAugment <- function (dataf, size.thresh, prop.dead) { 
  
  n.now <- sum(dataf$size > size.thresh)
  n.new.dead <- ceiling(prop.dead * n.now / (1 - prop.dead))
  new.size <- rnorm(n.new.dead,size.thresh, sd(dataf$size[dataf$size > size.thresh]))
  
  datanew <- data.frame(size = new.size, sizeNext = rep(NA, n.new.dead), surv = rep(0, n.new.dead), 
      covariate = rep(0, n.new.dead), covariateNext = rep(0, n.new.dead),
      fec = rep(NA, n.new.dead)) 
  
  dataf.new <- rbind(dataf,datanew)
  
  return(dataf.new)
  
}

# =============================================================================
# =============================================================================
# replace the growth object fit with a new, desired variance for predict
.alteredFit <- function(dummyFit = dummyFit, 
    newCoef = dummyFit$coefficients, 
    desiredSd = 1) {
  dummyFit$coefficients[] <- newCoef
  dummyFit$residuals <- rnorm(length(dummyFit$residuals), mean = 0, sd = desiredSd)	
  # need to use qr here to assign dummyFit so that there is no warning when decomposed for n
  # Error in rnorm(residDf, 0, sd = desiredSd) : object 'residDf' not found
  return(dummyFit)	
}

# =============================================================================
# =============================================================================
# Function to take a list of growth and survival objects and make a list of Pmatrices
#
# Parameters - growObjList - a list of growth objects
#            - survObjList - a list of survival objects
#            - nBigMatrix - the number of bins
#            - minSize - the minimum size
#            - maxSize - the maximum size
#            - cov - is a discrete covariate considered
#            - envMat - enviromental matrix for transition between
# 
# Returns    - a list of Pmatrices
.makeListPmatrix  <- function(growObjList,survObjList,
    nBigMatrix,minSize,maxSize, cov=FALSE, envMat=NULL,
    integrateType="midpoint",correction="none", discreteTransList=1) {
  
  if (length(growObjList)>length(survObjList)) { 
    survObjList <- sample(survObjList,size=length(growObjList),replace=TRUE)
  } else { 
    if (length(growObjList)<length(survObjList))  
      growObjList <- sample(growObjList,size=length(survObjList),replace=TRUE)
  }
  
  nsamp <- length(growObjList)
  
  if(length(discreteTransList)<nsamp ){
    # if(warn) warning('Length of discreteTrans list is less than the length of another vital rate object list, so some members of the discreteTrans list have been repeated.')
    discreteTransList <- sample(discreteTransList,size=nsamp,replace=T)
  }
  
  PmatrixList <- list()
  for ( k in 1:length(growObjList)) { 
    if (!cov) {
      PmatrixList[[k]] <- makeIPMPmatrix(nBigMatrix = nBigMatrix, minSize = minSize, 
          maxSize = maxSize, growObj = growObjList[[k]],discreteTrans=discreteTransList[[k]],
          survObj = survObjList[[k]],integrateType=integrateType, correction=correction) 
    } else {
      PmatrixList[[k]] <- makeCompoundPmatrix(nEnvClass = length(envMat[1,]),
          nBigMatrix = nBigMatrix, minSize = minSize, 
          maxSize = maxSize, envMatrix=envMat,discreteTrans=discreteTransList[[k]],
          growObj = growObjList[[k]],
          survObj = survObjList[[k]],integrateType=integrateType, correction=correction)    
    }
  }
  
  return(PmatrixList)
}

# =============================================================================
# =============================================================================
# Function to take a list of growth and survival objects and make a list of Fmatrices

.makeListFmatrix  <- function(fecObjList,nBigMatrix,minSize,maxSize, cov=FALSE, 
    envMat=NULL,integrateType="midpoint",correction="none") {
  
  nsamp <- max(length(fecObjList))
  if (length(fecObjList)<nsamp)  
    fecObjList <- sample(fecObjList,size=nsamp,replace=TRUE)
  
  FmatrixList <- list()
  for ( k in 1:nsamp) {
    if (!cov) { 
      FmatrixList[[k]] <- makeIPMFmatrix(nBigMatrix = nBigMatrix, minSize = minSize, 
          maxSize = maxSize, 
          fecObj=fecObjList[[k]],integrateType=integrateType, correction=correction)
    } else {
      FmatrixList[[k]] <- makeCompoundFmatrix(nEnvClass = length(envMat[1,]),
          nBigMatrix = nBigMatrix, minSize = minSize, 
          maxSize = maxSize, envMatrix=envMat,
          fecObj=fecObjList[[k]],integrateType=integrateType, correction=correction)
    }
    
    FmatrixList[[k]] <-  FmatrixList[[k]]
  }
  return(FmatrixList)
}


# =============================================================================
# =============================================================================
### new functions createGrowhtObj and createSurvObj which will make up their own data

.createGrowthObj <- function(Formula=sizeNext~size, coeff=c(1,1), sd=1){ 
  
  var.names <- all.vars(Formula)
  
  if (length(coeff)!=(length(var.names))) #length var.names, because intercept not named 
    stop("not enough coefficients supplied for the chosen Formula")
  
  dataf<- as.data.frame(matrix(rnorm(3*length(var.names)),3,length(var.names)))
  colnames(dataf) <- var.names
  
  fit <- lm(Formula, data=dataf)
  
  if (length(grep("sizeNext", as.character(Formula))) > 0) { 
    
    gr1 <- new("growthObj")
    gr1@fit <- fit
    gr1@fit$coefficients <- coeff
    gr1@sd <- sd
  }  
  
  if (length(grep("incr", as.character(Formula))) > 0) { 
    
    gr1 <- new("growthObjIncr")
    gr1@fit <- fit
    gr1@fit$coefficients <- coeff
    gr1@sd <- sd
  }  
  
  return(gr1)
  
}

# =============================================================================
# =============================================================================
.createSurvObj <- function(Formula=surv~size, coeff=c(1,1)){ 
  var.names <- all.vars(Formula)
  
  #not that although var.names will have one extra (cos of response variable
  # this will correspond to the intercept
  if (length(coeff)!=(length(var.names))) 
    stop("not enough coefficients supplied for the chosen Formula")
  
  dataf<- as.data.frame(matrix(rnorm(3*length(var.names)),3,length(var.names)))
  colnames(dataf) <- var.names
  dataf$surv <- sample(c(0,1),nrow(dataf), replace=TRUE)
  
  fit <- glm(Formula, data=dataf, family=binomial)
  
  sv1 <- new("survObj")
  sv1@fit <- fit
  
  return(sv1)
  
}

# =============================================================================
# =============================================================================
.createFecObj <- function(Formula=list(fec1~size,fec2~size+size2), 
    coeff=list(c(1,1),c(1,1,1)),
    Family = c("gaussian","binomial"),
    Transform = c("log","none"),
    meanOffspringSize = NA, sdOffspringSize = NA, 
    offspringSplitter = data.frame(continuous = 1), 
    vitalRatesPerOffspringType = data.frame(NA), 
    fecByDiscrete = data.frame(NA), 
    offspringSizeExplanatoryVariables = "1",
    fecConstants = data.frame(NA), 
    doOffspring=TRUE, 
    reproductionType="sexual"){ 
  var.names <- c()
  fecNames <- rep(NA,length(Formula))
  
  for (j in 1:length(Formula)) { 
    
    
    fecNames[j] <- all.vars(Formula[[j]])[1]
    var.names.here <- all.vars(Formula[[j]])
    
    if (length(coeff[[j]])!=(length(var.names.here))) 
      stop(paste("not enough coefficients supplied for the ",j, "th Formula", sep=""))
    
    var.names <- c(var.names,var.names.here)
  }
  
  var.names <- unique(var.names)
  
  #build a data-frame with all the right variables
  dataf<- as.data.frame(matrix(rnorm(3*length(var.names)),3,length(var.names)))
  colnames(dataf) <- var.names
  dataf$surv <- sample(c(0,1),nrow(dataf), replace=TRUE)
  
  if (sum(Transform=="log")>0) dataf[,fecNames[which(Transform=="log")]] <- pmax(dataf[,fecNames[which(Transform=="log")]],1)
  if (sum(Family=="binomial")>0) dataf[,fecNames[which(Family=="binomial")]] <- rbinom(nrow(dataf),1,0.5)
  
  fv1 <- makeFecObj(dataf=dataf, 
      fecConstants = fecConstants, 
      Formula = Formula, 
      Family = Family, 
      Transform = Transform, 
      meanOffspringSize = meanOffspringSize, 
      sdOffspringSize = sdOffspringSize, offspringSplitter = offspringSplitter, 
      vitalRatesPerOffspringType = vitalRatesPerOffspringType, 
      fecByDiscrete = fecByDiscrete, 
      offspringSizeExplanatoryVariables = offspringSizeExplanatoryVariables, 
      coeff=NULL, doOffspring=doOffspring, 
      reproductionType=reproductionType)
  
  #now over-write with the desired coefficients!
  for (j in 1:length(Formula)) { 
    fv1@fitFec[[j]]$coefficients <- coeff[[j]]
  }
  
  return(fv1)
  
}

# =============================================================================
# =============================================================================

Try the IPMpack package in your browser

Any scripts or data that you put into this service are public.

IPMpack documentation built on May 29, 2017, 3:45 p.m.