.Rhistory

else {
return(p + netcontributingunits * p / maxunits)
}
}
GetBinaryORFacInhResponse <- function(netcontributingunits, maxunits, p){
# OR response relating number of faciliatory and/or inhibitory species present to probability of focal species presence
if(netcontributingunits > 0){
stepfactor <- 1
}
else if(netcontributingunits < 0){
stepfactor <- -1
}
else{
stepfactor <- 0
}
if(p > 0.5){
return(p + stepfactor * (1 - p))
}
else {
return(p + stepfactor * p)
}
}
GetBinaryNonLinearFacInhResponse <- function(netcontributingunits, posEdges, negEdges, maxunits, p, a, b){
# nonlinear response (beta function with parameters a and b) relating number of faciliatory and/or inhibitory species present to probability of focal species presence
if(p > 0.5){
positivecontribution <- (1 - p) * (pbeta(posEdges / maxunits, a, b))
negativecontribution <- (1 - p) * (pbeta(negEdges / maxunits, a, b))
return(p + positivecontribution - negativecontribution)
}
else {
positivecontribution <- p * (pbeta(posEdges / maxunits, a, b))
negativecontribution <- p * (pbeta(negEdges / maxunits, a, b))
return(p + positivecontribution - negativecontribution)
}
}
TicksResponse <- function(netcontributingunits, posEdges, negEdges, maxunits, p, offset, a, b){
# nonlinear response (beta function with parameters a and b) relating number of faciliatory and/or inhibitory species present to probability of focal species presence
if(p > 0.5){
if(posEdges==(maxunits/2)){
returnvalue <- p
}
else if (posEdges > (maxunits/2)){
positivecontribution <- (p - offset) * pbeta((posEdges - (maxunits/2))/(maxunits/2), a, b)
if((p + positivecontribution) >= 1){
returnvalue <- 1
}
else {
returnvalue <- p + positivecontribution
}
}
else if (posEdges < (maxunits/2)){
positivecontribution <- (p - offset) * pbeta(((maxunits/2) - posEdges)/(maxunits/2), a, b)
returnvalue <- p - positivecontribution
}
else {print("ERROR in TicksResponse()")}
positivecontribution <- p * (pbeta(posEdges / maxunits, a, b))
return(returnvalue)
}
else {
if(posEdges==(maxunits/2)){
returnvalue <- p
}
else if (posEdges > (maxunits/2)){
positivecontribution <- (p - offset) * pbeta((posEdges - (maxunits/2))/(maxunits/2), a, b)
returnvalue <- p + positivecontribution
}
else if (posEdges < (maxunits/2)){
positivecontribution <- (p - offset) * pbeta(((maxunits/2) - posEdges)/(maxunits/2), a, b)
returnvalue <- p - positivecontribution
}
else {print("ERROR in TicksResponse()")}
positivecontribution <- p * (pbeta(posEdges / maxunits, a, b))
return(returnvalue)
}
}
############### BUILD TABLE FOR gRain
getNames <- function(prey){
# get rownames of a state table defined by which species are present and which are not
n <- length(prey)
ch <- matrix("", 2^n, n)
status <- rep("", 2^n)
for(i in 0:((2^n) - 1)){
x <- i
for(j in 1:n){
if(x%%2 == 0){
ch[i+1, j] <- as.character(prey[j])
x <- floor(x/2)
} else {
ch[i+1, j] <- "-"
x <- floor(x/2)
}
status[i + 1] <- paste(status[i + 1], " ", ch[i + 1, j], sep='')
}
}
return(status)
}
BuildTable <- function(M, Pi, offest, model, alpha=1, beta=1){
# form input for getMarginals_gRain()  including state tables for each species
use_sf <- 1  # USER-DEFINED OPTION
if (use_sf==1){significant_figures <- 7}
PiVector <- Pi
S <- length(PiVector)
Table <- list()
Table$PiVector <- PiVector
Table$M <- M
for (i in 1:S){
variable_name <- paste("V", i, sep='')
prey <- which(M[, i] != 0)
n <- length(prey)
if (n==0){
# species is basal
X <- matrix(0, 1, 1)
if (use_sf==1){X[1, 1] <- signif(PiVector[i], digits=significant_figures)}
else{X[1, 1] <- PiVector[i]}
Table[[variable_name]] <- X
}
else{
# species is a consumer
X <- matrix(0, 2^sum(M[, i] != 0), 1)
rownames(X) <- getNames(prey)  # top-to-bottom: all influencing species absent to all present
if(model=="linear" || model=="nonlinear"){
for(j in 1:(2^n)){
tmp <- j - 1
v <- numeric(n)  # create a vector of length n (entry-1 if species present in state)
for(k in 1:n){
if((tmp %% 2) == 0){
v[k] = 0
tmp <- floor(tmp/2)
} else {
v[k] = 1
tmp <- floor(tmp/2)
}
}
# Fraction resource lost; works for binary and flow
frac_loss <- 1 - (M[M[, i]!=0, i] %*% v / (M[M[, i]!=0, i] %*% rep(1, n)))
if(model=="linear"){
if (use_sf==1){X[j, 1] <- signif(GetLinearResponse(frac_loss, PiVector[i]), digits=significant_figures)}
else{X[j, 1] <- GetLinearResponse(frac_loss, PiVector[i])}
}
if(model=="nonlinear"){
if(use_sf==1){X[j, 1] <- signif(GetNonLinearResponse(frac_loss, PiVector[i], alpha, beta), digits=significant_figures)}
else{X[j, 1] <- GetNonLinearResponse(frac_loss, PiVector[i], alpha, beta)}
}
}
}
else if(model=="topo"){
for(j in 1:(2^n)){
tmp <- j-1
v <- numeric(n)
X[j, 1] <- 1
for(k in 1:n){
if(tmp%%2==0){
v[k] = 0
tmp <- floor(tmp/2)
} else{
v[k] = 1
X[j, 1] <- PiVector[i]
break
}
}
}
}
else if(model == "binaryplusminus_linear"){
posEdges <- sum(M[,i] > 0)  # Total number of facilitatory dependencies for i
negEdges <- sum(M[,i] < 0)  # Total number of inhibitory dependencies for i
maxunits <- max(posEdges, negEdges)
for(j in 1:(2^n)){
tmp <- j - 1
v <- numeric(n)  # create a vector of length n
## Determine which resources to consider for case j
for(k in 1:n){
if((tmp %% 2) == 0){
v[k] = 1
tmp <- floor(tmp/2)
} else {
v[k] = 0
tmp <- floor(tmp/2)
}
}
netcontributingunits <- M[M[, i] != 0, i]%*%v  # Total number of facilitatory minus inhibitory dependencies considered in case j
if(use_sf==1){
X[j, 1] <- signif(GetBinaryLinearFacInhResponse(netcontributingunits, maxunits, PiVector[i]), digits=significant_figures)
}
else{
X[j, 1] <- GetBinaryLinearFacInhResponse(netcontributingunits, maxunits, PiVector[i])
}
}
}
else if(model == "binaryplusminus_OR"){
posEdges <- sum(M[,i] > 0)  # Total number of facilitatory dependencies for i
negEdges <- sum(M[,i] < 0)  # Total number of inhibitory dependencies for i
maxunits <- max(posEdges, negEdges)
for(j in 1:(2^n)){
tmp <- j - 1
v <- numeric(n)  # create a vector of length n
## Determine which influencing species to consider for case j
for(k in 1:n){
if((tmp %% 2) == 0){
v[k] = 1
tmp <- floor(tmp/2)
} else {
v[k] = 0
tmp <- floor(tmp/2)
}
}
netcontributingunits <- M[M[, i] != 0, i]%*%v  # Total number of facilitatory minus dependencies considered in case j
if(use_sf==1){
X[j, 1] <- signif(GetBinaryORFacInhResponse(netcontributingunits, maxunits, PiVector[i]), digits=significant_figures)
}
else{
X[j, 1] <- GetBinaryORFacInhResponse(netcontributingunits, maxunits, PiVector[i])
}
}
}
else if(model == "binaryplusminus_nonlinear"){
posEdges <- sum(M[,i] > 0)  # Total number of facilitatory dependencies for i
negEdges <- sum(M[,i] < 0)  # Total number of inhibitory dependencies for i
maxunits <- max(posEdges, negEdges)
for(j in 1:(2^n)){
tmp <- j - 1
v <- numeric(n)  # create a vector of length n
## Determine which dependencies to consider for case j
for(k in 1:n){
if((tmp %% 2) == 0){
v[k] = 1
tmp <- floor(tmp/2)
} else {
v[k] = 0
tmp <- floor(tmp/2)
}
}
netcontributingunits <- M[M[, i] != 0, i]%*%v  # Total number of facilitatory minus inhibitory dependencies considered in case j
tempvector <- M[M[, i] != 0, i]
tempvector[tempvector < 1] <- 0
posEdges_state <- tempvector%*%v
tempvector <- M[M[, i] != 0, i]
tempvector[tempvector > -1] <- 0
negEdges_state <- abs(tempvector%*%v)
if(use_sf==1){
X[j, 1] <- signif(GetBinaryNonLinearFacInhResponse(netcontributingunits, posEdges_state, negEdges_state, maxunits, PiVector[i], alpha, beta), digits=significant_figures)
}
else{
X[j, 1] <- GetBinaryNonLinearFacInhResponse(netcontributingunits, posEdges_state, negEdges_state, maxunits, PiVector[i], alpha, beta)
}
}
}
else if(model == "ticks"){
posEdges <- sum(M[,i])  # Total number of facilitatory dependencies for i
negEdges <- 0  # Total number of inhibitory dependencies for i
maxunits <- max(posEdges, negEdges)
for(j in 1:(2^n)){
tmp <- j - 1
v <- numeric(n)  # create a vector of length n
## Determine which dependencies to consider for case j
for(k in 1:n){
if((tmp %% 2) == 0){
v[k] = 1
tmp <- floor(tmp/2)
} else {
v[k] = 0
tmp <- floor(tmp/2)
}
}
netcontributingunits <- M[M[, i] != 0, i]%*%v  # Total number of facilitatory minus inhibitory dependencies considered in case j
tempvector <- M[M[, i] != 0, i]
tempvector[tempvector < 1] <- 0
posEdges_state <- tempvector%*%v
tempvector <- M[M[, i] != 0, i]
tempvector[tempvector > -1] <- 0
negEdges_state <- abs(tempvector%*%v)
if(use_sf==1){
X[j, 1] <- signif(TicksResponse(netcontributingunits, posEdges_state, negEdges_state, maxunits, PiVector[i], offset[i], alpha, beta), digits=significant_figures)
}
else{
X[j, 1] <- TicksResponse(netcontributingunits, posEdges_state, negEdges_state, maxunits, PiVector[i], offset[i], alpha, beta)
}
}
}
else {
stop("You must specify a response form: binaryplusminus_linear")
}
Table[[variable_name]] <- X
}
}
Table$model <- model
return(Table)
}
############## Run gRain AND GET MARGINAL PROBS
getPreys <- function(v){
# get string of names of dependencies of a species with "+" in between
n <- length(v)
p <- sum(v > 0)
prey_str <- ''
for(i in 1:n){
if(p > 1){
if(v[i] > 0){
prey_str <- paste(prey_str, 'V', i, "+", sep='')
p <- p - 1
}
} else{
break
}
}
for(j in i:n){
if(v[j] > 0){
prey_str <- paste(prey_str, 'V', j, sep='')
}
}
return(prey_str)
}
getMarginals_gRain <- function(Table, file="test.R"){
# automatically produce an R script file for running gRain
write("require(gRain)", file=file)
write("temp_fn <- function(){", file, append=TRUE)
sp_levels <- "levels=c('absent', 'present'))"
S <- length(Table$PiVector)
for(i in 1:S){
sp_name <- paste("V", i, sep='')
if(dim(Table[[sp_name]])[1]==1){  # basal species
X <- Table[[sp_name]]
command <- paste("V", i, "<-cptable(~V", i, ", values=c(", 1 - X[1,1], ",", X[1,1], "), ", sp_levels, sep='')
write(command, file=file, append=TRUE)
} else {  # consumers
X <- Table[[sp_name]]
nX <- dim(X)[1]
command <- paste("V", i, "<-cptable(~V", i, "|", sep='')
prey_str <- getPreys(abs(Table$M[, i]))
command <- paste(command, prey_str, ", values=c(", sep='')
for(j in 1:(nX - 1)){
command <- paste(command, 1 - X[nX - j + 1, 1], ",", X[nX - j + 1, 1], ",", sep='')  # l-to-r: all prey present to all extinct
}
command <- paste(command, 1 - X[1, 1], ",", X[1, 1], "), ", sp_levels, sep='')
write(command, file, append=TRUE)
}
}
command2 <- paste("plist<-compileCPT(list(", sep='')
for(i in 1:(S-1)){
command2 <- paste(command2, "V", i, ",", sep='')
}
command2 <- paste(command2, "V", S, "))", sep='')
write(command2, file, append=TRUE)
write("BN<-grain(plist)", file, append=TRUE);
write("Results_tmp<-querygrain(BN)", file, append=TRUE)
write("return(Results_tmp)", file, append=TRUE)
write("}", file, append=TRUE)
# Run the script
source(file)
r <- temp_fn()
results <- list()
# Re-orgainse the sequence in the output
for(i in 1:S){
sp_name <- paste("V", i, sep='')
results[[sp_name]] <- r[[sp_name]]
}
return(results)
}
############## Example of using code
#M <- as.matrix(read.csv("blakeDAG_CD_BN_mod.csv", header=TRUE,row.names = 1, sep=","))Â # read BN
#M[is.na(M)] <- 0
# Tick toy
M <- matrix(c(rep(0,42),1,1,1,1,1,5,0), nrow=7, ncol=7)
Pi <- rep(0.7, dim(M)[1])
#Pi <- c(0,0,0,0,0,0,0.7)
offset <- rep(0.1, dim(M)[1])
mytable <- BuildTable(M, Pi, offset, model="ticks", alpha=1, beta=1)
# look at print(mytable) to see how changing alpha and beta affect probabilities
myresults <- getMarginals_gRain(mytable,"tempfile.R")
Posterior <- unlist(myresults)[seq(2, 2 * dim(M)[1], 2)]
r1 <- raster(matrix(nrow=3,ncol=3,c(rnorm(9,0.7,0.05))))
r2 <- raster(matrix(nrow=3,ncol=3,c(rnorm(9,0.7,0.05))))
r3 <- raster(matrix(nrow=3,ncol=3,c(rnorm(9,0.7,0.05))))
r4 <- raster(matrix(nrow=3,ncol=3,c(rnorm(9,0.7,0.05))))
r5 <- raster(matrix(nrow=3,ncol=3,c(rnorm(9,0.7,0.05))))
r6 <- raster(matrix(nrow=3,ncol=3,c(rnorm(9,0.7,0.05))))
r7 <- raster(matrix(nrow=3,ncol=3,c(rnorm(9,0.1,0.05))))
rstack <- stack(r1,r2,r3,r4,r5,r6,r7)
rastery <- function(stack.prior, M) {
stack.post <- stack.prior
for(i in 1:nrow(stack.prior)) {
for(j in 1:nrow(stack.prior)) {
Pi <- as.vector(rstack[i,j])
#Pi <- c(0,0,0,0,0,0,0.7)
offset <- rep(0.1, dim(M)[1])
mytable <- BuildTable(M, Pi, offset, model="ticks", alpha=1, beta=1)
# look at print(mytable) to see how changing alpha and beta affect probabilities
myresults <- getMarginals_gRain(mytable,"tempfile.R")
Posterior <- unlist(myresults)[seq(2, 2 * dim(M)[1], 2)]
stack.post[i,j] <- Posterior
}
print(i/nrow(stack.prior)*100)
}
return(stack.post)
}
plot(rastery(rstack,M))
BuildTable(M, Pi, offset, model="ticks", alpha=1, beta=1)
mytable <- BuildTable(M, Pi, offset, model="ticks", alpha=1, beta=1)
getMarginals_gRain(mytable,"tempfile.R")
unlist(myresults)[seq(2, 2 * dim(M)[1], 2)]
dim(myresults)
length(myresults)
library(raster)
library(rgdal)
REGIONS <- readOGR(layer="GBD-Regions-2018",dsn='C:/Users/cjcar/Dropbox/ZikaRecoil2018/Grids/Raw layers/GBD-Regions-2018')
rcp85 <- raster(file.choose())
rcp85.clip <- clip(rcp85, REGIONS)
?clip
raster::overlay(rcp85, REGIONS)
raster::overlay(REGIONS, rcp85)
fasterize(rcp85, REGIONS)
clip <- extract(rcp85, REGIONS)
library(fasterize)
clip <- fasterize(rcp85, REGIONS)
clip <- fasterize(as.sf(REGIONS), rcp85)
clip <- fasterize(REGIONS, rcp85)
messina <- raster(file.choose())
messina <- messina - messina
plot(messina)
rcp85 <- rcp85 + messina
plot(rcp85)
writeRaster(rcp85, 'C:/Users/cjcar/Dropbox/ZikaRecoil2018/Grids/rcp85clip')
install.packages('codependent')
setwd('C:/Users/cjcar/Documents/Github/codependent')
library(devtools)
virus <- read.csv("~/Github/brevity/olival nature 2017/associations.csv")
virus <- virus[,c(1:2)]
names(virus) <- c('Host','Parasite')
copredict(5128,virus,20,1)
virus <- read.csv("~/Github/brevity/olival nature 2017/associations.csv")
virus <- virus[,c(1:2)]
names(virus) <- c('Host','Parasite')
copredict(5128,virus,20,1)
14.93607*892
14.93607*894
14.93607*897
library(codependent)
data(sch2010)
library(codependent)
library(knitr)
devtools::install_github("cjcarlson/codependent")
library(spatExtinct)
data("ludo.map")
ludo.map
plot(ludo.map)
?spatExtinct
plot(spat.OLE(ludo,ludo.map))
??nls
binera.50 <-  function(assoc.df, iter, plots=FALSE) {
n.host <- n.par <- host <- pred <- 0
c <- curve.50(assoc.df, iter)
model1 <- stats::nls(n.par~b*n.host^z,start = list(b = 1, z = 0.5),data=c)
pred.df <- data.frame(pred = predict(model1), host = c$n.host)
if(plots==TRUE) {
g <- ggplot2::ggplot(c, aes(n.host, n.par)) + xlim(0,max(c$n.host)*1.05) + ylim(0,max(c$n.par)*1.05) + xlab('Hosts') + ylab('Affiliates') +
geom_point(shape = 16, size = 2.5, show.legend = FALSE, alpha = .15, color='darkturquoise') + theme_bw() +
geom_line(color='black',lwd=1,data = pred.df, aes(x=host, y=pred))
print(g)
}
return(model1)
}
setwd('C:/Users/cjcar/Documents/Github/codependent')
devtools::check()
devtools::check()
load_all
load_all()
library(devtools)
load_all()
codependent::binera()
?codependent::binera()
codependent::binera.50()
?codependent
??codependent
binera(rob1929,plots=TRUE)
binera(rob1929,iter=1,plots=TRUE)
binera(rob1929,iter=5,plots=TRUE)
devtools::check()
devtools::use_package("stats", "Import")
devtools::use_package("graphics", "Import")
document_all()
load_all()
load_all()
devtools::check()
qt
??qt
?qt
importFrom(stats, nls)
load_all()
devtools::check()
?devtools::use_package
#' @importFrom stats nls
load_all()
load_all()
devtools::check()
document()
document()
document()
#' @importFrom(stats,nls)
document()
document()
document()
document()
document()
document9)
document()
document()
document()
document()
devtools::check()
document()
copredict(100, rob1929, 10, 10)
??codependent
?codependent
document()