-
Notifications
You must be signed in to change notification settings - Fork 0
/
TPcorrection.r
66 lines (54 loc) · 2.88 KB
/
TPcorrection.r
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
library(MASS)
library(ggplot2)
library(RColorBrewer)
setwd("E:/Richard/Global_Hypsometry/GIS/sinusoidal/T_P_correction")
# NORMAL CORRECTION ###########################################
Cvars <- c("P",'T','tetrapods')
Cdata <- read.table("BID_for_lm.txt", header = FALSE, sep = ",", col.names = Cvars)
# Tetrapods
ClinMod <- lm(Cdata$tetrapods ~ Cdata$T + Cdata$P, data= Cdata)
Clm_results <- summary(ClinMod)
# Empirical expected tetrapod richness at every T/P pair
CTETemp <- Clm_results$coefficients[[1,1]] + Cdata$T * Clm_results$coefficients[[2,1]] + Cdata$P * Clm_results$coefficients[[3,1]]
# correct tetrapods
CTETcorr <- Cdata$tetrapods - CTETemp
# write(CTETcorr, file = "BIDcorr_for_matlab.txt", sep = ",")
# write.table(CTETcorr, "BIDcorr_for_matlab.txt", append = FALSE, sep = ",", dec = ".",
# row.names = F, col.names = F)
# THRESHOLD CORRECTION ########################################
###############################################################
Cvars <- c("P",'T','tetrapods')
Cdata <- read.table("MAM_for_lm.txt", header = FALSE, sep = ",", col.names = Cvars)
Cdata$P.thres = Cdata$P
Cdata$P.thres[Cdata$P > 3500] = 3500
# Tetrapods
ClinModT <- lm(Cdata$tetrapods ~ Cdata$T + Cdata$P.thres, data= Cdata)
Clm_resultsT <- summary(ClinModT)
# Empirical expected tetrapod richness at every T/P pair
CTETempT <- Clm_resultsT$coefficients[[1,1]] + Cdata$T * Clm_resultsT$coefficients[[2,1]] + Cdata$P.thres * Clm_resultsT$coefficients[[3,1]]
# correct tetrapods
CTETcorrT <- Cdata$tetrapods - CTETempT
# write(CTETcorrT, file = "BIDcorrThreshold_for_matlab.txt", sep = ",")
# write.table(CTETcorrT, "MAMcorrThreshold_for_matlab.txt", append = FALSE, sep = ",", dec = ".",
# row.names = F, col.names = F)
# PLOT BI-RELATIONSSHIPS ######################################
# T- BID ######################################################
TF <-ggplot(data = Cdata, aes(x=Cdata$T, y=Cdata$tetrapods)) +
stat_density_2d(aes(fill = ..level..), geom = "polygon", colour="white")
geom_bin2d(bins = 100)+geom_contour()
scale_fill_distiller(palette= "Spectral", direction=1) +
scale_x_continuous(expand = c(0, 0)) +
scale_y_continuous(expand = c(0, 0))
theme_bw()
TF + labs(title = "T-correlation",x = "T",y= "Tetrapod richness")
# P- BID ######################################################
P <-ggplot(data = Cdata, aes(x=Cdata$P, y=Cdata$tetrapods)) +
geom_density_2d()
# stat_density_2d(aes(fill = ..level..), geom = "polygon", colour="white")
# geom_bin2d(bins = 200) + geom_contour(data = Cdata)+
# stat_density_2d(aes(fill = ..density..), geom = "raster", contour = T) +
# scale_fill_distiller(palette= "Spectral", direction=1) +
# scale_x_continuous(expand = c(0, 0)) +
# scale_y_continuous(expand = c(0, 0)) +
theme_bw()
P + labs(title = "P-correlation",x = "P",y= "Tetrapod richness")