# R/explosion.R In Rramas: Matrix Population Models

#### Documented in explosion

explosion <-
function(rmas, bootsp=1000){
x<-rmas # TODO: change x for rmas in all function
if(class(x)[1]!="rmas") stop("explosion requires an rmas object
(i.e. a trajectory simulation from projectn)")
if(length(names(x))>0){ # bifurcacion provisional para separar rmas de projectn y de projectn2
x<- x\$vn
}

abundances <- sapply(x, function(rmas) apply(rmas,2,sum))
abundances.min <-round(apply(abundances[-1,],2,max))
abminbot <- as.list(1:bootsp)
# bootstrap maximum abundances
abminbot <- lapply(abminbot,function(x) x <-sample(abundances.min, replace=TRUE))
thresholds <- sort(unique(abundances.min))

# inner function to compute explosion probabilities
decl.prob <- function(abundances.min, thresholds){
cf <-NULL
for(i in 1:length(thresholds)){
# how many times has ocuured an abundance >= treshold[i]?
cf <- c(cf, sum(abundances.min >= thresholds[i]))
}
cf <- cf/length(abundances.min)
cf <- data.frame(Threshold=thresholds, Probability=cf)
return(cf)
}

cf.obs <- decl.prob(abundances.min, thresholds=thresholds)
cf.boot <- lapply(abminbot, decl.prob, thresholds=thresholds)
result<- list(cf.obs=cf.obs, cf.boot=cf.boot,abminbot=abminbot,
main="Explosion/Increase")
class(result)<- c("rmas.risk", class(result))
return(result)
}

## Try the Rramas package in your browser

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

Rramas documentation built on May 2, 2019, 6:52 a.m.