Nothing
# 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(fecObj@offspringRel)))
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(fecObj@offspringRel)))
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(fecObj@offspringRel)>1) {
# if (length(grep("logsize", fecObj@offspringRel$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(fecObj@fecNames,names(fecObj@fecConstants))
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(discreteTrans@discreteTrans)[k],colnames(discreteTrans@discreteTrans)[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(discreteTrans@discreteTrans))
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)
}
# =============================================================================
# =============================================================================
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.