R/IPMpack-Util.r
In wpetry/IPMpack2: Builds and analyses Integral Projection Models (IPMs).
Defines functions
invLogit
sampleIPMOutput
sampleIPM
sampleVitalRateObj
coerceSurvObj
coerceGrowthObj
addPdfGrowthPic
plotSurvModelComp
plotGrowthModelComp
survModelComp
growthModelComp
convertIncrement
simulateCarlina
sampleSequentialIPMs
generateData
picGrow
picSurv
Documented in
addPdfGrowthPic
coerceGrowthObj
coerceSurvObj
convertIncrement
generateData
growthModelComp
invLogit
picGrow
picSurv
plotGrowthModelComp
plotSurvModelComp
sampleIPM
sampleIPMOutput
sampleSequentialIPMs
sampleVitalRateObj
simulateCarlina
survModelComp
# =============================================================================
# =============================================================================
## FUNCTION FOR MAKING PICTURE OF SURVIVAL DATA AND FITTED SURVIVAL OBJECT ####
## and PICTURE GROWTH DATA AND FITTED GROWTH OBJECT
# ** note that growth may need redefining if new growth objects are defined
# parameters - dataf - the dataframe
# - survObj - a survival object
# - ncuts - the number of cuts used for binning survival
# returns -
picSurv <- function(dataf, survObj, ncuts = 20, makeTitle = "Survival", ...) {
#organize data and plot mean of ncut successive sizes, so trend more obvious
os<-order(dataf$size); os.surv<-(dataf$surv)[os]; os.size<-(dataf$size)[os];
psz<-tapply(os.size,as.numeric(cut(os.size,ncuts)),mean,na.rm=TRUE); #print(psz)
ps<-tapply(os.surv,as.numeric(cut(os.size,ncuts)), mean, na.rm = TRUE);#print(ps)
if (length(grep("covariate",names(survObj@fit$model)))==0) {
#plot data
plot(as.numeric(psz),as.numeric(ps),pch=19,
xlab="Size at t", ylab = "Survival to t+1", main = makeTitle, ...)
#Plot fitted models
points(dataf$size[order(dataf$size)],surv(dataf$size[order(dataf$size)],data.frame(covariate=1),survObj),type="l",col=2)
} else {
plot(as.numeric(psz),as.numeric(ps),
type="n",pch=19,xlab="Size at t", ylab="Survival to t+1",main="Survival",...)
dataf$covariate <- as.factor(dataf$covariate)
levels(dataf$covariate) <- 1:length(unique(dataf$covariate))
os.cov<-(dataf$covariate)[os]
sizes <- dataf$size[!is.na(dataf$size)]; sizes <- sizes[order(sizes)]
ud <- unique(dataf$covariate); ud <- ud[!is.na(ud)]
for (k in 1:length(ud)) {
tp <- os.cov==ud[k]
psz<-tapply(os.size[tp], as.numeric(cut(os.size[tp], ncuts)), mean, na.rm = TRUE); #print(psz)
ps<-tapply(os.surv[tp],as.numeric(cut(os.size[tp],ncuts)),mean,na.rm=TRUE);#print(ps)
points(as.numeric(psz),as.numeric(ps),pch=19,col=k)
newd <- data.frame(size=sizes,size2=sizes^2,size3=sizes^3,
covariate=rep(as.factor(ud[k]),length(sizes)))
if(length(grep("expsize",survObj@fit$formula))==1)
newd$expsize=exp(sizes)
if(length(grep("logsize",survObj@fit$formula))==1)
newd$logsize=log(sizes)
if(length(grep("logsize2",survObj@fit$formula))==1)
newd$logsize=(log(sizes))^2
pred.surv <- predict(survObj@fit,newd,type="response")
points(newd$size,pred.surv,type="l",col=k)
}
}
}
# =============================================================================
# =============================================================================
# Function to make pic of growth object fit with data
#
# parameters - dataf - the dataframe
# - growObj - a growth object
#
# returns -
#
picGrow <- function(dataf, growObj, mainTitle = "Growth",...) {
predVar <- attr(growObj@fit$terms,"predvars")[[2]] #jess quick fix. this function does not work with changingVar either at the moment
if (class(growObj)=="growthObjTruncIncr") {
predVar <- "incr"
} else {
predVar <- attr(growObj@fit$terms,"predvars")[[2]]
}
if(predVar == "sizeNext") {
plot(dataf$size, dataf$sizeNext,pch=19, xlab="Size at t", ylab="Size at t+1", main = mainTitle,...)
abline(a = 0, b = 1)
}else{
dataf$incr <- dataf$sizeNext - dataf$size
plot(dataf$size, dataf$incr, pch = 19, xlab = "Size at t", ylab="Size increment", main = mainTitle,...)
abline(a = 0, b = 0)
}
if (length(grep("covariate", names(growObj@fit$model))) > 0) {
#convert to 1:n for indexing later and to relate to discrete
dataf$covariate <- as.factor(dataf$covariate)
levels(dataf$covariate) <- 1:length(unique(dataf$covariate))
ud <- unique(dataf$covariate); ud <- ud[!is.na(ud)]
for (k in 1:length(ud)) {
tp <- dataf$covariate == ud[k]
points(dataf$size[tp], dataf$sizeNext[tp], pch = 19, col = k)
}
ud <- as.factor(ud)
} else {
ud <- 0
}
sizes <- dataf$size[!is.na(dataf$size)]; sizes <- sizes[order(sizes)]
for (k in 1:length(ud)) {
newd <- data.frame(size = sizes, size2 = sizes ^ 2,size3 = sizes ^ 3,
covariate = as.factor(rep(ud[k],length(sizes))))
if(length(grep("expsize", names(growObj@fit$coefficients))) == 1)
newd$expsize = exp(sizes)
if(length(grep("logsize", names(growObj@fit$coefficients))) == 1)
newd$logsize = log(sizes)
if(length(grep("logsize2", names(growObj@fit$coefficients))) == 1)
newd$logsize=(log(sizes))^2
if (length(grep("decline",tolower(as.character(class(growObj)))))>0 |
length(grep("trunc",tolower(as.character(class(growObj)))))>0) {
pred.size <- .predictMuX(growObj, newd, covPred = k)
} else {
pred.size <- predict(growObj@fit,newd,type = "response")
}
if (length(grep("incr", tolower(as.character(class(growObj))))) == 0) {
points(sizes, pred.size, type = "l", col = k + 1)
} else {
if (length(grep("logincr", tolower(as.character(class(growObj))))) > 0) {
points(sizes, sizes + exp(pred.size), type = "l", col = k + 1)
} else {
lines(sizes, pred.size, col = k + 1)
}
}
}
}
# =============================================================================
# =============================================================================
## FUNCTIONS FOR SIMULATING DATA ############################
## Return a data-frame with the headings used in the 'makeObject' functions
## Generate a simple data-frame for only continuous covariates
# with a total of 1000 measurements with columns called
# size, sizeNext, surv, covariate, covariateNext, fec,
#
#
generateData <- function(nSamp=1000, type="simple"){
if (nSamp<100) stop("Error: nSamp in generateData should be at least 100")
if (type=="simple"){
covariate <- sample(0:1, size=nSamp, replace=TRUE, prob = c(0.2, 0.8))
covariateNext <- sample(0:1, size=nSamp, replace=TRUE, prob = c(0.8, 0.2))
size <- rnorm(nSamp,5,2)
#size <- exp(rnorm(1000, -1, 1.1))
sizeNext <- 1+0.8*size-0.9*covariate+rnorm(nSamp,0,1)
seedlings <- sample(1:nSamp,size=100,replace=TRUE)
size[seedlings] <- NA; sizeNext[seedlings] <- rnorm(100,2,0.5)
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 <- fec*10
stage <- stageNext <- rep("continuous",nSamp)
stage[is.na(size)] <- NA
stageNext[is.na(sizeNext)] <- "dead"
dataf <- data.frame(size=size,sizeNext=sizeNext,surv=surv,
covariate=covariate,covariateNext=covariateNext,
fec=fec, stage=stage,stageNext=stageNext)
dataf$sizeNext[dataf$surv==0] <- NA
}
if (type=="discrete") dataf <-.generateDataDiscrete(nSamp=nSamp)
if (type=="stochastic") dataf <-.generateDataStoch(nSamp=nSamp)
if (type!="simple" & type!="stochastic" & type!="discrete"){
stop("unsupported type: supported types include simple, stochastic, or discrete")
}
return(dataf)
}
# =============================================================================
# =============================================================================
## FUNCTION FOR MAKING A LIST OF IPMS ############################################
# to do for stoch env with a single discrete covariate. ##########################
sampleSequentialIPMs <- function(dataf, nBigMatrix=10, minSize=-2,maxSize=10,
integrateType="midpoint", correction="none",
explSurv=surv~size+size2+covariate,
explGrow=sizeNext~size+size2+covariate,
regType="constantVar",explFec=fec~size,Family="gaussian",
Transform="none",fecConstants=data.frame(NA)) {
#convert to 1:n for indexing later
dataf$covariate <- as.factor(dataf$covariate)
levels(dataf$covariate) <- 1:length(unique(dataf$covariate))
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$covariate)
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(covariate=as.factor(k)),
fecObj = fv1,integrateType=integrateType, correction=correction)
tpS <- makeIPMPmatrix(nBigMatrix = nBigMatrix, minSize = minSize,
maxSize = maxSize, chosenCov = data.frame(covariate=as.factor(k)),
growObj = gr1, survObj = sv1,
integrateType=integrateType, correction=correction)
IPM.list[[k]] <- tpF+tpS
}
return(IPM.list)
}
# =============================================================================
# =============================================================================
### simulation Carlina ######################################################
## Function to simulate something a bit like Carlina
##
## Parameters - nSamp - starting pop sizes
# - nYrs - no years in simulation
# - nSampleYrs - no years to extract for the data
# - ... bunch of parameters
# - meanYear - means for year effects
# - matVarYear - variance covariances for year effects
# - densDep - model density dependence in seedling establishment or not
#
# Returns - list including dataf - data-frame of data
# - various of the simulation parameters for convenience
#
# PARAMETERS FROM TABLE 2 IN REES AND ELLNER 2009
simulateCarlina <- function(nSamp=200,nYrs=1000,nSampleYrs=15,
m0=-1.37,ms=0.59,
b0=-12.05,bs=3.64,
A=-1,B=2,
ag=1.14,bg=0.74,sig=0.29,
mean.kids=3.16,sd.kids=0.5,
meanYear=c(0,0,0),
matVarYear=matrix(c(1.03,0,0,0,0.037,0.041,0,0.041,0.075),3,3),
varA=0,varB=0,
densDep=TRUE,maxPerYr=1000,maxStoreSeedlingsPerYr=200,sizes=c()) {
require(mvtnorm)
#initiate and set up year index
if (length(sizes)==0) sizes <- rnorm(nSamp,3,0.5)
startYr <- (nYrs-nSampleYrs)
#recruits
nrec <- c(20,42,12,17,8,19,58,45,44,2,56,25,75,92,94,6,4,34,104)
#total plants
totpl <- c(21,57,47,25,25,33,88,94,97,26,85,80,122,175,160,10,6,189)
#set up dataframe
dataf <-matrix(NA,(maxPerYr+maxStoreSeedlingsPerYr)*nYrs,13)
maxPop <- nrow(dataf)
n.per.yr <- rep(NA,nYrs)
trueGrow <- rep(NA,nYrs)
count <- 0
#stoch sims
if (sum(matVarYear)!=0) {
tmp <- rmvnorm(nYrs,mean=meanYear,sigma=matVarYear)
} else {
tmp <- matrix(0,nYrs,3)
}
if (varA!=0) tmpA <- rnorm(nYrs,mean=A,sd=sqrt(varA))
if (varB!=0) tmpB <- rnorm(nYrs,mean=B,sd=sqrt(varB))
for (t in 1:nYrs) {
#yr effects
nSeedlings <- sample(nrec,size=1,replace=FALSE)
#print(tmp)
n.per.yr[t] <- length(sizes)
m.year <- tmp[t,1]
cg.year <- tmp[t,2]
b.year <- tmp[t,3]
if (varA!=0) A.year <- tmpA[t] else A.year <- A
if (varB!=0) B.year <- tmpB[t] else B.year <- B
#print(c(m.year,cg.year,b.year))
if (n.per.yr[t]>0) {
#survival
sx <- 1*(invLogit(m0+ms*sizes+m.year)>runif(n.per.yr[t]))
#flowering
fx <- 1*(invLogit(b0+bs*sizes)>runif(n.per.yr[t]))
#fertility
seedsx <- exp(A.year+B.year*sizes)*fx*sx
#seedsx[seedsx>0] <- rpois(sum(seedsx>0),seedsx[seedsx>0])
} else {
sx <- fx <- seedsx <- c()
print(c("extinct in year ", t))
}
if (densDep) pEst <- min(nSeedlings/max(sum(seedsx, na.rm=TRUE),1),1) else pEst <- 1
babies <- rnorm(ceiling(pEst*max(sum(seedsx,na.rm=TRUE),1)),
mean.kids+b.year,sd.kids)
#if (count>0) print(c("yr",t-startYr,pEst*max(sum(seedsx,na.rm=TRUE),1),length(babies)))
#will end up with nrec babies at least in density dependent case
if (length(babies)<nSeedlings & densDep) nSeedlings <- length(babies)
#growth
sizeNext <- rnorm(length(sizes),ag+bg*sizes+cg.year,sig)
#remove dead or flowered
fx[sx==0] <- NA
sizeNext[sx==0 | fx==1] <- NA
#storage
if (t > startYr) {
#print(t)
#don't 'do it at all if too big (so as not to do only adults, etc)
if ((count+length(sizes))>maxPop) { print("large pop size, breaking");break()}
if ((count+length(sizes)+length(babies))>maxPop) { print("large pop size, breaking");break()}
chs <- (count+1):(count+length(sizes))
dataf[chs,1] <- sizes
dataf[chs,2] <- sizeNext
dataf[chs,3] <- sx
dataf[chs,4] <- fx
dataf[chs,5] <- seedsx
dataf[chs,6] <- t
dataf[chs,7] <- nSeedlings
dataf[chs,8] <- m.year
dataf[chs,9] <- cg.year
dataf[chs,10] <- b.year
dataf[chs,12] <- A.year
dataf[chs,13] <- B.year
count <- count + length(sizes)
nbabes.store <- min(length(babies),maxStoreSeedlingsPerYr)
chs <- (count+1):(count+nbabes.store)
dataf[chs,2] <- babies[1:nbabes.store]
dataf[chs,6] <- t
dataf[chs,7] <- nSeedlings
dataf[chs,8] <- m.year
dataf[chs,9] <- cg.year
dataf[chs,10] <- b.year
dataf[chs,11] <- "sexual"
dataf[chs,12] <- A.year
dataf[chs,13] <- B.year
count <- count + nbabes.store
#print(nbabes.store)
}
#new pop - i.e. those that grew, survive, and didn't flower; as well as babies
sizes <- c(sizeNext[sx==1 & fx==0 & !is.na(fx)],babies)
if (length(sizes)==0) print("extinct")
#get true growth rate
trueGrow[t] <- log(length(sizes)/n.per.yr[t])
#cull size at t t down to maxPerYr, if not density dependent
if (!densDep) sizes <- sample(sizes,size=maxPerYr, replace=TRUE)
}
list.par <- list(m0=m0,ms=ms,
b0=b0,bs=bs,
A=A,B=B,varB=varB,
ag=ag,bg=bg,sig=sig,
mean.kids=mean.kids,sd.kids=sd.kids,
meanYear=meanYear,
matVarYear=matVarYear,
nrec=nrec)
dataf <- dataf[1:count,]
dataf <- data.frame(dataf,stringsAsFactors = FALSE)
#print(table(dataf[,6]))
#plot(cumsum(trueGrow[startYr:nYrs])/(1:length(trueGrow[startYr:nYrs])),type="l")
vartrueGrow <- var(trueGrow[startYr:nYrs],na.rm=TRUE)
meantrueGrow <- mean(exp(trueGrow[startYr:nYrs]),na.rm=TRUE)
trueGrow <- mean(trueGrow[startYr:nYrs],na.rm=TRUE)
#abline(h=trueGrow)
colnames(dataf) <- c("size","sizeNext","surv","flower","fec",
"year","nSeedlings","m.year","cg.year","b.year","offspringNext",
"A.year","B.year")
dataf$size <- as.numeric(dataf$size)
dataf$sizeNext <- as.numeric(dataf$sizeNext)
dataf$surv <- as.numeric(dataf$surv)
dataf$flower <- as.numeric(dataf$flower)
dataf$fec <- as.numeric(dataf$fec)
dataf$nSeedlings <- as.numeric(dataf$nSeedlings)
dataf$year <- as.numeric(dataf$year)
dataf$cg.year <- as.numeric(dataf$cg.year)
dataf$m.year <- as.numeric(dataf$m.year)
dataf$b.year <- as.numeric(dataf$b.year)
dataf$A.year <- as.numeric(dataf$A.year)
dataf$B.year <- as.numeric(dataf$B.year)
#print(table(dataf$year))
dataf$fec[dataf$fec==0] <- NA
return(list(dataf=dataf,meanYear=meanYear,matVarYear=matVarYear,
list.par=list.par,trueGrow=trueGrow,
vartrueGrow=vartrueGrow,meantrueGrow=meantrueGrow))
}
# =============================================================================
# =============================================================================
## Convert Increment - where exact dates of census vary but some multiplier of yearly increments
## are desired; this function takes a data-frame (with columns size, sizeNext,
## and, importantly exactDate, exactDatel)
## and returns a data-frame with sizeNext modified proportional to the time
## elapsed in the desired yaerly increments, adding an additional column denoted 'increment'.
#
# Parameters - dataf - a dataframe with headings size, sizeNext, exactDate,exactDatel
# - nYrs - the number of years between censuses desired (e.g. for CTFS data, 5 years)
#
# Returns - a dataframe with the same headings
convertIncrement <- function(dataf, nYrs=1) {
incr <- dataf$sizeNext - dataf$size
if (class(dataf$exactDatel) == "Date") {
timeElapsed <- (difftime(dataf$exactDatel,dataf$exactDate)[[1]])/(365*nYrs)
} else {
timeElapsed <- (dataf$exactDatel-dataf$exactDate)/(365*nYrs)
}
incrNew <- incr / timeElapsed
dataf$sizeNext <- dataf$size + incrNew
dataf$incr <- dataf$incrNew
return(dataf)
}
# =============================================================================
# =============================================================================
# Function to compare model fits for growth and survival objects built with different linear combinations of covariates.
# Growth can have multiple response forms. Uses .makeCovDf. Separate plot functions.
#
#
# Returns - list with a summary table of covariates and scores, and a sub-list of growth or survival objects.
#
growthModelComp <- function(dataf,
expVars = c(sizeNext~1, sizeNext~size, sizeNext~size + size2),
regressionType = "constantVar",
testType = "AIC",
makePlot = FALSE,
mainTitle = "",
plotLegend = TRUE,
legendPos = "topright", ...) {
varN <- length(expVars)
typeN <- length(regressionType)
treatN <- varN * typeN
summaryTable <- data.frame()
grObj <- vector("list", length = treatN)
i <- 1
for(v in 1:varN) {
for(t in 1:typeN) {
grObj[[i]] <- makeGrowthObj(dataf = dataf, regType = regressionType[t], Formula = expVars[[v]])
if (length(grep("decline",tolower(as.character(class(grObj[[i]])))))>0) {
summaryTable <- rbind(summaryTable, cbind(as.character(unlist(expVars))[v], regressionType[t], match.fun(testType)(grObj[[i]]@fit$fit)))
} else {
summaryTable <- rbind(summaryTable, cbind(as.character(unlist(expVars))[v], regressionType[t], match.fun(testType)(grObj[[i]]@fit)))
}
i <- i + 1
}
}
summaryTable <- as.data.frame(summaryTable)
names(summaryTable) <- c("Exp. Vars", "Reg. Type", testType)
outputList <- list(summaryTable = summaryTable, growthObjects = grObj)
# PLOT SECTION #
if(makePlot == TRUE) {
plotGrowthModelComp(grObj = grObj, summaryTable = summaryTable, dataf = dataf, expVars = expVars, testType = "AIC", plotLegend = plotLegend, mainTitle = mainTitle, legendPos, ...)
}
return(outputList)
}
# =============================================================================
# =============================================================================
survModelComp <- function(dataf,
expVars = c(surv~1, surv~size, surv~size + size2),
testType = "AIC",
makePlot = FALSE,
mainTitle = "", ncuts = 20,
plotLegend = TRUE,
legendPos = "bottomleft", ...) {
varN <- length(expVars)
treatN <- varN
summaryTable <- data.frame()
svObj <- vector("list", length = treatN)
i <- 1
for(v in 1:varN) {
svObj[[i]] <- makeSurvObj(dataf = dataf, Formula = expVars[[v]])
summaryTable <- rbind(summaryTable, cbind(as.character(unlist(expVars))[v], match.fun(testType)(svObj[[i]]@fit)))
i <- i + 1
}
summaryTable <- as.data.frame(summaryTable)
names(summaryTable) <- c("Exp. Vars", testType)
outputList <- list(summaryTable = summaryTable, survObjects = svObj)
# PLOT SECTION #
if(makePlot == TRUE) {
## this is the surv picture
plotSurvModelComp(svObj = svObj, summaryTable = summaryTable, dataf = dataf, expVars = expVars, testType = "AIC", plotLegend = plotLegend, mainTitle = mainTitle, ncuts = ncuts, legendPos = legendPos, ...)
}
return(outputList)
}
# =============================================================================
# =============================================================================
# Plot functions for model comparison. Plots the series of fitted models for growth and survival objects.
# Can plot a legend with the model covariates and model test criterion scores (defaults to AIC).
plotGrowthModelComp <- function(grObj, summaryTable, dataf, expVars, testType = "AIC", plotLegend = TRUE, mainTitle = "", legendPos = "topright", ...) {
treatN <- length(grObj)
sizeSorted <- unique(sort(dataf$size))
if(length(grep("sizeNext", unlist(as.character(expVars))))) {
y.lab <- "Size at t + 1"
dataSizeNext <- dataf$sizeNext
}
if(length(grep("incr", unlist(as.character(expVars))))) {
y.lab <- "Growth"
dataSizeNext <- dataf$sizeNext - dataf$size
}
if(length(grep("logincr", unlist(as.character(expVars))))){
y.lab <- "log(growth)"
dataSizeNext <- log(dataf$sizeNext - dataf$size)
}
plot(dataf$size, dataSizeNext, pch = 19, xlab = "Size at t", ylab = y.lab, main = mainTitle, cex = 0.8,...)
for(p in 1:treatN) {
#PROBLEM HERE
expVarsHere <- paste(attr(terms((expVars[[p]])),"term.labels"),collapse="+")
newd <- .makeCovDf(sizeSorted, expVarsHere)
if (length(grep("decline",tolower(as.character(class(grObj[[p]])))))>0) {
pred.size <- .predictMuX(grObj[[p]], newd)
} else { pred.size <- predict(grObj[[p]]@fit, newd, type = "response")}
lines(sizeSorted, pred.size, type = "l", col = (p + 1))
}
if(plotLegend) {
legend(legendPos, legend = sprintf("%s: %s = %.1f", expVars, testType, as.numeric(as.character(summaryTable[,3]))), col = c(2:(p + 1)), lty = 1, xjust = 1, bg = "white")
}
}
# =============================================================================
# =============================================================================
plotSurvModelComp <- function(svObj, summaryTable, dataf, expVars, testType = "AIC", plotLegend = TRUE, mainTitle = "", ncuts = 20, legendPos = "bottomleft", ...) {
treatN <- length(svObj)
#ncuts <- 20 # survival bins
os <- order(dataf$size) # order size
osSurv <- (dataf$surv)[os] # order survival data according to size
osSize<-(dataf$size)[os] # ordered size data
binnedSize <- tapply(osSize, as.numeric(cut(osSize, breaks=ncuts)), mean, na.rm = TRUE); # bin Size data
binnedSurv <- tapply(osSurv, as.numeric(cut(osSize, breaks=ncuts)), mean, na.rm = TRUE) #bin Survival probabilities
plot(binnedSize, binnedSurv, pch = 19, xlab = "Size at t", ylab = "Survival to t + 1", main = mainTitle, cex = 0.8,...)
for(p in 1:treatN) {
expVarsHere <- paste(attr(terms(expVars[[p]]),"term.labels"),collapse="+")
newd <- .makeCovDf(osSize, expVarsHere[p])
lines(dataf$size[order(dataf$size)], surv(dataf$size[os], data.frame(covariate=1), svObj[[p]]), col = (p + 1))
}
if(plotLegend) {
legend(legendPos, legend = sprintf("%s: %s = %.1f", expVars, testType, as.numeric(as.character(summaryTable[,2]))), col = c(2:(p + 1)), lty = 1, xjust = 1, bg = "white")
}
}
# =============================================================================
# =============================================================================
## Function to add pdf to model comparison pics
addPdfGrowthPic <- function(respType = "sizeNext",
sizesPlotAt=c(20,50,60),
sizeRange=c(20,400),
incrRange=c(-10,50),
scalar=100,
growthObjList,
cols=1:5,
cov=data.frame(covariate=1),
minShow=1e-2,
jitt=2, #how far apart should pdfs be if you are plotting several
...){
nval <- length(growthObjList)
sizes <- seq(sizeRange[1],sizeRange[2],length=500)
incr <- seq(incrRange[1],incrRange[2],length=500)
logincr <- log(incr); logincr[!is.finite(logincr)] <- NA
for (j in 1:nval) {
#print(j)
for (k in 1:length(sizesPlotAt)) {
if (respType=="sizeNext") {
pred <- growth(sizesPlotAt[k],sizes,cov,growthObjList[[j]])*scalar
pred[pred<minShow] <- NA
points(sizesPlotAt[k]+pred+jitt*(j-1),
sizes,type="l", col=cols[j],...)
}
if (respType=="incr") {
pred <- growth(sizesPlotAt[k],sizesPlotAt[k]+incr,cov,growthObjList[[j]])*scalar
pred[pred<minShow] <- NA
points(sizesPlotAt[k]+pred+jitt*(j-1),
incr,type="l", col=cols[j],...)
}
if (respType=="logincr") {
pred <- growth(sizesPlotAt[k],sizesPlotAt[k]+logincr,cov,growthObjList[[j]])*scalar
pred[pred<minShow] <- NA
#print(range(pred,na.rm=TRUE))
points(sizesPlotAt[k]+pred+jitt*(j-1),
logincr,type="l", col=cols[j],...)
}
}}
}
# =============================================================================
# =============================================================================
## Function to coerce Growth object to parameters and variance desired
coerceGrowthObj <- function(growthObj, coeff, sd){
if (length(growthObj@fit$coefficients) !=length(coeff)) {
print("warning: number of desired coefficients to not match number of growth object coefficients")
}
growthObj@fit$coefficients[] <- as.numeric(coeff)
growthObj@sd[] <- as.numeric(sd)
return(growthObj)
}
# =============================================================================
# =============================================================================
## Function to coerce Survival object to parameters desired
coerceSurvObj <- function(survObj, coeff){
if (length(survObj@fit$coefficients) !=length(coeff)) print("warning: number of desired coefficients to not match number of growth object coefficients")
survObj@fit$coefficients[] <- as.numeric(coeff)
return(survObj)
}
# =============================================================================
# =============================================================================
## Parametric boostrap of vital rate objects
# created by cory on 6-4-13 as a modification of getListRegObjs()
sampleVitalRateObj <- function(Obj,nSamp=100,nDiscreteGrowthTransitions=NULL,nDiscreteOffspringTransitions=NULL,nOffspring=NULL) {
require(mvtnorm)
require(MCMCpack)
objList <- list()
#generate new set parameters from mvn
# what is the difference between offspring splitter transitions and discrete
# for discrete transition matrices objects
if( sum(class(Obj)==c("discreteTrans","discreteTransInteger")) ) {
if(is.null(nDiscreteGrowthTransitions)){
stop('Error: Please specify the total number of transitions that were used to estimate the discrete transitions in Obj@discreteTrans to allow the function to estimate the appropriate variance for resampling them')
}
if(!is.null(nDiscreteGrowthTransitions)){
for (j in 1:nSamp) {
objList[[j]] <- Obj
for (k in 1:ncol(Obj@discreteTrans)){
objList[[j]]@discreteTrans[,k] <- rdirichlet(1, nDiscreteGrowthTransitions * as.numeric(Obj@discreteTrans[,k]))
}
}
}
}
# for growth and survival objects
if( !sum(class(Obj)==c("fecObj","fecObjInteger","discreteTrans", "discreteTransInteger")) ) {
newpar <- rmvnorm(nSamp, mean = Obj@fit$coefficients, sigma = vcov(Obj@fit))
for (j in 1:nSamp) {
objList[[j]] <- Obj
objList[[j]]@fit$coefficients <- newpar[j,]
}
}
# for fecundity objects
if( sum(class(Obj)==c("fecObj","fecObjInteger")) ) {
for (j in 1:nSamp) {
# sample fecundity regression parameters
for (k in 1:length(Obj@fitFec)) {
newpar <- rmvnorm(1, mean = Obj@fitFec[[k]]$coefficients,
sigma = vcov(Obj@fitFec[[k]]))
objList[[j]] <- Obj
objList[[j]]@fitFec[[k]]$coefficients <- newpar
}
# sample discrete transitions from dirichlet
if(ncol(Obj@offspringSplitter)>1){ # only if there are discrete fec stages
if(is.null(nDiscreteOffspringTransitions)){
stop('Error: Please specify the total number of discrete transitions that were used to estimate the fecundity transition in Obj@offspringSplitter to allow the function to estimate the appropriate variance for resampling them')
}
if(!is.null(nDiscreteOffspringTransitions)){
objList[[j]]@offspringSplitter[] <- rdirichlet(1, nDiscreteOffspringTransitions * as.numeric(Obj@offspringSplitter))
}
}
# sample the offspring sizes
# TODO: OffspringRel currently does not store all the information from the regression, which it should so that this works with offspring models that have more than an intercept
if(is.null(nOffspring)){
stop('Error: Please specify the total number of offspring that were used to estimate the offspring size distribtuion in Obj@offspringRel to allow the function to estimate the appropriate variance for resampling them')
}
# TODO: Make work with the new OffspringObjs
if(!is.null(nOffspring)){
newpar3 <- rnorm(1,coefficients(Obj@offspringRel)[1], Obj@sdOffspringSize/sqrt(nOffspring))
objList[[j]]@offspringRel$coefficients[] <- newpar3
newpar4 <- rnorm(1,Obj@sdOffspringSize, sqrt(2*Obj@sdOffspringSize^4)/(nOffspring-1))
objList[[j]]@sdOffspringSize <- newpar4
}
}
}
return(objList)
}
# =============================================================================
# =============================================================================
## Use list of vital rate objects to make List of IPMS
# created by cory on 6-4-13 as a modification of getListIPMs()
# TODO: MAKE WORK FOR OFFSPRING OBJS
sampleIPM<- function(
growObjList=NULL,survObjList=NULL,fecObjList=NULL,
offspringObjList=NULL, discreteTransList=1,
nBigMatrix,minSize,maxSize,
covariates=FALSE,envMat=NULL,
integrateType="midpoint",correction="none",warn=TRUE) {
if(!is.null(offspringObjList)){
stop('Sorry, lists of offspring objects are not yet supported. However, specify offspring size distributions via makeFecObj() is supported. Please set offspringObjList=NULL, or modify the function to accept offspring objects and email it to us.')
}
# make the same number of samples for each vital rate object
n.samples <- max(c(length(growObjList),length(survObjList),length(fecObjList)))
if(!is.null(growObjList) & length(growObjList)<n.samples){
if(warn) warning('Length of growth object list is less than the length of another vital rate object list, so some members of the growth object list have been repeated.')
growthObjList <- sample(growObjList,size=n.samples,replace=T)
}
if(!is.null(survObjList) & length(survObjList)<n.samples ){
if(warn) warning('Length of survival object list is less than the length of another vital rate object list, so some members of the survival object list have been repeated.')
survObjList <- sample(survObjList,size=n.samples,replace=T)
}
if(!is.null(fecObjList) & length(fecObjList)<n.samples ){
if(warn) warning('Length of fec object list is less than the length of another vital rate object list, so some members of the fec object list have been repeated.')
fecObjList <- sample(fecObjList,size=n.samples,replace=T)
}
if(!is.null(growObjList) & length(discreteTransList)<n.samples ){
# 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=n.samples,replace=T)
}
if(!is.null(growObjList)){
PmatrixList <- .makeListPmatrix(growObjList=growObjList, survObjList=survObjList, nBigMatrix=nBigMatrix, minSize=minSize, maxSize=maxSize, cov=covariates, envMat=envMat,discreteTransList=discreteTransList, integrateType=integrateType, correction=correction)
matrixList <- PmatrixList
}
if(!is.null(fecObjList)){
FmatrixList <- .makeListFmatrix(fecObjList, nBigMatrix=nBigMatrix, minSize=minSize, maxSize=maxSize, cov=covariates, envMat=envMat, integrateType=integrateType, correction=correction)
matrixList <- FmatrixList
}
if(!is.null(growObjList) & !is.null(fecObjList)){
for(i in 1:n.samples){
matrixList[[i]] <- PmatrixList[[i]] # to get all the Pmatrix slots
matrixList[[i]]@.Data <- PmatrixList[[i]]+FmatrixList[[i]]
}
}
return(matrixList)
}
# =============================================================================
# =============================================================================
## Use list of IPMs to make List of IPM output
# created by cory on 6-4-13 as a modification of getIPMOutput() and getIPMOutputDirect()
sampleIPMOutput <- function(IPMList=NULL,PMatrixList=NULL,passageTimeTargetSize=c(), sizeToAgeStartSize=c(),sizeToAgeTargetSize=c(),warn=TRUE){
# removed from old version of getIPMOutputDirect
# # cov=FALSE, envMat=NULL,
# # chosenCov = data.frame(covariate = 1),
# # onlyLowerTriGrowth=FALSE)
# if only a Pmatrix is supplied, just call it IPMList for simplicity
if(is.null(IPMList)) IPMList=PMatrixList
if ( is.null(passageTimeTargetSize) ) {
if(warn) print("no target size for passage time provided; taking meshpoint median")
passageTimeTargetSize <- median(IPMList[[1]]@meshpoints)
}
if ( is.null(sizeToAgeStartSize) ) {
if(warn) print("no starting size for size to age provided; taking minimum size")
sizeToAgeStartSize <- min(IPMList[[1]]@meshpoints)
}
if ( is.null(sizeToAgeTargetSize) ) {
if(warn) print("no target size for size to age provided; taking meshpoint values")
sizeToAgeTargetSize <- IPMList[[1]]@meshpoints
}
nSamps <- length(IPMList)
stableStage <- LE <- pTime <- matrix(NA,nSamps,length(IPMList[[1]]@.Data[1,]))
lambda <- rep(NA,nSamps)
resAge <- resSize <- matrix(NA,nSamps,length(sizeToAgeTargetSize))
h1 <- diff(IPMList[[1]]@meshpoints)[1]
for (k in 1:nSamps) {
IPM <- IPMList[[k]]
LE[k,]<-meanLifeExpect(IPM)
pTime[k,]<-passageTime(passageTimeTargetSize,IPM)
if (is.null(PMatrixList)) {
lambda[k] <- Re(eigen(IPM)$value[1])
stableStage[k,] <- eigen(IPM)$vector[,1]
#normalize stable size distribution
stableStage[k,] <- stableStage[k,]/sum(stableStage[k,])
}
# get size to age results
res2 <- sizeToAge(Pmatrix=IPM,startingSize=sizeToAgeStartSize,
targetSize=sizeToAgeTargetSize)
resAge[k,] <- res2$timeInYears
resSize[k,] <- res2$targetSize
}
return(list(LE=LE,pTime=pTime,lambda=lambda,stableStage=stableStage,
resAge=resAge,resSize=resSize,meshpoints=IPM@meshpoints))
}
# =============================================================================
# =============================================================================
invLogit <- function(x) {
u <- exp(pmin(x, 50))
return(u / (1 + u))
}
# =============================================================================
# =============================================================================
# For a single Pmatrix (!not compound and no discrete stages!), this functions makes a series
# of diagnostic plots - this is defined for growthObj,
# growthObjIncr - modification required
# if other objects used
#
# Parameters - the Pmatrix
# - growObj - growth object used to build it
# - survObj - survival object used to build it
# - dff - the data from which it was built
# - integrateType - "midpoint", or "cumul" - should
# correspond to what the IPM was built with
# - do you want to implement the corrections for integration?
# Returns -
#
diagnosticsPmatrix <- function (Pmatrix, growObj, survObj, dff=NULL,
integrateType = "midpoint",
correction = "none", cov = data.frame(covariate = 1),
sizesToPlot=c(),extendSizeRange=c()) {
print("Range of Pmatrix is ")
print(range(c(Pmatrix)))
if (Pmatrix@meshpoints[1] > 0)
new.min <- Pmatrix@meshpoints[1]/2
else new.min <- Pmatrix@meshpoints[1] * 1.5
new.max <- 1.5 * max(Pmatrix@meshpoints)
if (length(extendSizeRange)>0 & length(extendSizeRange)!=2) print("require two values for extendSizeRange, reflecting upper and lower limits")
if (length(extendSizeRange)>0) { new.min <- extendSizeRange[1]; new.max <- extendSizeRange[2]}
#colours for 1) current; 2) bigger size range; 3) bigger no bins
cols <- c("black","tomato","darkblue")
ltys <- c(1,1,3)
#matrix with bigger size range
Pmatrix1 <- makeIPMPmatrix(nEnvClass = 1, nBigMatrix = length(Pmatrix@meshpoints),
minSize = new.min, maxSize = new.max,
chosenCov = cov, growObj = growObj, survObj = survObj,
integrateType = integrateType, correction = correction)
if (sum(is.na(Pmatrix1)) > 0) {
print("Pmatrix with extended size range returns NAs; changing these to 0, and putting survival value onto diagonal for columns that sum to zero")
Pmatrix1[is.na(Pmatrix1)] <- 0
bad <- which(colSums(Pmatrix1) == 0, arr.ind = TRUE)
if (length(bad) > 0)
Pmatrix1[cbind(bad, bad)] <- surv(size = Pmatrix1@meshpoints[bad],
cov = cov, survObj = survObj)
}
#matrix with bigger number of bins
Pmatrix2 <- makeIPMPmatrix(nEnvClass = 1, nBigMatrix = floor(length(Pmatrix@meshpoints) *
1.5), minSize =Pmatrix@meshpoints[1], maxSize = max(Pmatrix@meshpoints),
chosenCov = cov, growObj = growObj, survObj = survObj,
integrateType = integrateType, correction = correction)
if (sum(is.na(Pmatrix2)) > 0) {
print("Pmatrix with extended number of bins returns NAs; changing these to 0, and putting survival value onto diagonal for columns that sum to zero")
Pmatrix2[is.na(Pmatrix2)] <- 0
bad <- which(colSums(Pmatrix2) == 0, arr.ind = TRUE)
if (length(bad) > 0)
Pmatrix2[cbind(bad, bad)] <- surv(size = Pmatrix2@meshpoints[bad],
cov = cov, survObj = survObj)
}
#start plots - put original Pmatrix in black
#par(mfrow = c(3, 3), bty = "l")
#save plot characters
old.par <- par(no.readonly = TRUE)
par(mfrow = c(1, 3), bty = "l",pty="s", mar=c(5,4,4,1))
if (!is.null(dff)) {
xlims <- range(c(Pmatrix@meshpoints, Pmatrix1@meshpoints,
dff$size, dff$sizeNext), na.rm = TRUE)
a1 <- hist(c(dff$size, dff$sizeNext), xlab = "Sizes observed", axes=FALSE,
ylab = "Frequency", main = "", xlim = xlims, col="lightgrey", border="lightgrey",plot=TRUE)
axis(1); axis(2)
} else {
a1 <- list(); a1$counts <- 1:100
xlims <- range(c(Pmatrix@meshpoints, Pmatrix1@meshpoints))
plot(1:100,type="n",xlab = "Sizes", ylab="", ylim=range(a1$counts), xlim=xlims)
}
lcs <- c(0.8,0.6,0.4)
lcs.x <- mean(xlims)
text(lcs.x ,max(a1$counts)*lcs[1],"Current",pos=3,col=cols[1])
arrows(Pmatrix@meshpoints[1], max(a1$counts)*lcs[1],Pmatrix@meshpoints[length(Pmatrix@meshpoints)], max(a1$counts)*lcs[1],col=cols[1], length=0.1, code=3,lty=ltys[1])
text(lcs.x ,max(a1$counts)*lcs[2],"Extended range",pos=3,col=cols[2])
arrows(Pmatrix1@meshpoints[1], max(a1$counts)*lcs[2],Pmatrix1@meshpoints[length(Pmatrix1@meshpoints)], max(a1$counts)*lcs[2],col=cols[2], length=0.1, code=3,lty=ltys[1])
text(lcs.x ,max(a1$counts)*lcs[3],"Increased no of bins",pos=3,col=cols[3])
arrows(Pmatrix2@meshpoints[1], max(a1$counts)*lcs[3],Pmatrix2@meshpoints[length(Pmatrix2@meshpoints)], max(a1$counts)*lcs[3],col=cols[3], length=0.1, code=3,lty=ltys[1])
#legend("topright", legend = "fitted range", col = "black", lty = 2, bty = "n") ##!!! change this
title("Size range")
#survival sums
lims <- range(c(colSums(Pmatrix), surv(Pmatrix@meshpoints, cov, survObj)))
plot(colSums(Pmatrix), surv(Pmatrix@meshpoints, cov, survObj),
type = "n", xlab = "Surviving in Pmatrix", ylab = "Should be surviving",col=cols[1],lty=ltys[1],
xlim=lims,ylim=lims)
abline(0, 1, lwd=2,col="grey")
points(colSums(Pmatrix), surv(Pmatrix@meshpoints, cov, survObj), type = "l", col = cols[1],lty=ltys[1])
points(colSums(Pmatrix1), surv(Pmatrix1@meshpoints, cov, survObj), type = "l", col = cols[2],lty=ltys[2])
points(colSums(Pmatrix2), surv(Pmatrix2@meshpoints, cov, survObj), type = "l", col = cols[3],lty=ltys[3])
title("Survival")
#text(lims[2],lims[2],"(0,1)",pos=1)
# mtext("Points should lie along the 0,1 line, shown in grey", side=1,outer=TRUE,line=-1)
LE <- meanLifeExpect(Pmatrix)
LE1 <- meanLifeExpect(Pmatrix1)
LE2 <- meanLifeExpect(Pmatrix2)
plot(Pmatrix@meshpoints, LE, type = "l", xlim = range(Pmatrix1@meshpoints),
ylim = range(c(LE, LE1)), xlab = "Sizes", ylab = "Life expectancy",col=cols[1],lty=ltys[1])
points(Pmatrix1@meshpoints, LE1, type = "l", col = cols[2],lty=ltys[2])
points(Pmatrix2@meshpoints, LE2, type = "l", col = cols[3],lty=ltys[3])
legend("topleft", legend = c("Current", "Extended range", "Increased no of bins"), col = cols, lty = ltys, bty = "n")
title("Life expectancy")
#mtext("Increasing size range (red) or number of bins (blue) should not alter life expectancy estimates", side=1,outer=TRUE,line=-1)
if (length(sizesToPlot)==0 & !is.null(dff)) sizesToPlot <- as.numeric(quantile(dff$size,c(0.25, 0.5,0.75),na.rm=TRUE))
if (length(sizesToPlot)==0 & is.null(dff)) sizesToPlot <- as.numeric(quantile(Pmatrix@meshpoints,c(0.25, 0.5,0.75),na.rm=TRUE))
loctest <- rep(NA,length(sizesToPlot))
print("Please hit any key for the next plot")
scan(quiet="TRUE")
par(mfrow = c(2, 3), bty = "l",pty="s")
for (kk in c(1,3)) {
if (kk==1) Pmat <- Pmatrix
if (kk==3) Pmat <- Pmatrix2
h <- diff(Pmat@meshpoints)[1]
testSizes <- seq(min(Pmat@meshpoints), max(Pmat@meshpoints),
length = 5000)
for (j in 1:3) {
loctest[j] <- which(abs(Pmat@meshpoints-sizesToPlot[j])==min(abs(Pmat@meshpoints-sizesToPlot[j])),arr.ind=TRUE)[1]
#print(loctest[j])
ps <- surv(Pmat@meshpoints[loctest[j]], cov, survObj)
newd <- data.frame(size = Pmat@meshpoints[loctest[j]],
size2 = Pmat@meshpoints[loctest[j]]^2,
size3 = Pmat@meshpoints[loctest[j]]^3,
covariate = Pmat@env.index[1])
if (length(grep("expsize", names(growObj@fit$coefficients))))
newd$expsize = log(Pmat@meshpoints[loctest[j]])
if (length(grep("logsize", names(growObj@fit$coefficients))))
newd$logsize = log(Pmat@meshpoints[loctest[j]])
if (length(growObj@fit$model$covariate) > 0)
if (is.factor(growObj@fit$model$covariate))
newd$covariate <- as.factor(newd$covariate)
if (length(grep("decline", tolower(as.character(class(growObj))))) > 0 |
length(grep("trunc", tolower(as.character(class(growObj))))) > 0) {
mux <- .predictMuX(growObj, newd)
} else {
mux <- predict(growObj@fit, newd, type = "response")
}
if (length(grep("incr", tolower(as.character(class(growObj))))) > 0 &
length(grep("logincr", tolower(as.character(class(growObj))))) == 0) {
mux <- Pmat@meshpoints[loctest[j]] + mux
}
if (length(grep("decline", tolower(as.character(class(growObj))))) == 0 &
length(grep("trunc", tolower(as.character(class(growObj))))) == 0) {
sigmax2 <- growObj@sd^2
} else {
sigmax2 <- growObj@fit$sigmax2
var.exp.coef <- growObj@fit$var.exp.coef
sigmax2 <- sigmax2 * exp(2 * (var.exp.coef * mux))
if (length(grep("trunc", tolower(as.character(class(growObj))))) > 0)
sigmax2 <- growObj@fit$sigmax2
}
# range.x <- range(Pmat@meshpoints[loctest[j]] + c(-3.5 * sqrt(sigmax2), +3.5 * sqrt(sigmax2)))
range.x <- range(mux + c(-3.5 * sqrt(sigmax2), +3.5 * sqrt(sigmax2)))
plot(Pmat@meshpoints, Pmat@.Data[, loctest[j]]/h/ps,
type = "n", xlim = range.x, xlab = "Size next", ylab = "Kernel")
#if (j == 1 & kk==1) title("Numerical resolution and growth")
if (j == 2 & kk==1) mtext("Numerical resolution and growth",3,line=4,font=2)
if (j ==1 & kk==1) title("Current Pmatrix")
if (j ==1 & kk==3) title("Increased bins")
for (k in 1:length(Pmat@meshpoints)) {
points(c(Pmat@meshpoints[k]) + c(-h/2, h/2), rep(Pmat@.Data[k,loctest[j]], 2)/h/ps, type = "l",col=cols[kk])
points(rep(Pmat@meshpoints[k] + h/2, 2), c(0, Pmat@.Data[k, loctest[j]]/h/ps), type = "l",
lty = 1,col=cols[kk])
points(rep(Pmat@meshpoints[k] - h/2, 2), c(0,Pmat@.Data[k, loctest[j]]/h/ps), type = "l",
lty = 1,col=cols[kk])
}
if (length(grep("logincr", tolower(as.character(class(growObj))))) == 0 &
length(grep("trunc", tolower(as.character(class(growObj))))) == 0) {
points(testSizes, dnorm(testSizes, mux, sqrt(sigmax2)), type = "l", col = 2)
} else {
if (length(grep("trunc", tolower(as.character(class(growObj))))) > 0) {
require(truncnorm)
points(testSizes, dtruncnorm(testSizes, a = Pmat@meshpoints[loctest[j]],
b = Inf, mean = mux, sd = sqrt(sigmax2)), type = "l",col = 2)
} else {
points(testSizes, dlnorm(testSizes - Pmat@meshpoints[loctest[j]],
mux, sqrt(sigmax2)), type = "l", col = 2)
}
}
legend("topright",legend=paste("size=",round(Pmat@meshpoints[loctest[j]],2)), col = "white",
lty = 1, bty = "n")
}
}
#reset plot characters
on.exit(par(old.par))
}
# =============================================================================
# =============================================================================
# show the news file
IPMpackNews <- function ()
{
newsfile <- file.path(system.file(package = "IPMpack"),
"IPMpackNews")
file.show(newsfile)
}
wpetry/IPMpack2 documentation built on Sept. 29, 2022, 9:41 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Defines functions invLogit sampleIPMOutput sampleIPM sampleVitalRateObj coerceSurvObj coerceGrowthObj addPdfGrowthPic plotSurvModelComp plotGrowthModelComp survModelComp growthModelComp convertIncrement simulateCarlina sampleSequentialIPMs generateData picGrow picSurv
Documented in addPdfGrowthPic coerceGrowthObj coerceSurvObj convertIncrement generateData growthModelComp invLogit picGrow picSurv plotGrowthModelComp plotSurvModelComp sampleIPM sampleIPMOutput sampleSequentialIPMs sampleVitalRateObj simulateCarlina survModelComp
# =============================================================================
# =============================================================================
## FUNCTION FOR MAKING PICTURE OF SURVIVAL DATA AND FITTED SURVIVAL OBJECT ####
## and PICTURE GROWTH DATA AND FITTED GROWTH OBJECT
# ** note that growth may need redefining if new growth objects are defined
# parameters - dataf - the dataframe
# - survObj - a survival object
# - ncuts - the number of cuts used for binning survival
# returns -
picSurv <- function(dataf, survObj, ncuts = 20, makeTitle = "Survival", ...) {
#organize data and plot mean of ncut successive sizes, so trend more obvious
os<-order(dataf$size); os.surv<-(dataf$surv)[os]; os.size<-(dataf$size)[os];
psz<-tapply(os.size,as.numeric(cut(os.size,ncuts)),mean,na.rm=TRUE); #print(psz)
ps<-tapply(os.surv,as.numeric(cut(os.size,ncuts)), mean, na.rm = TRUE);#print(ps)
if (length(grep("covariate",names(survObj@fit$model)))==0) {
#plot data
plot(as.numeric(psz),as.numeric(ps),pch=19,
xlab="Size at t", ylab = "Survival to t+1", main = makeTitle, ...)
#Plot fitted models
points(dataf$size[order(dataf$size)],surv(dataf$size[order(dataf$size)],data.frame(covariate=1),survObj),type="l",col=2)
} else {
plot(as.numeric(psz),as.numeric(ps),
type="n",pch=19,xlab="Size at t", ylab="Survival to t+1",main="Survival",...)
dataf$covariate <- as.factor(dataf$covariate)
levels(dataf$covariate) <- 1:length(unique(dataf$covariate))
os.cov<-(dataf$covariate)[os]
sizes <- dataf$size[!is.na(dataf$size)]; sizes <- sizes[order(sizes)]
ud <- unique(dataf$covariate); ud <- ud[!is.na(ud)]
for (k in 1:length(ud)) {
tp <- os.cov==ud[k]
psz<-tapply(os.size[tp], as.numeric(cut(os.size[tp], ncuts)), mean, na.rm = TRUE); #print(psz)
ps<-tapply(os.surv[tp],as.numeric(cut(os.size[tp],ncuts)),mean,na.rm=TRUE);#print(ps)
points(as.numeric(psz),as.numeric(ps),pch=19,col=k)
newd <- data.frame(size=sizes,size2=sizes^2,size3=sizes^3,
covariate=rep(as.factor(ud[k]),length(sizes)))
if(length(grep("expsize",survObj@fit$formula))==1)
newd$expsize=exp(sizes)
if(length(grep("logsize",survObj@fit$formula))==1)
newd$logsize=log(sizes)
if(length(grep("logsize2",survObj@fit$formula))==1)
newd$logsize=(log(sizes))^2
pred.surv <- predict(survObj@fit,newd,type="response")
points(newd$size,pred.surv,type="l",col=k)
}
}
}
# =============================================================================
# =============================================================================
# Function to make pic of growth object fit with data
#
# parameters - dataf - the dataframe
# - growObj - a growth object
#
# returns -
#
picGrow <- function(dataf, growObj, mainTitle = "Growth",...) {
predVar <- attr(growObj@fit$terms,"predvars")[[2]] #jess quick fix. this function does not work with changingVar either at the moment
if (class(growObj)=="growthObjTruncIncr") {
predVar <- "incr"
} else {
predVar <- attr(growObj@fit$terms,"predvars")[[2]]
}
if(predVar == "sizeNext") {
plot(dataf$size, dataf$sizeNext,pch=19, xlab="Size at t", ylab="Size at t+1", main = mainTitle,...)
abline(a = 0, b = 1)
}else{
dataf$incr <- dataf$sizeNext - dataf$size
plot(dataf$size, dataf$incr, pch = 19, xlab = "Size at t", ylab="Size increment", main = mainTitle,...)
abline(a = 0, b = 0)
}
if (length(grep("covariate", names(growObj@fit$model))) > 0) {
#convert to 1:n for indexing later and to relate to discrete
dataf$covariate <- as.factor(dataf$covariate)
levels(dataf$covariate) <- 1:length(unique(dataf$covariate))
ud <- unique(dataf$covariate); ud <- ud[!is.na(ud)]
for (k in 1:length(ud)) {
tp <- dataf$covariate == ud[k]
points(dataf$size[tp], dataf$sizeNext[tp], pch = 19, col = k)
}
ud <- as.factor(ud)
} else {
ud <- 0
}
sizes <- dataf$size[!is.na(dataf$size)]; sizes <- sizes[order(sizes)]
for (k in 1:length(ud)) {
newd <- data.frame(size = sizes, size2 = sizes ^ 2,size3 = sizes ^ 3,
covariate = as.factor(rep(ud[k],length(sizes))))
if(length(grep("expsize", names(growObj@fit$coefficients))) == 1)
newd$expsize = exp(sizes)
if(length(grep("logsize", names(growObj@fit$coefficients))) == 1)
newd$logsize = log(sizes)
if(length(grep("logsize2", names(growObj@fit$coefficients))) == 1)
newd$logsize=(log(sizes))^2
if (length(grep("decline",tolower(as.character(class(growObj)))))>0 |
length(grep("trunc",tolower(as.character(class(growObj)))))>0) {
pred.size <- .predictMuX(growObj, newd, covPred = k)
} else {
pred.size <- predict(growObj@fit,newd,type = "response")
}
if (length(grep("incr", tolower(as.character(class(growObj))))) == 0) {
points(sizes, pred.size, type = "l", col = k + 1)
} else {
if (length(grep("logincr", tolower(as.character(class(growObj))))) > 0) {
points(sizes, sizes + exp(pred.size), type = "l", col = k + 1)
} else {
lines(sizes, pred.size, col = k + 1)
}
}
}
}
# =============================================================================
# =============================================================================
## FUNCTIONS FOR SIMULATING DATA ############################
## Return a data-frame with the headings used in the 'makeObject' functions
## Generate a simple data-frame for only continuous covariates
# with a total of 1000 measurements with columns called
# size, sizeNext, surv, covariate, covariateNext, fec,
#
#
generateData <- function(nSamp=1000, type="simple"){
if (nSamp<100) stop("Error: nSamp in generateData should be at least 100")
if (type=="simple"){
covariate <- sample(0:1, size=nSamp, replace=TRUE, prob = c(0.2, 0.8))
covariateNext <- sample(0:1, size=nSamp, replace=TRUE, prob = c(0.8, 0.2))
size <- rnorm(nSamp,5,2)
#size <- exp(rnorm(1000, -1, 1.1))
sizeNext <- 1+0.8*size-0.9*covariate+rnorm(nSamp,0,1)
seedlings <- sample(1:nSamp,size=100,replace=TRUE)
size[seedlings] <- NA; sizeNext[seedlings] <- rnorm(100,2,0.5)
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 <- fec*10
stage <- stageNext <- rep("continuous",nSamp)
stage[is.na(size)] <- NA
stageNext[is.na(sizeNext)] <- "dead"
dataf <- data.frame(size=size,sizeNext=sizeNext,surv=surv,
covariate=covariate,covariateNext=covariateNext,
fec=fec, stage=stage,stageNext=stageNext)
dataf$sizeNext[dataf$surv==0] <- NA
}
if (type=="discrete") dataf <-.generateDataDiscrete(nSamp=nSamp)
if (type=="stochastic") dataf <-.generateDataStoch(nSamp=nSamp)
if (type!="simple" & type!="stochastic" & type!="discrete"){
stop("unsupported type: supported types include simple, stochastic, or discrete")
}
return(dataf)
}
# =============================================================================
# =============================================================================
## FUNCTION FOR MAKING A LIST OF IPMS ############################################
# to do for stoch env with a single discrete covariate. ##########################
sampleSequentialIPMs <- function(dataf, nBigMatrix=10, minSize=-2,maxSize=10,
integrateType="midpoint", correction="none",
explSurv=surv~size+size2+covariate,
explGrow=sizeNext~size+size2+covariate,
regType="constantVar",explFec=fec~size,Family="gaussian",
Transform="none",fecConstants=data.frame(NA)) {
#convert to 1:n for indexing later
dataf$covariate <- as.factor(dataf$covariate)
levels(dataf$covariate) <- 1:length(unique(dataf$covariate))
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$covariate)
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(covariate=as.factor(k)),
fecObj = fv1,integrateType=integrateType, correction=correction)
tpS <- makeIPMPmatrix(nBigMatrix = nBigMatrix, minSize = minSize,
maxSize = maxSize, chosenCov = data.frame(covariate=as.factor(k)),
growObj = gr1, survObj = sv1,
integrateType=integrateType, correction=correction)
IPM.list[[k]] <- tpF+tpS
}
return(IPM.list)
}
# =============================================================================
# =============================================================================
### simulation Carlina ######################################################
## Function to simulate something a bit like Carlina
##
## Parameters - nSamp - starting pop sizes
# - nYrs - no years in simulation
# - nSampleYrs - no years to extract for the data
# - ... bunch of parameters
# - meanYear - means for year effects
# - matVarYear - variance covariances for year effects
# - densDep - model density dependence in seedling establishment or not
#
# Returns - list including dataf - data-frame of data
# - various of the simulation parameters for convenience
#
# PARAMETERS FROM TABLE 2 IN REES AND ELLNER 2009
simulateCarlina <- function(nSamp=200,nYrs=1000,nSampleYrs=15,
m0=-1.37,ms=0.59,
b0=-12.05,bs=3.64,
A=-1,B=2,
ag=1.14,bg=0.74,sig=0.29,
mean.kids=3.16,sd.kids=0.5,
meanYear=c(0,0,0),
matVarYear=matrix(c(1.03,0,0,0,0.037,0.041,0,0.041,0.075),3,3),
varA=0,varB=0,
densDep=TRUE,maxPerYr=1000,maxStoreSeedlingsPerYr=200,sizes=c()) {
require(mvtnorm)
#initiate and set up year index
if (length(sizes)==0) sizes <- rnorm(nSamp,3,0.5)
startYr <- (nYrs-nSampleYrs)
#recruits
nrec <- c(20,42,12,17,8,19,58,45,44,2,56,25,75,92,94,6,4,34,104)
#total plants
totpl <- c(21,57,47,25,25,33,88,94,97,26,85,80,122,175,160,10,6,189)
#set up dataframe
dataf <-matrix(NA,(maxPerYr+maxStoreSeedlingsPerYr)*nYrs,13)
maxPop <- nrow(dataf)
n.per.yr <- rep(NA,nYrs)
trueGrow <- rep(NA,nYrs)
count <- 0
#stoch sims
if (sum(matVarYear)!=0) {
tmp <- rmvnorm(nYrs,mean=meanYear,sigma=matVarYear)
} else {
tmp <- matrix(0,nYrs,3)
}
if (varA!=0) tmpA <- rnorm(nYrs,mean=A,sd=sqrt(varA))
if (varB!=0) tmpB <- rnorm(nYrs,mean=B,sd=sqrt(varB))
for (t in 1:nYrs) {
#yr effects
nSeedlings <- sample(nrec,size=1,replace=FALSE)
#print(tmp)
n.per.yr[t] <- length(sizes)
m.year <- tmp[t,1]
cg.year <- tmp[t,2]
b.year <- tmp[t,3]
if (varA!=0) A.year <- tmpA[t] else A.year <- A
if (varB!=0) B.year <- tmpB[t] else B.year <- B
#print(c(m.year,cg.year,b.year))
if (n.per.yr[t]>0) {
#survival
sx <- 1*(invLogit(m0+ms*sizes+m.year)>runif(n.per.yr[t]))
#flowering
fx <- 1*(invLogit(b0+bs*sizes)>runif(n.per.yr[t]))
#fertility
seedsx <- exp(A.year+B.year*sizes)*fx*sx
#seedsx[seedsx>0] <- rpois(sum(seedsx>0),seedsx[seedsx>0])
} else {
sx <- fx <- seedsx <- c()
print(c("extinct in year ", t))
}
if (densDep) pEst <- min(nSeedlings/max(sum(seedsx, na.rm=TRUE),1),1) else pEst <- 1
babies <- rnorm(ceiling(pEst*max(sum(seedsx,na.rm=TRUE),1)),
mean.kids+b.year,sd.kids)
#if (count>0) print(c("yr",t-startYr,pEst*max(sum(seedsx,na.rm=TRUE),1),length(babies)))
#will end up with nrec babies at least in density dependent case
if (length(babies)<nSeedlings & densDep) nSeedlings <- length(babies)
#growth
sizeNext <- rnorm(length(sizes),ag+bg*sizes+cg.year,sig)
#remove dead or flowered
fx[sx==0] <- NA
sizeNext[sx==0 | fx==1] <- NA
#storage
if (t > startYr) {
#print(t)
#don't 'do it at all if too big (so as not to do only adults, etc)
if ((count+length(sizes))>maxPop) { print("large pop size, breaking");break()}
if ((count+length(sizes)+length(babies))>maxPop) { print("large pop size, breaking");break()}
chs <- (count+1):(count+length(sizes))
dataf[chs,1] <- sizes
dataf[chs,2] <- sizeNext
dataf[chs,3] <- sx
dataf[chs,4] <- fx
dataf[chs,5] <- seedsx
dataf[chs,6] <- t
dataf[chs,7] <- nSeedlings
dataf[chs,8] <- m.year
dataf[chs,9] <- cg.year
dataf[chs,10] <- b.year
dataf[chs,12] <- A.year
dataf[chs,13] <- B.year
count <- count + length(sizes)
nbabes.store <- min(length(babies),maxStoreSeedlingsPerYr)
chs <- (count+1):(count+nbabes.store)
dataf[chs,2] <- babies[1:nbabes.store]
dataf[chs,6] <- t
dataf[chs,7] <- nSeedlings
dataf[chs,8] <- m.year
dataf[chs,9] <- cg.year
dataf[chs,10] <- b.year
dataf[chs,11] <- "sexual"
dataf[chs,12] <- A.year
dataf[chs,13] <- B.year
count <- count + nbabes.store
#print(nbabes.store)
}
#new pop - i.e. those that grew, survive, and didn't flower; as well as babies
sizes <- c(sizeNext[sx==1 & fx==0 & !is.na(fx)],babies)
if (length(sizes)==0) print("extinct")
#get true growth rate
trueGrow[t] <- log(length(sizes)/n.per.yr[t])
#cull size at t t down to maxPerYr, if not density dependent
if (!densDep) sizes <- sample(sizes,size=maxPerYr, replace=TRUE)
}
list.par <- list(m0=m0,ms=ms,
b0=b0,bs=bs,
A=A,B=B,varB=varB,
ag=ag,bg=bg,sig=sig,
mean.kids=mean.kids,sd.kids=sd.kids,
meanYear=meanYear,
matVarYear=matVarYear,
nrec=nrec)
dataf <- dataf[1:count,]
dataf <- data.frame(dataf,stringsAsFactors = FALSE)
#print(table(dataf[,6]))
#plot(cumsum(trueGrow[startYr:nYrs])/(1:length(trueGrow[startYr:nYrs])),type="l")
vartrueGrow <- var(trueGrow[startYr:nYrs],na.rm=TRUE)
meantrueGrow <- mean(exp(trueGrow[startYr:nYrs]),na.rm=TRUE)
trueGrow <- mean(trueGrow[startYr:nYrs],na.rm=TRUE)
#abline(h=trueGrow)
colnames(dataf) <- c("size","sizeNext","surv","flower","fec",
"year","nSeedlings","m.year","cg.year","b.year","offspringNext",
"A.year","B.year")
dataf$size <- as.numeric(dataf$size)
dataf$sizeNext <- as.numeric(dataf$sizeNext)
dataf$surv <- as.numeric(dataf$surv)
dataf$flower <- as.numeric(dataf$flower)
dataf$fec <- as.numeric(dataf$fec)
dataf$nSeedlings <- as.numeric(dataf$nSeedlings)
dataf$year <- as.numeric(dataf$year)
dataf$cg.year <- as.numeric(dataf$cg.year)
dataf$m.year <- as.numeric(dataf$m.year)
dataf$b.year <- as.numeric(dataf$b.year)
dataf$A.year <- as.numeric(dataf$A.year)
dataf$B.year <- as.numeric(dataf$B.year)
#print(table(dataf$year))
dataf$fec[dataf$fec==0] <- NA
return(list(dataf=dataf,meanYear=meanYear,matVarYear=matVarYear,
list.par=list.par,trueGrow=trueGrow,
vartrueGrow=vartrueGrow,meantrueGrow=meantrueGrow))
}
# =============================================================================
# =============================================================================
## Convert Increment - where exact dates of census vary but some multiplier of yearly increments
## are desired; this function takes a data-frame (with columns size, sizeNext,
## and, importantly exactDate, exactDatel)
## and returns a data-frame with sizeNext modified proportional to the time
## elapsed in the desired yaerly increments, adding an additional column denoted 'increment'.
#
# Parameters - dataf - a dataframe with headings size, sizeNext, exactDate,exactDatel
# - nYrs - the number of years between censuses desired (e.g. for CTFS data, 5 years)
#
# Returns - a dataframe with the same headings
convertIncrement <- function(dataf, nYrs=1) {
incr <- dataf$sizeNext - dataf$size
if (class(dataf$exactDatel) == "Date") {
timeElapsed <- (difftime(dataf$exactDatel,dataf$exactDate)[[1]])/(365*nYrs)
} else {
timeElapsed <- (dataf$exactDatel-dataf$exactDate)/(365*nYrs)
}
incrNew <- incr / timeElapsed
dataf$sizeNext <- dataf$size + incrNew
dataf$incr <- dataf$incrNew
return(dataf)
}
# =============================================================================
# =============================================================================
# Function to compare model fits for growth and survival objects built with different linear combinations of covariates.
# Growth can have multiple response forms. Uses .makeCovDf. Separate plot functions.
#
#
# Returns - list with a summary table of covariates and scores, and a sub-list of growth or survival objects.
#
growthModelComp <- function(dataf,
expVars = c(sizeNext~1, sizeNext~size, sizeNext~size + size2),
regressionType = "constantVar",
testType = "AIC",
makePlot = FALSE,
mainTitle = "",
plotLegend = TRUE,
legendPos = "topright", ...) {
varN <- length(expVars)
typeN <- length(regressionType)
treatN <- varN * typeN
summaryTable <- data.frame()
grObj <- vector("list", length = treatN)
i <- 1
for(v in 1:varN) {
for(t in 1:typeN) {
grObj[[i]] <- makeGrowthObj(dataf = dataf, regType = regressionType[t], Formula = expVars[[v]])
if (length(grep("decline",tolower(as.character(class(grObj[[i]])))))>0) {
summaryTable <- rbind(summaryTable, cbind(as.character(unlist(expVars))[v], regressionType[t], match.fun(testType)(grObj[[i]]@fit$fit)))
} else {
summaryTable <- rbind(summaryTable, cbind(as.character(unlist(expVars))[v], regressionType[t], match.fun(testType)(grObj[[i]]@fit)))
}
i <- i + 1
}
}
summaryTable <- as.data.frame(summaryTable)
names(summaryTable) <- c("Exp. Vars", "Reg. Type", testType)
outputList <- list(summaryTable = summaryTable, growthObjects = grObj)
# PLOT SECTION #
if(makePlot == TRUE) {
plotGrowthModelComp(grObj = grObj, summaryTable = summaryTable, dataf = dataf, expVars = expVars, testType = "AIC", plotLegend = plotLegend, mainTitle = mainTitle, legendPos, ...)
}
return(outputList)
}
# =============================================================================
# =============================================================================
survModelComp <- function(dataf,
expVars = c(surv~1, surv~size, surv~size + size2),
testType = "AIC",
makePlot = FALSE,
mainTitle = "", ncuts = 20,
plotLegend = TRUE,
legendPos = "bottomleft", ...) {
varN <- length(expVars)
treatN <- varN
summaryTable <- data.frame()
svObj <- vector("list", length = treatN)
i <- 1
for(v in 1:varN) {
svObj[[i]] <- makeSurvObj(dataf = dataf, Formula = expVars[[v]])
summaryTable <- rbind(summaryTable, cbind(as.character(unlist(expVars))[v], match.fun(testType)(svObj[[i]]@fit)))
i <- i + 1
}
summaryTable <- as.data.frame(summaryTable)
names(summaryTable) <- c("Exp. Vars", testType)
outputList <- list(summaryTable = summaryTable, survObjects = svObj)
# PLOT SECTION #
if(makePlot == TRUE) {
## this is the surv picture
plotSurvModelComp(svObj = svObj, summaryTable = summaryTable, dataf = dataf, expVars = expVars, testType = "AIC", plotLegend = plotLegend, mainTitle = mainTitle, ncuts = ncuts, legendPos = legendPos, ...)
}
return(outputList)
}
# =============================================================================
# =============================================================================
# Plot functions for model comparison. Plots the series of fitted models for growth and survival objects.
# Can plot a legend with the model covariates and model test criterion scores (defaults to AIC).
plotGrowthModelComp <- function(grObj, summaryTable, dataf, expVars, testType = "AIC", plotLegend = TRUE, mainTitle = "", legendPos = "topright", ...) {
treatN <- length(grObj)
sizeSorted <- unique(sort(dataf$size))
if(length(grep("sizeNext", unlist(as.character(expVars))))) {
y.lab <- "Size at t + 1"
dataSizeNext <- dataf$sizeNext
}
if(length(grep("incr", unlist(as.character(expVars))))) {
y.lab <- "Growth"
dataSizeNext <- dataf$sizeNext - dataf$size
}
if(length(grep("logincr", unlist(as.character(expVars))))){
y.lab <- "log(growth)"
dataSizeNext <- log(dataf$sizeNext - dataf$size)
}
plot(dataf$size, dataSizeNext, pch = 19, xlab = "Size at t", ylab = y.lab, main = mainTitle, cex = 0.8,...)
for(p in 1:treatN) {
#PROBLEM HERE
expVarsHere <- paste(attr(terms((expVars[[p]])),"term.labels"),collapse="+")
newd <- .makeCovDf(sizeSorted, expVarsHere)
if (length(grep("decline",tolower(as.character(class(grObj[[p]])))))>0) {
pred.size <- .predictMuX(grObj[[p]], newd)
} else { pred.size <- predict(grObj[[p]]@fit, newd, type = "response")}
lines(sizeSorted, pred.size, type = "l", col = (p + 1))
}
if(plotLegend) {
legend(legendPos, legend = sprintf("%s: %s = %.1f", expVars, testType, as.numeric(as.character(summaryTable[,3]))), col = c(2:(p + 1)), lty = 1, xjust = 1, bg = "white")
}
}
# =============================================================================
# =============================================================================
plotSurvModelComp <- function(svObj, summaryTable, dataf, expVars, testType = "AIC", plotLegend = TRUE, mainTitle = "", ncuts = 20, legendPos = "bottomleft", ...) {
treatN <- length(svObj)
#ncuts <- 20 # survival bins
os <- order(dataf$size) # order size
osSurv <- (dataf$surv)[os] # order survival data according to size
osSize<-(dataf$size)[os] # ordered size data
binnedSize <- tapply(osSize, as.numeric(cut(osSize, breaks=ncuts)), mean, na.rm = TRUE); # bin Size data
binnedSurv <- tapply(osSurv, as.numeric(cut(osSize, breaks=ncuts)), mean, na.rm = TRUE) #bin Survival probabilities
plot(binnedSize, binnedSurv, pch = 19, xlab = "Size at t", ylab = "Survival to t + 1", main = mainTitle, cex = 0.8,...)
for(p in 1:treatN) {
expVarsHere <- paste(attr(terms(expVars[[p]]),"term.labels"),collapse="+")
newd <- .makeCovDf(osSize, expVarsHere[p])
lines(dataf$size[order(dataf$size)], surv(dataf$size[os], data.frame(covariate=1), svObj[[p]]), col = (p + 1))
}
if(plotLegend) {
legend(legendPos, legend = sprintf("%s: %s = %.1f", expVars, testType, as.numeric(as.character(summaryTable[,2]))), col = c(2:(p + 1)), lty = 1, xjust = 1, bg = "white")
}
}
# =============================================================================
# =============================================================================
## Function to add pdf to model comparison pics
addPdfGrowthPic <- function(respType = "sizeNext",
sizesPlotAt=c(20,50,60),
sizeRange=c(20,400),
incrRange=c(-10,50),
scalar=100,
growthObjList,
cols=1:5,
cov=data.frame(covariate=1),
minShow=1e-2,
jitt=2, #how far apart should pdfs be if you are plotting several
...){
nval <- length(growthObjList)
sizes <- seq(sizeRange[1],sizeRange[2],length=500)
incr <- seq(incrRange[1],incrRange[2],length=500)
logincr <- log(incr); logincr[!is.finite(logincr)] <- NA
for (j in 1:nval) {
#print(j)
for (k in 1:length(sizesPlotAt)) {
if (respType=="sizeNext") {
pred <- growth(sizesPlotAt[k],sizes,cov,growthObjList[[j]])*scalar
pred[pred<minShow] <- NA
points(sizesPlotAt[k]+pred+jitt*(j-1),
sizes,type="l", col=cols[j],...)
}
if (respType=="incr") {
pred <- growth(sizesPlotAt[k],sizesPlotAt[k]+incr,cov,growthObjList[[j]])*scalar
pred[pred<minShow] <- NA
points(sizesPlotAt[k]+pred+jitt*(j-1),
incr,type="l", col=cols[j],...)
}
if (respType=="logincr") {
pred <- growth(sizesPlotAt[k],sizesPlotAt[k]+logincr,cov,growthObjList[[j]])*scalar
pred[pred<minShow] <- NA
#print(range(pred,na.rm=TRUE))
points(sizesPlotAt[k]+pred+jitt*(j-1),
logincr,type="l", col=cols[j],...)
}
}}
}
# =============================================================================
# =============================================================================
## Function to coerce Growth object to parameters and variance desired
coerceGrowthObj <- function(growthObj, coeff, sd){
if (length(growthObj@fit$coefficients) !=length(coeff)) {
print("warning: number of desired coefficients to not match number of growth object coefficients")
}
growthObj@fit$coefficients[] <- as.numeric(coeff)
growthObj@sd[] <- as.numeric(sd)
return(growthObj)
}
# =============================================================================
# =============================================================================
## Function to coerce Survival object to parameters desired
coerceSurvObj <- function(survObj, coeff){
if (length(survObj@fit$coefficients) !=length(coeff)) print("warning: number of desired coefficients to not match number of growth object coefficients")
survObj@fit$coefficients[] <- as.numeric(coeff)
return(survObj)
}
# =============================================================================
# =============================================================================
## Parametric boostrap of vital rate objects
# created by cory on 6-4-13 as a modification of getListRegObjs()
sampleVitalRateObj <- function(Obj,nSamp=100,nDiscreteGrowthTransitions=NULL,nDiscreteOffspringTransitions=NULL,nOffspring=NULL) {
require(mvtnorm)
require(MCMCpack)
objList <- list()
#generate new set parameters from mvn
# what is the difference between offspring splitter transitions and discrete
# for discrete transition matrices objects
if( sum(class(Obj)==c("discreteTrans","discreteTransInteger")) ) {
if(is.null(nDiscreteGrowthTransitions)){
stop('Error: Please specify the total number of transitions that were used to estimate the discrete transitions in Obj@discreteTrans to allow the function to estimate the appropriate variance for resampling them')
}
if(!is.null(nDiscreteGrowthTransitions)){
for (j in 1:nSamp) {
objList[[j]] <- Obj
for (k in 1:ncol(Obj@discreteTrans)){
objList[[j]]@discreteTrans[,k] <- rdirichlet(1, nDiscreteGrowthTransitions * as.numeric(Obj@discreteTrans[,k]))
}
}
}
}
# for growth and survival objects
if( !sum(class(Obj)==c("fecObj","fecObjInteger","discreteTrans", "discreteTransInteger")) ) {
newpar <- rmvnorm(nSamp, mean = Obj@fit$coefficients, sigma = vcov(Obj@fit))
for (j in 1:nSamp) {
objList[[j]] <- Obj
objList[[j]]@fit$coefficients <- newpar[j,]
}
}
# for fecundity objects
if( sum(class(Obj)==c("fecObj","fecObjInteger")) ) {
for (j in 1:nSamp) {
# sample fecundity regression parameters
for (k in 1:length(Obj@fitFec)) {
newpar <- rmvnorm(1, mean = Obj@fitFec[[k]]$coefficients,
sigma = vcov(Obj@fitFec[[k]]))
objList[[j]] <- Obj
objList[[j]]@fitFec[[k]]$coefficients <- newpar
}
# sample discrete transitions from dirichlet
if(ncol(Obj@offspringSplitter)>1){ # only if there are discrete fec stages
if(is.null(nDiscreteOffspringTransitions)){
stop('Error: Please specify the total number of discrete transitions that were used to estimate the fecundity transition in Obj@offspringSplitter to allow the function to estimate the appropriate variance for resampling them')
}
if(!is.null(nDiscreteOffspringTransitions)){
objList[[j]]@offspringSplitter[] <- rdirichlet(1, nDiscreteOffspringTransitions * as.numeric(Obj@offspringSplitter))
}
}
# sample the offspring sizes
# TODO: OffspringRel currently does not store all the information from the regression, which it should so that this works with offspring models that have more than an intercept
if(is.null(nOffspring)){
stop('Error: Please specify the total number of offspring that were used to estimate the offspring size distribtuion in Obj@offspringRel to allow the function to estimate the appropriate variance for resampling them')
}
# TODO: Make work with the new OffspringObjs
if(!is.null(nOffspring)){
newpar3 <- rnorm(1,coefficients(Obj@offspringRel)[1], Obj@sdOffspringSize/sqrt(nOffspring))
objList[[j]]@offspringRel$coefficients[] <- newpar3
newpar4 <- rnorm(1,Obj@sdOffspringSize, sqrt(2*Obj@sdOffspringSize^4)/(nOffspring-1))
objList[[j]]@sdOffspringSize <- newpar4
}
}
}
return(objList)
}
# =============================================================================
# =============================================================================
## Use list of vital rate objects to make List of IPMS
# created by cory on 6-4-13 as a modification of getListIPMs()
# TODO: MAKE WORK FOR OFFSPRING OBJS
sampleIPM<- function(
growObjList=NULL,survObjList=NULL,fecObjList=NULL,
offspringObjList=NULL, discreteTransList=1,
nBigMatrix,minSize,maxSize,
covariates=FALSE,envMat=NULL,
integrateType="midpoint",correction="none",warn=TRUE) {
if(!is.null(offspringObjList)){
stop('Sorry, lists of offspring objects are not yet supported. However, specify offspring size distributions via makeFecObj() is supported. Please set offspringObjList=NULL, or modify the function to accept offspring objects and email it to us.')
}
# make the same number of samples for each vital rate object
n.samples <- max(c(length(growObjList),length(survObjList),length(fecObjList)))
if(!is.null(growObjList) & length(growObjList)<n.samples){
if(warn) warning('Length of growth object list is less than the length of another vital rate object list, so some members of the growth object list have been repeated.')
growthObjList <- sample(growObjList,size=n.samples,replace=T)
}
if(!is.null(survObjList) & length(survObjList)<n.samples ){
if(warn) warning('Length of survival object list is less than the length of another vital rate object list, so some members of the survival object list have been repeated.')
survObjList <- sample(survObjList,size=n.samples,replace=T)
}
if(!is.null(fecObjList) & length(fecObjList)<n.samples ){
if(warn) warning('Length of fec object list is less than the length of another vital rate object list, so some members of the fec object list have been repeated.')
fecObjList <- sample(fecObjList,size=n.samples,replace=T)
}
if(!is.null(growObjList) & length(discreteTransList)<n.samples ){
# 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=n.samples,replace=T)
}
if(!is.null(growObjList)){
PmatrixList <- .makeListPmatrix(growObjList=growObjList, survObjList=survObjList, nBigMatrix=nBigMatrix, minSize=minSize, maxSize=maxSize, cov=covariates, envMat=envMat,discreteTransList=discreteTransList, integrateType=integrateType, correction=correction)
matrixList <- PmatrixList
}
if(!is.null(fecObjList)){
FmatrixList <- .makeListFmatrix(fecObjList, nBigMatrix=nBigMatrix, minSize=minSize, maxSize=maxSize, cov=covariates, envMat=envMat, integrateType=integrateType, correction=correction)
matrixList <- FmatrixList
}
if(!is.null(growObjList) & !is.null(fecObjList)){
for(i in 1:n.samples){
matrixList[[i]] <- PmatrixList[[i]] # to get all the Pmatrix slots
matrixList[[i]]@.Data <- PmatrixList[[i]]+FmatrixList[[i]]
}
}
return(matrixList)
}
# =============================================================================
# =============================================================================
## Use list of IPMs to make List of IPM output
# created by cory on 6-4-13 as a modification of getIPMOutput() and getIPMOutputDirect()
sampleIPMOutput <- function(IPMList=NULL,PMatrixList=NULL,passageTimeTargetSize=c(), sizeToAgeStartSize=c(),sizeToAgeTargetSize=c(),warn=TRUE){
# removed from old version of getIPMOutputDirect
# # cov=FALSE, envMat=NULL,
# # chosenCov = data.frame(covariate = 1),
# # onlyLowerTriGrowth=FALSE)
# if only a Pmatrix is supplied, just call it IPMList for simplicity
if(is.null(IPMList)) IPMList=PMatrixList
if ( is.null(passageTimeTargetSize) ) {
if(warn) print("no target size for passage time provided; taking meshpoint median")
passageTimeTargetSize <- median(IPMList[[1]]@meshpoints)
}
if ( is.null(sizeToAgeStartSize) ) {
if(warn) print("no starting size for size to age provided; taking minimum size")
sizeToAgeStartSize <- min(IPMList[[1]]@meshpoints)
}
if ( is.null(sizeToAgeTargetSize) ) {
if(warn) print("no target size for size to age provided; taking meshpoint values")
sizeToAgeTargetSize <- IPMList[[1]]@meshpoints
}
nSamps <- length(IPMList)
stableStage <- LE <- pTime <- matrix(NA,nSamps,length(IPMList[[1]]@.Data[1,]))
lambda <- rep(NA,nSamps)
resAge <- resSize <- matrix(NA,nSamps,length(sizeToAgeTargetSize))
h1 <- diff(IPMList[[1]]@meshpoints)[1]
for (k in 1:nSamps) {
IPM <- IPMList[[k]]
LE[k,]<-meanLifeExpect(IPM)
pTime[k,]<-passageTime(passageTimeTargetSize,IPM)
if (is.null(PMatrixList)) {
lambda[k] <- Re(eigen(IPM)$value[1])
stableStage[k,] <- eigen(IPM)$vector[,1]
#normalize stable size distribution
stableStage[k,] <- stableStage[k,]/sum(stableStage[k,])
}
# get size to age results
res2 <- sizeToAge(Pmatrix=IPM,startingSize=sizeToAgeStartSize,
targetSize=sizeToAgeTargetSize)
resAge[k,] <- res2$timeInYears
resSize[k,] <- res2$targetSize
}
return(list(LE=LE,pTime=pTime,lambda=lambda,stableStage=stableStage,
resAge=resAge,resSize=resSize,meshpoints=IPM@meshpoints))
}
# =============================================================================
# =============================================================================
invLogit <- function(x) {
u <- exp(pmin(x, 50))
return(u / (1 + u))
}
# =============================================================================
# =============================================================================
# For a single Pmatrix (!not compound and no discrete stages!), this functions makes a series
# of diagnostic plots - this is defined for growthObj,
# growthObjIncr - modification required
# if other objects used
#
# Parameters - the Pmatrix
# - growObj - growth object used to build it
# - survObj - survival object used to build it
# - dff - the data from which it was built
# - integrateType - "midpoint", or "cumul" - should
# correspond to what the IPM was built with
# - do you want to implement the corrections for integration?
# Returns -
#
diagnosticsPmatrix <- function (Pmatrix, growObj, survObj, dff=NULL,
integrateType = "midpoint",
correction = "none", cov = data.frame(covariate = 1),
sizesToPlot=c(),extendSizeRange=c()) {
print("Range of Pmatrix is ")
print(range(c(Pmatrix)))
if (Pmatrix@meshpoints[1] > 0)
new.min <- Pmatrix@meshpoints[1]/2
else new.min <- Pmatrix@meshpoints[1] * 1.5
new.max <- 1.5 * max(Pmatrix@meshpoints)
if (length(extendSizeRange)>0 & length(extendSizeRange)!=2) print("require two values for extendSizeRange, reflecting upper and lower limits")
if (length(extendSizeRange)>0) { new.min <- extendSizeRange[1]; new.max <- extendSizeRange[2]}
#colours for 1) current; 2) bigger size range; 3) bigger no bins
cols <- c("black","tomato","darkblue")
ltys <- c(1,1,3)
#matrix with bigger size range
Pmatrix1 <- makeIPMPmatrix(nEnvClass = 1, nBigMatrix = length(Pmatrix@meshpoints),
minSize = new.min, maxSize = new.max,
chosenCov = cov, growObj = growObj, survObj = survObj,
integrateType = integrateType, correction = correction)
if (sum(is.na(Pmatrix1)) > 0) {
print("Pmatrix with extended size range returns NAs; changing these to 0, and putting survival value onto diagonal for columns that sum to zero")
Pmatrix1[is.na(Pmatrix1)] <- 0
bad <- which(colSums(Pmatrix1) == 0, arr.ind = TRUE)
if (length(bad) > 0)
Pmatrix1[cbind(bad, bad)] <- surv(size = Pmatrix1@meshpoints[bad],
cov = cov, survObj = survObj)
}
#matrix with bigger number of bins
Pmatrix2 <- makeIPMPmatrix(nEnvClass = 1, nBigMatrix = floor(length(Pmatrix@meshpoints) *
1.5), minSize =Pmatrix@meshpoints[1], maxSize = max(Pmatrix@meshpoints),
chosenCov = cov, growObj = growObj, survObj = survObj,
integrateType = integrateType, correction = correction)
if (sum(is.na(Pmatrix2)) > 0) {
print("Pmatrix with extended number of bins returns NAs; changing these to 0, and putting survival value onto diagonal for columns that sum to zero")
Pmatrix2[is.na(Pmatrix2)] <- 0
bad <- which(colSums(Pmatrix2) == 0, arr.ind = TRUE)
if (length(bad) > 0)
Pmatrix2[cbind(bad, bad)] <- surv(size = Pmatrix2@meshpoints[bad],
cov = cov, survObj = survObj)
}
#start plots - put original Pmatrix in black
#par(mfrow = c(3, 3), bty = "l")
#save plot characters
old.par <- par(no.readonly = TRUE)
par(mfrow = c(1, 3), bty = "l",pty="s", mar=c(5,4,4,1))
if (!is.null(dff)) {
xlims <- range(c(Pmatrix@meshpoints, Pmatrix1@meshpoints,
dff$size, dff$sizeNext), na.rm = TRUE)
a1 <- hist(c(dff$size, dff$sizeNext), xlab = "Sizes observed", axes=FALSE,
ylab = "Frequency", main = "", xlim = xlims, col="lightgrey", border="lightgrey",plot=TRUE)
axis(1); axis(2)
} else {
a1 <- list(); a1$counts <- 1:100
xlims <- range(c(Pmatrix@meshpoints, Pmatrix1@meshpoints))
plot(1:100,type="n",xlab = "Sizes", ylab="", ylim=range(a1$counts), xlim=xlims)
}
lcs <- c(0.8,0.6,0.4)
lcs.x <- mean(xlims)
text(lcs.x ,max(a1$counts)*lcs[1],"Current",pos=3,col=cols[1])
arrows(Pmatrix@meshpoints[1], max(a1$counts)*lcs[1],Pmatrix@meshpoints[length(Pmatrix@meshpoints)], max(a1$counts)*lcs[1],col=cols[1], length=0.1, code=3,lty=ltys[1])
text(lcs.x ,max(a1$counts)*lcs[2],"Extended range",pos=3,col=cols[2])
arrows(Pmatrix1@meshpoints[1], max(a1$counts)*lcs[2],Pmatrix1@meshpoints[length(Pmatrix1@meshpoints)], max(a1$counts)*lcs[2],col=cols[2], length=0.1, code=3,lty=ltys[1])
text(lcs.x ,max(a1$counts)*lcs[3],"Increased no of bins",pos=3,col=cols[3])
arrows(Pmatrix2@meshpoints[1], max(a1$counts)*lcs[3],Pmatrix2@meshpoints[length(Pmatrix2@meshpoints)], max(a1$counts)*lcs[3],col=cols[3], length=0.1, code=3,lty=ltys[1])
#legend("topright", legend = "fitted range", col = "black", lty = 2, bty = "n") ##!!! change this
title("Size range")
#survival sums
lims <- range(c(colSums(Pmatrix), surv(Pmatrix@meshpoints, cov, survObj)))
plot(colSums(Pmatrix), surv(Pmatrix@meshpoints, cov, survObj),
type = "n", xlab = "Surviving in Pmatrix", ylab = "Should be surviving",col=cols[1],lty=ltys[1],
xlim=lims,ylim=lims)
abline(0, 1, lwd=2,col="grey")
points(colSums(Pmatrix), surv(Pmatrix@meshpoints, cov, survObj), type = "l", col = cols[1],lty=ltys[1])
points(colSums(Pmatrix1), surv(Pmatrix1@meshpoints, cov, survObj), type = "l", col = cols[2],lty=ltys[2])
points(colSums(Pmatrix2), surv(Pmatrix2@meshpoints, cov, survObj), type = "l", col = cols[3],lty=ltys[3])
title("Survival")
#text(lims[2],lims[2],"(0,1)",pos=1)
# mtext("Points should lie along the 0,1 line, shown in grey", side=1,outer=TRUE,line=-1)
LE <- meanLifeExpect(Pmatrix)
LE1 <- meanLifeExpect(Pmatrix1)
LE2 <- meanLifeExpect(Pmatrix2)
plot(Pmatrix@meshpoints, LE, type = "l", xlim = range(Pmatrix1@meshpoints),
ylim = range(c(LE, LE1)), xlab = "Sizes", ylab = "Life expectancy",col=cols[1],lty=ltys[1])
points(Pmatrix1@meshpoints, LE1, type = "l", col = cols[2],lty=ltys[2])
points(Pmatrix2@meshpoints, LE2, type = "l", col = cols[3],lty=ltys[3])
legend("topleft", legend = c("Current", "Extended range", "Increased no of bins"), col = cols, lty = ltys, bty = "n")
title("Life expectancy")
#mtext("Increasing size range (red) or number of bins (blue) should not alter life expectancy estimates", side=1,outer=TRUE,line=-1)
if (length(sizesToPlot)==0 & !is.null(dff)) sizesToPlot <- as.numeric(quantile(dff$size,c(0.25, 0.5,0.75),na.rm=TRUE))
if (length(sizesToPlot)==0 & is.null(dff)) sizesToPlot <- as.numeric(quantile(Pmatrix@meshpoints,c(0.25, 0.5,0.75),na.rm=TRUE))
loctest <- rep(NA,length(sizesToPlot))
print("Please hit any key for the next plot")
scan(quiet="TRUE")
par(mfrow = c(2, 3), bty = "l",pty="s")
for (kk in c(1,3)) {
if (kk==1) Pmat <- Pmatrix
if (kk==3) Pmat <- Pmatrix2
h <- diff(Pmat@meshpoints)[1]
testSizes <- seq(min(Pmat@meshpoints), max(Pmat@meshpoints),
length = 5000)
for (j in 1:3) {
loctest[j] <- which(abs(Pmat@meshpoints-sizesToPlot[j])==min(abs(Pmat@meshpoints-sizesToPlot[j])),arr.ind=TRUE)[1]
#print(loctest[j])
ps <- surv(Pmat@meshpoints[loctest[j]], cov, survObj)
newd <- data.frame(size = Pmat@meshpoints[loctest[j]],
size2 = Pmat@meshpoints[loctest[j]]^2,
size3 = Pmat@meshpoints[loctest[j]]^3,
covariate = Pmat@env.index[1])
if (length(grep("expsize", names(growObj@fit$coefficients))))
newd$expsize = log(Pmat@meshpoints[loctest[j]])
if (length(grep("logsize", names(growObj@fit$coefficients))))
newd$logsize = log(Pmat@meshpoints[loctest[j]])
if (length(growObj@fit$model$covariate) > 0)
if (is.factor(growObj@fit$model$covariate))
newd$covariate <- as.factor(newd$covariate)
if (length(grep("decline", tolower(as.character(class(growObj))))) > 0 |
length(grep("trunc", tolower(as.character(class(growObj))))) > 0) {
mux <- .predictMuX(growObj, newd)
} else {
mux <- predict(growObj@fit, newd, type = "response")
}
if (length(grep("incr", tolower(as.character(class(growObj))))) > 0 &
length(grep("logincr", tolower(as.character(class(growObj))))) == 0) {
mux <- Pmat@meshpoints[loctest[j]] + mux
}
if (length(grep("decline", tolower(as.character(class(growObj))))) == 0 &
length(grep("trunc", tolower(as.character(class(growObj))))) == 0) {
sigmax2 <- growObj@sd^2
} else {
sigmax2 <- growObj@fit$sigmax2
var.exp.coef <- growObj@fit$var.exp.coef
sigmax2 <- sigmax2 * exp(2 * (var.exp.coef * mux))
if (length(grep("trunc", tolower(as.character(class(growObj))))) > 0)
sigmax2 <- growObj@fit$sigmax2
}
# range.x <- range(Pmat@meshpoints[loctest[j]] + c(-3.5 * sqrt(sigmax2), +3.5 * sqrt(sigmax2)))
range.x <- range(mux + c(-3.5 * sqrt(sigmax2), +3.5 * sqrt(sigmax2)))
plot(Pmat@meshpoints, Pmat@.Data[, loctest[j]]/h/ps,
type = "n", xlim = range.x, xlab = "Size next", ylab = "Kernel")
#if (j == 1 & kk==1) title("Numerical resolution and growth")
if (j == 2 & kk==1) mtext("Numerical resolution and growth",3,line=4,font=2)
if (j ==1 & kk==1) title("Current Pmatrix")
if (j ==1 & kk==3) title("Increased bins")
for (k in 1:length(Pmat@meshpoints)) {
points(c(Pmat@meshpoints[k]) + c(-h/2, h/2), rep(Pmat@.Data[k,loctest[j]], 2)/h/ps, type = "l",col=cols[kk])
points(rep(Pmat@meshpoints[k] + h/2, 2), c(0, Pmat@.Data[k, loctest[j]]/h/ps), type = "l",
lty = 1,col=cols[kk])
points(rep(Pmat@meshpoints[k] - h/2, 2), c(0,Pmat@.Data[k, loctest[j]]/h/ps), type = "l",
lty = 1,col=cols[kk])
}
if (length(grep("logincr", tolower(as.character(class(growObj))))) == 0 &
length(grep("trunc", tolower(as.character(class(growObj))))) == 0) {
points(testSizes, dnorm(testSizes, mux, sqrt(sigmax2)), type = "l", col = 2)
} else {
if (length(grep("trunc", tolower(as.character(class(growObj))))) > 0) {
require(truncnorm)
points(testSizes, dtruncnorm(testSizes, a = Pmat@meshpoints[loctest[j]],
b = Inf, mean = mux, sd = sqrt(sigmax2)), type = "l",col = 2)
} else {
points(testSizes, dlnorm(testSizes - Pmat@meshpoints[loctest[j]],
mux, sqrt(sigmax2)), type = "l", col = 2)
}
}
legend("topright",legend=paste("size=",round(Pmat@meshpoints[loctest[j]],2)), col = "white",
lty = 1, bty = "n")
}
}
#reset plot characters
on.exit(par(old.par))
}
# =============================================================================
# =============================================================================
# show the news file
IPMpackNews <- function ()
{
newsfile <- file.path(system.file(package = "IPMpack"),
"IPMpackNews")
file.show(newsfile)
}
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.