.Rhistory

# vector of T
t = seq(from=5, to = 40, by=0.5)
# survival
p = exp(-1/(-4.4+1.31*t-0.03*t^2))
plot(t,p)
# MEI
a = 0.5
mei = (a^2*p^(111/(t-16)))/(-log(p))
# plot
plot(t,mei,type='l')
### other metric in McCord paper for sporogony: r
num = 0.06044*(t/296.65)*exp((17545/1.987)*(1/296.65 - 1/t))
den = 1 + exp((-142843/1.987)*(1/288.85 - 1/t)) + exp((110980/1.987)*(1/306.90 - 1/t))
r = num/den
plot(t,r) # hmmmmm can't be right
mei_alt = (a^2*p^r)/(-log(p))
plot(t,mei_alt)
library(raster)
alt <- raster('C:/Users/cjcar/Dropbox/Cbot_global2019/GIS/scaled/alt_resample.asc')
load('C:/Users/cjcar/Downloads/NEGridList.Rdata')
blank <- raster('D:/ICEBERG/UniversalBlank.tif')
values(blank) <- NEGridList$Futures1$OverlapSum
alt <- resample(alt, blank)
d <- data.frame(alt = alt@data@values,
new = blank@data@values)
p <- ggplot(d, aes(alt, new)) + geom_smooth(method = "gam")
library(ggplot2)
p <- ggplot(d, aes(alt, new)) + geom_smooth(method = "gam")
p
plot(alt)
plot(blank)
p <- ggplot(d, aes(alt, new)) + geom_smooth(method = "gam") + geom_rug()
p
p <- ggplot(d, aes(alt, new)) + geom_smooth(method = "gam") + geom_point()
p <- ggplot(d, aes(alt, new)) + geom_smooth(method = "gam") + geom_point()
p <- ggplot(d, aes(alt, new)) + geom_hex()
p
p <- ggplot(d, aes(alt, new)) + geom_smooth(method = "auto")
p
writeRaster(alt,'altitude.tif')
plot(alt)
plot(alt[alt==NA])
iron <- raster('C:/Users/cjcar/Dropbox/PlagueWNA/Layers/Top5_Fe.tif')
library(raster)
iron <- raster('C:/Users/cjcar/Dropbox/PlagueWNA/Layers/Top5_Fe.tif')
plot(iron)
iron
r <- readOGR('C:/Users/cjcar/Downloads',
'Fe.kmz')
library(rgdal)
r <- readOGR('C:/Users/cjcar/Downloads',
'Fe.kmz')
r <- readOGR(dsn='C:/Users/cjcar/Downloads',
layer='Fe.kmz')
list.files(C:/Users/cjcar/Downloads)
list.files('C:/Users/cjcar/Downloads')
r <- readOGR(dsn='C:/Users/cjcar/Downloads',
layer='Fe.kmz', 'Fe Top 0- to 5-cm layer')
library(raster)
library(rgdal)
global.par <- readOGR(dsn='C:/Users/cjcar/Dropbox/continents/gadm28_adm0',
layer='gadm28_adm0')
# Chunk 1
knitr::opts_chunk$set(
collapse = TRUE,
comment = "#>"
)
suppressMessages(library(tidyverse))
# Chunk 2: setup
library(embarcadero)
library(velox)
set.seed(12345)
# Chunk 3
data(covs)
cov.big <- bigstack(covs, 10)
# Chunk 4
data(ticks)
head(ticks)
nrow(ticks)
# Chunk 5
mod <- SpatialPointsDataFrame(ticks[,3:4],data.frame(ticks[,1]))
names(mod@data) <- 'Presence'
# Rasterizing makes unique points to the grid cell
tmp=rasterize(mod, covs[[1]], field="Presence", fun="min")
pts.sp1=rasterToPoints(tmp, fun=function(x){x>0})
nrow(pts.sp1)
pres.cov <- raster::extract(covs, pts.sp1[,1:2])
head(pres.cov)
# Chunk 6
#Generate the data
absence <- randomPoints(covs,nrow(ticks))
abs.cov <- raster::extract(covs, absence)
#Code the response
pres.cov <- data.frame(pres.cov); pres.cov$tick <- 1
abs.cov <- data.frame(abs.cov); abs.cov$tick <- 0
# And one to bind them
all.cov <- rbind(pres.cov, abs.cov)
head(all.cov)
# Let's just clean it up a little bit
all.cov <- all.cov[complete.cases(all.cov),]
# Chunk 7
first.model <- bart(all.cov[,1:11], all.cov[,'tick'], keeptrees=TRUE)
first.model$fit
first.model$fit$sampleTreesFromPrior
first.model$fit$printTrees
first.model$fit$printTrees()
first.model$varcount
nrow(first.model$varcount)
??predict
library(dbarts)
??predict
?dbarts
library(rgdal)
load('C:/Users/cjcar/Downloads/gadm36_IND_2_sp.rds')
loadRDS('C:/Users/cjcar/Downloads/gadm36_IND_2_sp.rds')
??RDS
load("C:/Users/cjcar/Downloads/gadm36_IND_2_sp.rds")
load("C:/Users/cjcar/Downloads/gadm36_IND_2_sf.rds")
library(tidyverse)
library(rgdal)
gbd <- readOGR('C:/Users/cjcar/Dropbox/ZikaRecoil2018/Grids/Raw layers/GBD-Regions-2018/GBD-Regions-2018.shp')
ssp5 <- read.csv('C:/Users/cjcar/Dropbox/ZikaRecoil2018/Analysis/Any month risk_SSP5.csv')
ssp5 %>% filter(rcp == 'rcp85') %>%
group_by(Region) %>% summarize(First = mean(Firsts)) %>%
right_join(gbd@data, by='Region') -> gbd@data
gbd %>% st_as_sf -> gbdsf
plot(st_geometry(gbdsf), col="darkseagreen3", border="darkseagreen4",
bg = "lightblue1", lwd = 0.5)
library(sf)
gbd <- readOGR('C:/Users/cjcar/Dropbox/ZikaRecoil2018/Grids/Raw layers/GBD-Regions-2018/GBD-Regions-2018.shp')
ssp5 <- read.csv('C:/Users/cjcar/Dropbox/ZikaRecoil2018/Analysis/Any month risk_SSP5.csv')
ssp5 %>% filter(rcp == 'rcp85') %>%
group_by(Region) %>% summarize(First = mean(Firsts)) %>%
right_join(gbd@data, by='Region') -> gbd@data
gbd %>% st_as_sf -> gbdsf
plot(st_geometry(gbdsf), col="darkseagreen3", border="darkseagreen4",
bg = "lightblue1", lwd = 0.5)
devtools::install_github('cjcarlson/embarcadero')
install.packages('ggpubr')
devtools::install_github('cjcarlson/embarcadero')
install.packages("installr")
library(installr)
updateR()
library(rworldmap)
install.packages('rworldmap')
library(rworldmap)
data(countryExData)
countryExData
View(countryExData)
sPDF <- joinCountryData2Map( countryExData
,joinCode = "ISO3"
,nameJoinColumn = "ISO3V10")
mapDevice() #create world map shaped window
mapCountryData(sPDF
,nameColumnToPlot='BIODIVERSITY')
## 'Not run' code you need rjags (and JAGS) to be installed.
require(mgcv)
set.seed(2) ## simulate some data...
n <- 400
dat <- gamSim(1,n=n,dist="normal",scale=2)
## regular gam fit for comparison...
b0 <- gam(y~s(x0)+s(x1) + s(x2)+s(x3),data=dat,method="REML")
## Set directory and file name for file containing jags code.
## In real use you would *never* use tempdir() for this. It is
## only done here to keep CRAN happy, and avoid any chance of
## an accidental overwrite. Instead you would use
## setwd() to set an appropriate working directory in which
## to write the file, and just set the file name to what you
## want to call it (e.g. "test.jags" here).
jags.file <- paste(tempdir(),"/test.jags",sep="")
## Set up JAGS code and data. In this one might want to diagonalize
## to use conjugate samplers. Usually call 'setwd' first, to set
## directory in which model file ("test.jags") will be written.
jd <- jagam(y~s(x0)+s(x1)+s(x2)+s(x3),data=dat,file=jags.file,
sp.prior="gamma",diagonalize=TRUE)
plot(jd)
jd
summary(jd)
library(BBS.SDM)
?BBS.SDM
View(BBS.SDM::RunMod())
View(BBS.SDM::RunMod)
?jagam
vignette(BBS.SDM)
vignette("BBS.SDM")
vignette("bbs-sdm")
library(BBS.tenstop)
bbs <- get_BBS10()
target_dir <- here::here("C:/Users/cjcar/Desktop/hole 2")
GetCorrData(bbs_raw = bbs, alpha = "FICR", path = target_dir,
start.year = 1972, end.year = 2014)
setwd("C:/Users/cjcar/Desktop/hole 2")
target_dir <- here::here(".")
GetCorrData(bbs_raw = bbs, alpha = "FICR", path = target_dir,
start.year = 1972, end.year = 2014)
bba
bbs
View(bbs)
View(bbs$counts)
GetBioVars("FICR", path = target_dir)
GetInits(alpha = "FICR", path = target_dir)
target_dir
target_dir <- here::here(".")
target_dir
target_dir <- getwd('')
target_dir <- getwd()
GetInits(alpha = "FICR", path = target_dir)
GetCorrData(bbs_raw = bbs, alpha = "FICR", path = target_dir,
start.year = 1972, end.year = 2014)
GetBioVars("FICR", path = target_dir)
GetInits(alpha = "FICR", path = target_dir)
target_dir
target_dir <- here::here("output") # change 'output' to the desired directory (relative to the current WD)
target_dir
getwd()
target_dir <- here::here("output") # change 'output' to the desired directory (relative to the current WD)
setwd('C:/Users/cjcar/Desktop/r hole')
GetCorrData(bbs_raw = bbs, alpha = "FICR", path = target_dir,
start.year = 1972, end.year = 2014)
data <- read.csv('C:/Users/cjcar/Documents/GitHub/humanhelminths/Lit review data/cleaned lit review results.csv')
data$valid_name
?paste
paste("(",paste(data,collapse='\" or \"'),sep='')
paste("(",paste(data$valid_name,collapse='\" or \"'),sep='')
paste("(",paste(data$valid_name,collapse='\\" or \\"'),sep='')
paste("(",paste(data$valid_name,collapse='" or "'),sep='')
paste(paste("(",paste(data$valid_name,collapse='" or "'),sep=''),'")',sep='')
write.csv(x)
x <- paste(paste("(",paste(data$valid_name,collapse='" or "'),sep=''),'")',sep='')
write.csv(x)
write.csv(x,'test.csv')
library(raster)
cur <- raster('file:///C:/Users/cjcar/Dropbox/Global Algae Macroecology & Mapping Analysis (GAMMA)/OUTPUTS/allmaps/acanthococcus antarcticus current.tif')
cur <- raster('C:/Users/cjcar/Dropbox/Global Algae Macroecology & Mapping Analysis (GAMMA)/OUTPUTS/allmaps/acanthococcus antarcticus current.tif')
plot(cur)
cur
(1-0.23)*3234
(1-0.15)*37904
(1-0.07)*7250
(1-0.06)*19628
(1-0.16)*24144
(1-0.26)*17399
r0t <- function(T) {
a = 0.000203*T*(T-11.7)*((42.3-T)^0.5)
bc = -0.54*T*T + 25.2*T - 206
p = -0.000828*T*T + 0.0367*T + 0.522 # e^-mu
mu = -1*log(p)
PDR = 0.000111*T*(T-14.7)*((34.4-T)^0.5) # 1/EIP
pEA = -0.00924*T*T + 0.453*T - 4.77
MDR = 0.000111*T*(T-14.7)*((34-T)^0.5) # 1/tauEA
EFD = -0.153*T*T + 8.61*T - 97.7
R0 = (((a^2)*bc*(p^(1/PDR))*EFD*pEA*MDR)/(mu^3))^(1/2)
return(R0)
}
lapply(c(15:35), r0t)
temps <- lapply(c(15:35), r0t)
temps <- 15:35
r0 <- lapply(c(15:35), r0t)
r0 <- unlist(ro)
r0 <- unlist(r0)
plot(r0 ~ temps)
temps <- c(15:35, by=0.2)
temps
temps <- seq(15:35, by=0.2)
temps <- seq(15, 35, by=0.2)
temps
r0 <- sapply(c(15:35), r0t)
plot(r0 ~ temps)
r0
temps <- seq(15, 34, by=0.2)
r0 <- sapply(temps, r0t)
plot(r0 ~ temps)
temps <- raster('C:/Users/cjcar/Documents/GitHub/embarcadero/vignettes/covariates/bio1.asc')
library(raster)
temps <- raster('C:/Users/cjcar/Documents/GitHub/embarcadero/vignettes/covariates/bio1.asc')
plot(temps)
plot(r0t(temps))
r0t <- function(T) {
a = 0.000203*T*(T-11.7)*((42.3-T)^0.5)
bc = -0.54*T*T + 25.2*T - 206
p = -0.000828*T*T + 0.0367*T + 0.522 # e^-mu
mu = -1*log(p)
PDR = 0.000111*T*(T-14.7)*((34.4-T)^0.5) # 1/EIP
pEA = -0.00924*T*T + 0.453*T - 4.77
MDR = 0.000111*T*(T-14.7)*((34-T)^0.5) # 1/tauEA
EFD = -0.153*T*T + 8.61*T - 97.7
R0 = (((a^2)*bc*(p^(1/PDR))*EFD*pEA*MDR)/(mu^3))^(1/2)
if(is.nan(R0)){return(0)}
return(R0)
}
temps <- seq(14, 36, by=0.2)
r0 <- sapply(temps, r0t)
plot(r0 ~ temps)
temps <- raster('C:/Users/cjcar/Documents/GitHub/embarcadero/vignettes/covariates/bio1.asc')
plot(temps)
plot(r0t(temps))
r0t <- function(T) {
if(is.na(T)) return(NA)
a = 0.000203*T*(T-11.7)*((42.3-T)^0.5)
bc = -0.54*T*T + 25.2*T - 206
p = -0.000828*T*T + 0.0367*T + 0.522 # e^-mu
mu = -1*log(p)
PDR = 0.000111*T*(T-14.7)*((34.4-T)^0.5) # 1/EIP
pEA = -0.00924*T*T + 0.453*T - 4.77
MDR = 0.000111*T*(T-14.7)*((34-T)^0.5) # 1/tauEA
EFD = -0.153*T*T + 8.61*T - 97.7
R0 = (((a^2)*bc*(p^(1/PDR))*EFD*pEA*MDR)/(mu^3))^(1/2)
if(is.nan(R0)){return(0)}
return(R0)
}
r0T(NA)
r0t(NA)
r0t <- function(T) {
if(is.na(T)) return(NA)
a = 0.000203*T*(T-11.7)*((42.3-T)^0.5)
bc = -0.54*T*T + 25.2*T - 206
p = -0.000828*T*T + 0.0367*T + 0.522 # e^-mu
mu = -1*log(p)
PDR = 0.000111*T*(T-14.7)*((34.4-T)^0.5) # 1/EIP
pEA = -0.00924*T*T + 0.453*T - 4.77
MDR = 0.000111*T*(T-14.7)*((34-T)^0.5) # 1/tauEA
EFD = -0.153*T*T + 8.61*T - 97.7
R0 = (((a^2)*bc*(p^(1/PDR))*EFD*pEA*MDR)/(mu^3))^(1/2)
if(is.nan(R0)){return(0)}
return(R0)
}
plot(r0t(temps))
plot(c <- calc(temps, r0t))
max(c)
calc(c, max)
max(values(c))
max(values(c), na.omit=TRUE)
max(na.omit(values(c)))
plot(c/max(na.omit(values(c))))
plot(temps)
r = raster("tmean")
r = raster("tmean01")
w = getData('worldclim', var='tmean', res=0.5, lon=5, lat=45)
plot(w)
w = getData('worldclim', var='tmean', res=0.5)
?getData
list.files('C:/Users/cjcar/Downloads/wc2.0_2.5m_tavg', pattern='tif')
s <- stack(sapply(list.files('C:/Users/cjcar/Downloads/wc2.0_2.5m_tavg', pattern='tif'), raster))
s <- do.call(lapply(list.files('C:/Users/cjcar/Downloads/wc2.0_2.5m_tavg', pattern='tif'), raster), stack)
lapply(list.files('C:/Users/cjcar/Downloads/wc2.0_2.5m_tavg', pattern='tif'), raster)
s <- do.call(lapply(list.files('C:/Users/cjcar/Downloads/wc2.0_2.5m_tavg', pattern='tif', full.names = TRUE), raster), stack)
lapply(list.files('C:/Users/cjcar/Downloads/wc2.0_2.5m_tavg', pattern='tif', full.names = TRUE)
lapply(list.files('C:/Users/cjcar/Downloads/wc2.0_2.5m_tavg', pattern='tif', full.names = TRUE), raster)
s <- do.call(stack, lapply(list.files('C:/Users/cjcar/Downloads/wc2.0_2.5m_tavg', pattern='tif', full.names = TRUE), raster)
)
plot(s)
t <- mean(s)
plot(s[[1]])
plot(calc(s[[1]], r0t))
-1 <= 0
r0t <- function(T) {
if(T<=14) {return(0)}
if(is.na(T)) return(NA)
a = 0.000203*T*(T-11.7)*((42.3-T)^0.5)
bc = -0.54*T*T + 25.2*T - 206
p = -0.000828*T*T + 0.0367*T + 0.522 # e^-mu
mu = -1*log(p)
PDR = 0.000111*T*(T-14.7)*((34.4-T)^0.5) # 1/EIP
pEA = -0.00924*T*T + 0.453*T - 4.77
MDR = 0.000111*T*(T-14.7)*((34-T)^0.5) # 1/tauEA
EFD = -0.153*T*T + 8.61*T - 97.7
R0 = (((a^2)*bc*(p^(1/PDR))*EFD*pEA*MDR)/(mu^3))^(1/2)
if(is.nan(R0)){return(0)}
return(R0)
}
plot(calc(s[[1]], r0t))
r0t <- function(T) {
if(is.na(T)) return(NA)
if(T<=14) {return(0)}
a = 0.000203*T*(T-11.7)*((42.3-T)^0.5)
bc = -0.54*T*T + 25.2*T - 206
p = -0.000828*T*T + 0.0367*T + 0.522 # e^-mu
mu = -1*log(p)
PDR = 0.000111*T*(T-14.7)*((34.4-T)^0.5) # 1/EIP
pEA = -0.00924*T*T + 0.453*T - 4.77
MDR = 0.000111*T*(T-14.7)*((34-T)^0.5) # 1/tauEA
EFD = -0.153*T*T + 8.61*T - 97.7
R0 = (((a^2)*bc*(p^(1/PDR))*EFD*pEA*MDR)/(mu^3))^(1/2)
if(is.nan(R0)){return(0)}
return(R0)
}
plot(calc(s[[1]], r0t))
s2 <- calc(s, r0t)
s
s2 <- stackApply(s, c(1:12), r0t)
r0t <- function(T, na.rm=TRUE) {
if(is.na(T)) return(NA)
if(T<=14) {return(0)}
a = 0.000203*T*(T-11.7)*((42.3-T)^0.5)
bc = -0.54*T*T + 25.2*T - 206
p = -0.000828*T*T + 0.0367*T + 0.522 # e^-mu
mu = -1*log(p)
PDR = 0.000111*T*(T-14.7)*((34.4-T)^0.5) # 1/EIP
pEA = -0.00924*T*T + 0.453*T - 4.77
MDR = 0.000111*T*(T-14.7)*((34-T)^0.5) # 1/tauEA
EFD = -0.153*T*T + 8.61*T - 97.7
R0 = (((a^2)*bc*(p^(1/PDR))*EFD*pEA*MDR)/(mu^3))^(1/2)
if(is.nan(R0)){return(0)}
return(R0)
}
s2 <- stackApply(s, c(1:12), r0t)
plot(s2)
View(s2)
plot(calc(s2, mean))
library(available)
suggest(text = "Species distribution modeling with Bayesian additive regression trees")
suggest(text = "Species distribution modeling with BART")
suggest(text = "SDMs with BART")
suggest(text = "SDM with BART")
human <- velox('C:/Users/cjcar/Dropbox/CEPI Vaccine Stockpiling/Nipah/Code for Spillover Layer/human/humanfine.grd')
library(raster)
library(rgdal)
library(rgeos)
library(sp)
library(velox)
human <- velox('C:/Users/cjcar/Dropbox/CEPI Vaccine Stockpiling/Nipah/Code for Spillover Layer/human/humanfine.grd')
gadm <- readOGR(dsn='C:/Users/cjcar/Dropbox/CEPI Vaccine Stockpiling/MERS/Country_layers',
layer='middle_east')
human.pop <- human$extract(sp = gadm, fun = function(x) sum(x, na.rm = TRUE))
human
hostest <- function(codobj, rate) {(rate*codobj[[1]][[1]]/codobj[[2]][[1]])^(1/codobj[[2]][[2]])}
hostest(est.a.amp, 0.5)
##########
# RUN 03 - RATE OF DESCRIPTION.R FIRST
library(codependent)
library(tidyverse)
library(reshape2)
set.seed(8675309)
setwd('~/Github/helminths')
all.data <- read_delim('./data/associations cleaned v2.csv',delim=',')[,-1]
all.data$hostgroup[all.data$hostgroup=='Chondrostei'] <- 'Osteichthyes'
all.data$hostgroup[all.data$hostgroup=='Holostei'] <- 'Osteichthyes'
all.data$hostgroup[all.data$hostgroup=='Cladistei'] <- 'Osteichthyes'
all.data$hostgroup[all.data$hostgroup=='Teleostei'] <- 'Osteichthyes'
cest.df <- all.data[all.data$group=='Acanthocephalans',c('Parasite','hostgroup')]
cest.df$value =1
cest.df <- unique(na.omit(cest.df))
cest.df <- acast(cest.df , Parasite ~ hostgroup)
cest.df[is.na(cest.df)] <- 0
cest.df <- cest.df[,c('Amphibia', 'Aves','Mammalia','Reptilia','Osteichthyes','Chondrichthyes')] # ADD FISH BACK - BUT > 2 FISH GROUPS CURRENTLY
cest.df <- cest.df[rowSums(cest.df)>0,]
all.data %>% filter(group=='Acanthocephalans') %>% filter(hostgroup=='Amphibia') %>%
select(Host, Parasite) %>% data.frame() %>% unique() -> cest.amph
all.data %>% filter(group=='Acanthocephalans') %>% filter(hostgroup=='Aves') %>%
select(Host, Parasite) %>% data.frame() %>% unique() -> cest.aves
all.data %>% filter(group=='Acanthocephalans') %>% filter(hostgroup=='Mammalia') %>%
select(Host, Parasite) %>% data.frame() %>% unique() -> cest.mamm
all.data %>% filter(group=='Acanthocephalans') %>% filter(hostgroup=='Reptilia') %>%
select(Host, Parasite) %>% data.frame() %>% unique() -> cest.rept
all.data %>% filter(group=='Acanthocephalans') %>% filter(hostgroup=='Osteichthyes') %>%
select(Host, Parasite) %>% data.frame() %>% unique() -> cest.oste
all.data %>% filter(group=='Acanthocephalans') %>% filter(hostgroup=='Chondrichthyes') %>%
select(Host, Parasite) %>% data.frame() %>% unique() -> cest.chon
est.c.amp <- copredict(cest.amph, n=7302, iter=100)
est.c.ave <- copredict(cest.aves, n=10425, iter=100)
est.c.mam <- copredict(cest.mamm, n=5513, iter=100)
est.c.rep <- copredict(cest.rept, n=10038, iter=100)
est.c.ost <- copredict(cest.oste, n=28000, iter=100)
est.c.con <- copredict(cest.chon, n=1111, iter=100)
hostest(est.c.amp, 0.5) - length(unique(cest.amph$Host))
hostest(est.c.ave, 0.5) - length(unique(cest.aves$Host))
hostest(est.c.mam, 0.5) - length(unique(cest.mamm$Host))
hostest(est.c.rep, 0.5) - length(unique(cest.rept$Host))
hostest(est.c.ost, 0.5) - length(unique(cest.oste$Host))
hostest(est.c.con, 0.5) - length(unique(cest.chon$Host))
# RUN 03 - RATE OF DESCRIPTION.R FIRST
library(codependent)
library(tidyverse)
library(reshape2)
set.seed(8675309)
setwd('~/Github/helminths')
all.data <- read_delim('./data/associations cleaned v2.csv',delim=',')[,-1]
all.data$hostgroup[all.data$hostgroup=='Chondrostei'] <- 'Osteichthyes'
all.data$hostgroup[all.data$hostgroup=='Holostei'] <- 'Osteichthyes'
all.data$hostgroup[all.data$hostgroup=='Cladistei'] <- 'Osteichthyes'
all.data$hostgroup[all.data$hostgroup=='Teleostei'] <- 'Osteichthyes'
cest.df <- all.data[all.data$group=='Nematodes',c('Parasite','hostgroup')]
cest.df$value =1
cest.df <- unique(na.omit(cest.df))
cest.df <- acast(cest.df , Parasite ~ hostgroup)
cest.df[is.na(cest.df)] <- 0
cest.df <- cest.df[,c('Amphibia', 'Aves','Mammalia','Reptilia','Osteichthyes','Chondrichthyes')] # ADD FISH BACK - BUT > 2 FISH GROUPS CURRENTLY
cest.df <- cest.df[rowSums(cest.df)>0,]
all.data %>% filter(group=='Nematodes') %>% filter(hostgroup=='Amphibia') %>%
select(Host, Parasite) %>% data.frame() %>% unique() -> cest.amph
all.data %>% filter(group=='Nematodes') %>% filter(hostgroup=='Aves') %>%
select(Host, Parasite) %>% data.frame() %>% unique() -> cest.aves
all.data %>% filter(group=='Nematodes') %>% filter(hostgroup=='Mammalia') %>%
select(Host, Parasite) %>% data.frame() %>% unique() -> cest.mamm
all.data %>% filter(group=='Nematodes') %>% filter(hostgroup=='Reptilia') %>%
select(Host, Parasite) %>% data.frame() %>% unique() -> cest.rept
all.data %>% filter(group=='Nematodes') %>% filter(hostgroup=='Osteichthyes') %>%
select(Host, Parasite) %>% data.frame() %>% unique() -> cest.oste
all.data %>% filter(group=='Nematodes') %>% filter(hostgroup=='Chondrichthyes') %>%
select(Host, Parasite) %>% data.frame() %>% unique() -> cest.chon
est.c.amp <- copredict(cest.amph, n=7302, iter=100)
est.c.ave <- copredict(cest.aves, n=10425, iter=100)
est.c.mam <- copredict(cest.mamm, n=5513, iter=100)
est.c.rep <- copredict(cest.rept, n=10038, iter=100)
est.c.ost <- copredict(cest.oste, n=28000, iter=100)
est.c.con <- copredict(cest.chon, n=1111, iter=100)
hostest(est.c.amp, 0.5) - length(unique(cest.amph$Host))
hostest(est.c.ave, 0.5) - length(unique(cest.aves$Host))
hostest(est.c.mam, 0.5) - length(unique(cest.mamm$Host))
hostest(est.c.rep, 0.5) - length(unique(cest.rept$Host))
hostest(est.c.ost, 0.5) - length(unique(cest.oste$Host))
hostest(est.c.con, 0.5) - length(unique(cest.chon$Host))