biggest=co.env.spe2[top$V1[1],top$V2[1]]
smallest=co.env.spe2[top$V1[length],top$V2[length]]
# co.env.spe[which(co.env.spe>0 & co.env.spe<smallest)]=0
# co.env.spe[which(co.env.spe> -smallest & co.env.spe<0)]=0
#for extracting only minus
co.env.spe[which(co.env.spe>smallest & co.env.spe<0)]=0
co.env.spe3=abs(co.env.spe3)
co.env.spe3[which(co.env.spe3 == 1)]= 0
choose=(choose(spno,2)+spno)*0.025*2
top=arrayInd(sort.list(co.env.spe3,decreasing=T)[1:choose],dim(co.env.spe3))# top 5%
top=as.data.frame(top)
length=length(top$V1)
biggest=co.env.spe3[top$V1[1],top$V2[1]]
smallest=co.env.spe3[top$V1[length],top$V2[length]]
co.env.spe[which(co.env.spe>0 & co.env.spe<smallest)]=0
corrplot(co.env.spe,method = "ellipse",shade.col = NA, order = "alphabet",tl.col ="black",hclust.method = "centroid",tl.srt = 10,type="lower",tl.cex = 0.4,col=col2(200),number.font = 4)
x=an
type="R2_McFadden"
#x$species$R2_McFadden_shared
#modify by the original plot.anova from sjsdm package
internals = list()
df = data.frame(
env = ifelse(x$sites[[type]]$F_A<0, 0, x$sites[[type]]$F_A),
spa = ifelse(x$sites[[type]]$F_S<0, 0, x$sites[[type]]$F_S),
codist = ifelse(x$sites[[type]]$F_B<0, 0, x$sites[[type]]$F_B),
r2  = ifelse(x$sites[[type]]$Full<0, 0, x$sites[[type]]$Full)
)
internals[[1]] = df
names(internals)[1] = "Sites"
df = data.frame(
env = ifelse(x$species[[type]]$F_A<0, 0, x$species[[type]]$F_A),
spa = ifelse(x$species[[type]]$F_S<0, 0, x$species[[type]]$F_S),
codist = ifelse(x$species[[type]]$F_B<0, 0, x$species[[type]]$F_B),
r2  = ifelse(x$species[[type]]$Full<0, 0, x$species[[type]]$Full),
group=NA,
name=TP$species_name,
sp=model$species
)
df$sp<-as.character(df$sp)
if (waterspecies=="FALSE") {
df$group[grep("Acti*",df$sp)]="Fish"
df$group[grep("Amph*",df$sp)]="Amphibians"
df$group[grep("Mamm*",df$sp)]="Mammal"
df$group[grep("Aves*",df$sp)]="Perching birds"
df$group[grep("Wafl",df$sp)]="Waterfowl"
df$group[grep("Dmst",df$sp)]="Domestic animals"
df=df %>% select (-sp)
df$group=factor(df$group,levels = c('Fish','Amphibians','Perching birds','Domestic animals','Mammal','Waterfowl'))
}
if (waterspecies=="TRUE") {
df$group[grep("Acti*",df$sp)]="Fish"
df$group[grep("Amph*",df$sp)]="Amphibians"
df=df %>% select (-sp)
df$group=factor(df$group,levels = c('Amphibians','Fish'))
}
internals[[2]] = df
names(internals)[2] = "Species"
SP=data.frame(sp)
OTU=data.frame(Y)
Env=data.frame(env)
indexdata=cbind(Env,SP)
indexdata$richness=rowSums(Y)
indexdata.OTU=cbind(indexdata,OTU)
rm(SP,OTU,Env)
site.date=data.frame(internals[[1]])
indexdata.OTU2=cbind(indexdata.OTU,site.date)
#does r2 =0 is meaning for the sites?
#site.date=site.date %>% filter(r2!=0)
#indexdata.OTU2=indexdata.OTU2 %>% filter(r2!=0)
indexdata.OTU2$r2=site.date$r2
indexdata.OTU2$env=site.date$env
indexdata.OTU2$spa=site.date$spa
indexdata.OTU2$codist=site.date$codist
library(ggtern)
top = 1
r2max = ceiling(max(internals[[2]]$r2)*1e2)/1e2
internals.name=data.frame(name=internals[[2]]$name)
internals.name<-separate(internals.name,name,c("g","species"),sep= "_")
internals.name<-paste0(internals.name$g," ",internals.name$species)
internals[[2]]$name=internals.name
ggtern::ggtern(internals[[2]], ggplot2::aes(x = env, z = spa, y = codist, size = r2,color =group)) +
ggtern::scale_T_continuous(limits=c(0,1),
breaks=seq(0, 1,by=0.2),
labels=seq(0,1, by= 0.2)) +
ggtern::scale_L_continuous(limits=c(0,1),
breaks=seq(0, 1,by=0.2),
labels=seq(0, 1,by=0.2)) +
ggtern::scale_R_continuous(limits=c(0,1),
breaks=seq(0, 1,by=0.2),
labels=seq(0, 1,by=0.2)) +
ggplot2::labs(title = names(internals)[2],
x = "E",
xarrow = "Environment",
y = "C",
yarrow = "Co-Distribution",
z = "S",
zarrow = "Space") +
ggtern::theme_bw() +
ggtern::theme_showarrows() +
ggtern::theme_arrowlong() +
ggplot2::theme(
panel.grid = ggplot2::element_line(color = "darkgrey", linewidth = 0.3),
plot.tag = ggplot2::element_text(size = 11),
plot.title = ggplot2::element_text(size = 11, hjust = 0.1 , margin = ggplot2::margin(t = 10, b = -20)),
tern.axis.arrow = ggplot2::element_line(linewidth = 1),
tern.axis.arrow.text = ggplot2::element_text(size = 6),
axis.text = ggplot2::element_text(size = 4),
axis.title = ggplot2::element_text(size = 6),
legend.text = ggplot2::element_text(size = 6),
legend.title = ggplot2::element_text(size = 8),
strip.text = ggplot2::element_text(size = 8),
plot.margin = unit(c(top,1,1,1)*0.2, "cm"),
strip.background = ggplot2::element_rect(color = NA),
) +
ggplot2::theme(tern.axis.arrow.text = element_text(size = 7),legend.position = "bottom", legend.margin = margin(r = 30), legend.box="vertical")+
ggplot2::guides(size = ggplot2::guide_legend(title = expression(R^2), order = 1, nrow = 1, label.position = "bottom"), color = ggplot2::guide_legend(title = "Species grouping", title.position = "top", order = 2, barheight = 0.5, barwidth = 8))+
scale_color_viridis(option="D",discrete=TRUE)+
geom_point(alpha = 0.7,position= position_jitter_tern(x=0.01, y=0.01, z=0.00))+
geom_text(label=internals[[2]]$name,alpha = 0.6,size=1,col="black",check_overlap = F,show.legend = F,fontface = "italic")+ggtern::theme_nomask()
plots_internals = list()
#show the environmental deviance instead of species richness
Zscores = apply(indexdata.OTU2[,c(1:8)], 2, function(x) (x-mean(x))/sd(x))
env_deviance =rowSums(Zscores)
env_deviance_name="env_deviance"
#show the environmental deviance instead of species richness
#Zscores = apply(indexdata.OTU2[,c(1:8)], 2, function(x) (x-mean(x))/sd(x))
#env_deviance =rowSums(Zscores)
env_deviance_name="env_deviance"
#env_deviance =indexdata.OTU2$richness
Zscores = data.frame(indexdata.OTU2[,c(1:5)])
Zscores$HSI2_Pond_area[which(Zscores$HSI2_Pond_area == "0")]=0.01
Zscores$HSI2_Pond_area=log(Zscores$HSI2_Pond_area)
plots_internals = list()
#show the environmental deviance instead of species richness
Zscores = apply(indexdata.OTU2[,c(1:8)], 2, function(x) (x-mean(x))/sd(x))
env_deviance =rowSums(Zscores)
env_deviance_name="env_deviance"
# Zscores$HSI2_Pond_area[which(Zscores$HSI2_Pond_area == "0")]=0.01
# Zscores$HSI2_Pond_area=log(Zscores$HSI2_Pond_area)
# hist(Zscores$HSI2_Pond_area)
# Zscores=as.matrix(Zscores)
# Zscores=scales::rescale(Zscores, to=c(-1,1))
# hist(as.data.frame(Zscores)$HSI2_Pond_area)
# Zscores=abs(Zscores)
# env_deviance =rowSums(Zscores)
# hist(env_deviance)
############options for sites ternary#########
top = 7
r2max = ceiling(max(site.date$r2)*1e2)/1e2
ggtern::ggtern(site.date, ggplot2::aes(x = env, z = spa, y = codist, size = r2,color =env_deviance)) +
ggtern::scale_T_continuous(limits=c(0,1),
breaks=seq(0, 1,by=0.2),
labels=seq(0,1, by= 0.2)) +
ggtern::scale_L_continuous(limits=c(0,1),
breaks=seq(0, 1,by=0.2),
labels=seq(0, 1,by=0.2)) +
ggtern::scale_R_continuous(limits=c(0,1),
breaks=seq(0, 1,by=0.2),
labels=seq(0, 1,by=0.2)) +
ggplot2::labs(title = names(internals)[1],
x = "E",
xarrow = "Environment",
y = "C",
yarrow = "Co-Distribution",
z = "S",
zarrow = "Space") +
ggtern::theme_bw() +
ggtern::theme_showarrows() +
ggtern::theme_arrowlong() +
ggplot2::theme(
panel.grid = ggplot2::element_line(color = "darkgrey", linewidth = 0.3),
plot.tag = ggplot2::element_text(size = 11),
plot.title = ggplot2::element_text(size = 11, hjust = 0.1 , margin = ggplot2::margin(t = 10, b = -20)),
tern.axis.arrow = ggplot2::element_line(linewidth = 1),
tern.axis.arrow.text = ggplot2::element_text(size = 6),
axis.text = ggplot2::element_text(size = 4),
axis.title = ggplot2::element_text(size = 6),
legend.text = ggplot2::element_text(size = 6),
legend.title = ggplot2::element_text(size = 8),
strip.text = ggplot2::element_text(size = 8),
plot.margin = unit(c(top,1,1,1)*0.2, "cm"),
strip.background = ggplot2::element_rect(color = NA),
) +
ggplot2::theme(tern.axis.arrow.text = element_text(size = 7),legend.position = "bottom", legend.margin = margin(r = 30), legend.box="vertical") +
ggplot2::scale_size_continuous(range = c(1,10),limits = c(0, r2max), breaks = seq(0.0001, r2max, length.out=5))+
ggplot2::guides(size = ggplot2::guide_legend(title = expression(R^2), order = 1, nrow = 1, label.position = "bottom"),color = ggplot2::guide_colorbar(title = env_deviance_name, title.position = "top", order = 2, barheight = 0.5, barwidth = 8))+
geom_point(alpha = 0.8, position= position_jitter_tern(x=0.02, y=0.02, z=0.02)) +
scale_color_viridis_c(na.value = "#666666",begin = 0,end = 1)+
ggtern::theme_nomask()
source("code/function_forest.R")
x=internals[[2]]
forest.figure(x)
library(grid)
library(forestploter)
data=x
data$p.env =data$env/(data$env+data$spa+data$codist)
data$p.spatial =data$spa/(data$env+data$spa+data$codist)
data$p.codist=data$codist/(data$env+data$spa+data$codist)
data.fish=data %>% filter(group=="Fish")
data.Amphibians=data %>% filter(group=="Amphibians")
data.Perching=data %>% filter(group=="Perching birds")
data.Domestic=data %>% filter(group=="Domestic animals")
data.Mammal=data %>% filter(group=="Mammal")
data.Waterfowl=data %>% filter(group=="Waterfowl")
data.list=list(data.fish,data.Amphibians,data.Perching,data.Domestic,data.Mammal,data.Waterfowl)
rm(data.fish,data.Amphibians,data.Perching,data.Domestic,data.Waterfowl)
data2 = data.frame(group = c('Fish','  environmental','  biotic_covariance','  spatial','Amphibians','  environmental','  biotic_covariance','  spatial','Perching birds','  environmental','  biotic_covariance','  spatial','Domestic animals','  environmental','  biotic_covariance','  spatial','Mammal','  environmental','  biotic_covariance','  spatial','Waterfowl','  environmental','  biotic_covariance','  spatial'), ciu=NA,cid=NA,cim=NA)
env=which(data2$group=="  environmental")
spatial=which(data2$group=="  spatial")
codist=which(data2$group=="  biotic_covariance")
for (i in 1:6) {
data2$cid[env[i]]=data.list[[i]]$p.env[which.min(data.list[[i]]$p.env)]
data2$ciu[env[i]]=data.list[[i]]$p.env[which.max(data.list[[i]]$p.env)]
data2$cim[env[i]]=median(data.list[[i]]$p.env,na.rm =TRUE)
data2$cid[spatial[i]]=data.list[[i]]$p.spatial[which.min(data.list[[i]]$p.spatial)]
data2$ciu[spatial[i]]=data.list[[i]]$p.spatial[which.max(data.list[[i]]$p.spatial)]
data2$cim[spatial[i]]=median(data.list[[i]]$p.spatial,na.rm =TRUE)
data2$cid[codist[i]]=data.list[[i]]$p.codist[which.min(data.list[[i]]$p.codist)]
data2$ciu[codist[i]]=data.list[[i]]$p.codist[which.max(data.list[[i]]$p.codist)]
data2$cim[codist[i]]=median(data.list[[i]]$p.codist,na.rm =TRUE)
}
View(data2)
View(data)
source("code/function_forest.R")
x=internals[[2]]
forest.figure(x)
data=x
data$p.env =data$env/(data$env+data$spa+data$codist)
data$p.spatial =data$spa/(data$env+data$spa+data$codist)
data$p.codist=data$codist/(data$env+data$spa+data$codist)
if (waterspecies=="FALSE"){
data.fish=data %>% filter(group=="Fish")
data.Amphibians=data %>% filter(group=="Amphibians")
data.Perching=data %>% filter(group=="Perching birds")
data.Domestic=data %>% filter(group=="Domestic animals")
data.Mammal=data %>% filter(group=="Mammal")
data.Waterfowl=data %>% filter(group=="Waterfowl")
data.list=list(data.fish,data.Amphibians,data.Perching,data.Domestic,data.Mammal,data.Waterfowl)
rm(data.fish,data.Amphibians,data.Perching,data.Domestic,data.Waterfowl)
data2 = data.frame(group = c('Fish','  environmental','  biotic_covariance','  spatial','Amphibians','  environmental','  biotic_covariance','  spatial','Perching birds','  environmental','  biotic_covariance','  spatial','Domestic animals','  environmental','  biotic_covariance','  spatial','Mammal','  environmental','  biotic_covariance','  spatial','Waterfowl','  environmental','  biotic_covariance','  spatial'), ciu=NA,cid=NA,cim=NA)
}
if (waterspecies=="TRUE"){
data.fish=data %>% filter(group=="Fish")
data.Amphibians=data %>% filter(group=="Amphibians")
data.list=list(data.fish,data.Amphibians)
rm(data.fish,data.Amphibians)
data2 = data.frame(group = c('Fish','  environmental','  biotic_covariance','  spatial','Amphibians','  environmental','  biotic_covariance','  spatial'), ciu=NA,cid=NA,cim=NA)
}
View(data2)
env=which(data2$group=="  environmental")
spatial=which(data2$group=="  spatial")
codist=which(data2$group=="  biotic_covariance")
for (i in 1:6) {
data2$cid[env[i]]=data.list[[i]]$p.env[which.min(data.list[[i]]$p.env)]
data2$ciu[env[i]]=data.list[[i]]$p.env[which.max(data.list[[i]]$p.env)]
data2$cim[env[i]]=median(data.list[[i]]$p.env,na.rm =TRUE)
data2$cid[spatial[i]]=data.list[[i]]$p.spatial[which.min(data.list[[i]]$p.spatial)]
data2$ciu[spatial[i]]=data.list[[i]]$p.spatial[which.max(data.list[[i]]$p.spatial)]
data2$cim[spatial[i]]=median(data.list[[i]]$p.spatial,na.rm =TRUE)
data2$cid[codist[i]]=data.list[[i]]$p.codist[which.min(data.list[[i]]$p.codist)]
data2$ciu[codist[i]]=data.list[[i]]$p.codist[which.max(data.list[[i]]$p.codist)]
data2$cim[codist[i]]=median(data.list[[i]]$p.codist,na.rm =TRUE)
}
if (waterspecies=="TRUE"){
for (i in 1:2) {
data2$cid[env[i]]=data.list[[i]]$p.env[which.min(data.list[[i]]$p.env)]
data2$ciu[env[i]]=data.list[[i]]$p.env[which.max(data.list[[i]]$p.env)]
data2$cim[env[i]]=median(data.list[[i]]$p.env,na.rm =TRUE)
data2$cid[spatial[i]]=data.list[[i]]$p.spatial[which.min(data.list[[i]]$p.spatial)]
data2$ciu[spatial[i]]=data.list[[i]]$p.spatial[which.max(data.list[[i]]$p.spatial)]
data2$cim[spatial[i]]=median(data.list[[i]]$p.spatial,na.rm =TRUE)
data2$cid[codist[i]]=data.list[[i]]$p.codist[which.min(data.list[[i]]$p.codist)]
data2$ciu[codist[i]]=data.list[[i]]$p.codist[which.max(data.list[[i]]$p.codist)]
data2$cim[codist[i]]=median(data.list[[i]]$p.codist,na.rm =TRUE)
}
data2$`median (95% CI)` <- ifelse(is.na(data2$ciu), "",
sprintf("%.2f (%.2f to %.2f)",
data2$cim, data2$cid, data2$ciu))
data2$` ` <- paste(rep(" ", 6), collapse = " ")
data2$`  ` <- paste(rep(" ", 6), collapse = " ")
}
View(data2)
p=forest(data2[,c(1,6,5)],
est = data2$cim,
lower = data2$cid,
upper = data2$ciu,
ci_column = 2,
vert_line = c(0.05)
)
g <- edit_plot(p,
row = c(1,5,9,13,17,21),
gp = gpar(fontface = "bold"))
g <- edit_plot(g,
row = env,
col = 2,
which = "ci",
gp = gpar(col = "#1FA187FF"))
g <- edit_plot(g, row = env, col=1,gp = gpar(col = "#1FA187FF"))
g <- edit_plot(g,
row = codist,
col = 2,
which = "ci",
gp = gpar(col = "#D55E00"))
g <- edit_plot(g, row = codist, col=1,gp = gpar(col = "#D55E00"))
plot(g)
source("code/function_forest.R")
x=internals[[2]]
forest.figure(x)
x=indexdata.OTU2$richness
p.env=indexdata.OTU2$env/(indexdata.OTU2$env+indexdata.OTU2$spa+indexdata.OTU2$codist)
#test the other env_deviance
Zscores = data.frame(ENV)
Zscores$HSI2_Pond_area[which(Zscores$HSI2_Pond_area == "0")]=0.01
Zscores$HSI2_Pond_area=log(Zscores$HSI2_Pond_area)
hist(data.frame(Zscores)$HSI2_Pond_area)
Zscores=as.matrix(Zscores)
Zscores=scales::rescale(Zscores, to=c(0,1))
env_deviance=rowSums(Zscores)
plot(env_deviance~p.env)
summary(lm(env_deviance~p.env))
ggplot(internals[[2]],aes(x =env, y =r2, color = group))+geom_point()+
geom_text(label=internals[[2]]$name,alpha = 0.6,size=3,col="black",check_overlap = F,show.legend = F,fontface = "italic")#codist,spa
species="15sp"
load(paste0("data/model/model_",species,"_AUC_20230602.rdata"))
spnames <- colnames(Y)
auc<-data.frame(sp=spnames,e.aucs=train_aucs,p.aucs=test_aucs,animalclass=TP$group,incidence=TP$inc.FHT,r2=internals[[2]]$r2)
auc$sp2=substring(auc$sp,6,)
auc <- auc %>%
mutate(
animalclass = case_when(
animalclass == "Amphibian" ~ "Amphibians",
animalclass == "Other Bird" ~ "Perching birds",
animalclass == "Domestic" ~ "Domestic animals",
TRUE ~ as.character(animalclass)
)
)
aucplot<- function(auc) {
mean.e=signif(mean(auc$e.aucs,na.rm = T),digits = 2)
mean.p=signif(mean(auc$p.aucs,na.rm = T),digits = 2)
auc$animalclass=factor(auc$animalclass,levels = c('Amphibians','Fish','Perching birds','Domestic animals','Mammal','Waterfowl'))
ggplot(auc,aes(x =e.aucs, y =p.aucs, color = animalclass,size=incidence))+
scale_color_viridis(option="D",discrete=TRUE)+
xlab("explanatory AUC") +
ylab("predictive AUC") +
xlim(0.3,1)+
ylim(0.3,1)+
geom_abline(intercept=0,slope=1,colour='azure3')+
geom_hline(yintercept = 0.5,colour='azure3')+
geom_vline(xintercept = 0.5,colour='azure3')+
geom_hline(yintercept = mean.p,linetype="dotted")+
geom_vline(xintercept = mean.e,linetype="dotted")+
labs(fill="Species Index") +
geom_point()+
ggplot2::annotate("text",label=mean.p, x = 0.45, y = mean.p+0.02)+
ggplot2::annotate("text",label=mean.e, x = 0.75, y = 0.4)+
theme(panel.grid.major =element_blank(), panel.grid.minor  =element_blank(),panel.background = element_blank(),panel.border = element_blank())+
geom_text_repel(aes(e.aucs, p.aucs, label=sp2,size=20,fontface = "italic"),show.legend = FALSE)
}
aucplot(auc)
species="15sp"
load(paste0("data/model/model_",species,"_AUC_20230602.rdata"))
spnames <- colnames(Y)
auc<-data.frame(sp=spnames,e.aucs=train_aucs,p.aucs=test_aucs,animalclass=TP$group,incidence=TP$inc.FHT,r2=internals[[2]]$r2)
auc$sp2=substring(auc$sp,6,)
auc <- auc %>%
mutate(
animalclass = case_when(
animalclass == "Amphibian" ~ "Amphibians",
animalclass == "Other Bird" ~ "Perching birds",
animalclass == "Domestic" ~ "Domestic animals",
TRUE ~ as.character(animalclass)
)
)
aucplot<- function(auc) {
mean.e=signif(mean(auc$e.aucs,na.rm = T),digits = 2)
mean.p=signif(mean(auc$p.aucs,na.rm = T),digits = 2)
auc$animalclass=factor(auc$animalclass,levels = c('Amphibians','Fish','Perching birds','Domestic animals','Mammal','Waterfowl'))
ggplot(auc,aes(x =e.aucs, y =p.aucs, color = animalclass,size=incidence))+
scale_color_viridis(option="D",discrete=TRUE)+
xlab("explanatory AUC") +
ylab("predictive AUC") +
xlim(0.3,1)+
ylim(0.3,1)+
geom_abline(intercept=0,slope=1,colour='azure3')+
geom_hline(yintercept = 0.5,colour='azure3')+
geom_vline(xintercept = 0.5,colour='azure3')+
geom_hline(yintercept = mean.p,linetype="dotted")+
geom_vline(xintercept = mean.e,linetype="dotted")+
labs(fill="Species Index") +
geom_point()+
ggplot2::annotate("text",label=mean.p, x = 0.45, y = mean.p+0.02)+
ggplot2::annotate("text",label=mean.e, x = 0.75, y = 0.4)+
theme(panel.grid.major =element_blank(), panel.grid.minor  =element_blank(),panel.background = element_blank(),panel.border = element_blank())+
geom_text_repel(aes(e.aucs, p.aucs, label=sp2,size=20,fontface = "italic"),show.legend = FALSE)
}
aucplot(auc)
rm(list=ls())
species="15sp"
load(paste0("data/model/model_",species,"_AUC_20230602.rdata"))
dataselect="FHT"
load("data/FHT_320_48sp_20211110.rdata")
TP<-read.csv("data/FHT_DATA_TP_20211110.csv", header=TRUE)
waterspecies="TRUE"  #TRUE or FALSE
species="15sp"
if (waterspecies=="TRUE"){
Y=data.frame(Y)%>%select(starts_with("Acti") | starts_with("Amph"))
which(colSums(Y) == 0)
ex=which(rowSums(Y) == 0)
Y=Y[-ex,]
env=env[-ex,]
sp=sp[-ex,]
TP<-TP %>% filter(group =="Fish" | group =="Amphibian" )
}
Y=as.matrix(Y)
which(colSums(Y) == 0)
which(colSums(Y) < 5)
TP$species_name=substring(TP$species,6,)
spnames <- colnames(Y)
auc15<-data.frame(sp=spnames,e.aucs15=train_aucs,p.aucs15=test_aucs,animalclass=TP$group,incidence=TP$inc.FHT)
auc15$sp2=substring(auc15$sp,6,)
auc15 <- auc15 %>%
mutate(
animalclass = case_when(
animalclass == "Amphibian" ~ "Amphibians",
animalclass == "Other Bird" ~ "Perching birds",
animalclass == "Domestic" ~ "Domestic animals",
TRUE ~ as.character(animalclass)
)
)
species="48sp"
load(paste0("data/model/model_",species,"_AUC_20230602.rdata"))
dataselect="FHT"
load("data/FHT_320_48sp_20211110.rdata")
TP<-read.csv("data/FHT_DATA_TP_20211110.csv", header=TRUE)
waterspecies="FALSE"  #TRUE or FALSE
species="48sp"
Y=as.matrix(Y)
which(colSums(Y) == 0)
which(colSums(Y) < 5)
TP$species_name=substring(TP$species,6,)
spnames <- colnames(Y)
auc48<-data.frame(sp=spnames,e.aucs48=train_aucs,p.aucs48=test_aucs,animalclass=TP$group,incidence=TP$inc.FHT)
auc48$sp2=substring(auc48$sp,6,)
auc48 <- auc48 %>%
mutate(
animalclass = case_when(
animalclass == "Amphibian" ~ "Amphibians",
animalclass == "Other Bird" ~ "Perching birds",
animalclass == "Domestic" ~ "Domestic animals",
TRUE ~ as.character(animalclass)
)
)
auc=left_join(auc15,auc48,by="sp")
auc$animalclass=factor(auc$animalclass,levels = c('Amphibians','Fish','Perching birds','Domestic animals','Mammal','Waterfowl'))
rm(list=ls())
species="15sp"
load(paste0("data/model/model_",species,"_AUC_20230602.rdata"))
dataselect="FHT"
load("data/FHT_320_48sp_20211110.rdata")
TP<-read.csv("data/FHT_DATA_TP_20211110.csv", header=TRUE)
waterspecies="TRUE"  #TRUE or FALSE
species="15sp"
if (waterspecies=="TRUE"){
Y=data.frame(Y)%>%select(starts_with("Acti") | starts_with("Amph"))
which(colSums(Y) == 0)
ex=which(rowSums(Y) == 0)
Y=Y[-ex,]
env=env[-ex,]
sp=sp[-ex,]
TP<-TP %>% filter(group =="Fish" | group =="Amphibian" )
}
Y=as.matrix(Y)
which(colSums(Y) == 0)
which(colSums(Y) < 5)
TP$species_name=substring(TP$species,6,)
spnames <- colnames(Y)
auc15<-data.frame(sp=spnames,e.aucs15=train_aucs,p.aucs15=test_aucs,animalclass=TP$group,incidence=TP$inc.FHT)
auc15$sp2=substring(auc15$sp,6,)
auc15 <- auc15 %>%
mutate(
animalclass = case_when(
animalclass == "Amphibian" ~ "Amphibians",
animalclass == "Other Bird" ~ "Perching birds",
animalclass == "Domestic" ~ "Domestic animals",
TRUE ~ as.character(animalclass)
)
)
species="48sp"
load(paste0("data/model/model_",species,"_AUC_20230602.rdata"))
dataselect="FHT"
load("data/FHT_320_48sp_20211110.rdata")
TP<-read.csv("data/FHT_DATA_TP_20211110.csv", header=TRUE)
waterspecies="FALSE"  #TRUE or FALSE
species="48sp"
Y=as.matrix(Y)
which(colSums(Y) == 0)
which(colSums(Y) < 5)
TP$species_name=substring(TP$species,6,)
spnames <- colnames(Y)
auc48<-data.frame(sp=spnames,e.aucs48=train_aucs,p.aucs48=test_aucs,animalclass=TP$group,incidence=TP$inc.FHT)
auc48$sp2=substring(auc48$sp,6,)
auc48 <- auc48 %>%
mutate(
animalclass = case_when(
animalclass == "Amphibian" ~ "Amphibians",
animalclass == "Other Bird" ~ "Perching birds",
animalclass == "Domestic" ~ "Domestic animals",
TRUE ~ as.character(animalclass)
)
)
auc=left_join(auc15,auc48,by="sp")
auc$animalclass=factor(auc$animalclass,levels = c('Amphibians','Fish','Perching birds','Domestic animals','Mammal','Waterfowl'))
