diff --git a/code/Trait-Environment-Relationships1.R b/code/Trait-Environment-Relationships1.R
deleted file mode 100644
index 7e20b4e55785dc078385f0ba845ecd267ec5439a..0000000000000000000000000000000000000000
--- a/code/Trait-Environment-Relationships1.R
+++ /dev/null
@@ -1,1859 +0,0 @@
-getwd()
-###### Traits #####
-load("data/TRY.all.mean.sd.2.Rdata")
-str(TRY.all.mean.sd.2)
-# 40791 obs. of 38 variables
-# tests
-
-plot(SLA.mean~LDMC.mean, data=TRY.all.mean.sd.2)
-#summary(lm(SLA.mean~LDMC.mean, data=TRY.all.mean.sd.2))
-#abline(lm(SLA.mean~LDMC.mean, data=TRY.all.mean.sd.2),col="blue", lwd=2)
-#sqrt(0.3415)
-
-###### Backbone #####
-load("data/backbone.v.2.splot.try3.Rdata")
-# alternatively:
-#sPlot3 <- read.csv(paste("data/", "backbone.v.2.splot.try3.csv",
-# sep = ""), sep=",", quote = "\"'")
-str(backbone.splot.try3)
-# data.frame': 122901 obs. of 33 variables:
-
-###### Harmonized plot data ######
-DT <- read.csv("data/plots_without_unified_cover_scales_20160120a.csv")
-library(data.table)
-DT <- data.table(DT)
-colnames(DT)
-colnames(DT)[9] <- "Cover"
-setkey(DT,PlotObservationID)
-
-###### Header data #######
-splot.header <- read.csv("data/sPlot_14_04_2015_coord.csv")
-str(splot.header)
-'data.frame: 1117940 obs. of 3 variables:
-$ PlotObservationID: int 1 2 3 4 5 6 7 8 9 10 ...
-$ Longitude : num -154 -154 -154 -154 -154 ...
-$ Latitude : num 62.4 62.4 62.4 62.4 62.4 ...'
-
-##### Bioclim data ####
-splot.bioclim <- read.csv(paste("data/","sPlot.csv", sep = ""), sep = ";", na.strings=c("","NA"))
-str(splot.bioclim)
-#data.frame': 1117382 obs. of 52 variables:
-
-####################
-##### Matching #####
-index6 <- match(TRY.all.mean.sd.2$StandSpeciesName,backbone.splot.try3$names.sPlot.TRY)
-any(is.na(index6)) #T
-# not all species in the trait list are also in the backbone
-TRY.all.mean.sd.2$StandSpeciesName[is.na(index6)]
-# 3959 species in TRY are not in the backbone
-
-### data cleaning ###
-length(TRY.all.mean.sd.2$StandSpeciesName[is.na(TRY.all.mean.sd.2$StandSpeciesName)])
-# 1 record is NA
-TRY.all.mean.sd.2[which(is.na(TRY.all.mean.sd.2$StandSpeciesName)),]
-# record no. 7947 has StandSpeciesName = NA
-TRY.all.mean.sd.2 <- TRY.all.mean.sd.2[-which(is.na(TRY.all.mean.sd.2$StandSpeciesName)),]
-str(TRY.all.mean.sd.2) #40790 obs. of 38 variables
-# 1 record less
-'data.frame: 40790 obs. of 38 variables:
-$ StandSpeciesName : Factor w/ 40790 levels "Abarema adenophora",..: 965 3066 3128 3445 4239 7449 7467 7829 8665 12419 ...
-$ n : int 23 32 15 8 21 32 20 31 7 8 ...
-$ SLA.mean : num 1.42 1.83 1.92 2.7 2.31 ...
-$ PlantHeight.mean : num 1.282 2.92 0.585 0.488 1.1 ...
-$ SeedMass.mean : num 0.0656 2.7757 2.6427 -2.6794 -2.4111 ...
-$ LDMC.mean : num -0.765 -1.367 -1.009 -1.228 -0.963 ...
-$ StemDens.mean : num -0.298 -0.524 -0.322 -0.69 -0.592 ...
-$ LeafArea.mean : num 3.39 8.36 5.26 3.28 5.16 ...
-$ LeafN.mean : num 2.33 2.5 2.21 2.95 3.17 ...
-$ LeafP.mean : num 0.118 0.074 -0.127 0.792 0.421 ...
-$ LeafNperArea.mean : num 0.692 0.485 0.381 0.452 0.709 ...
-$ Leaffreshmass.mean : num -2.325 0.112 -1.916 -3.777 -2.714 ...
-$ LeafNPratio.mean : num 2.23 2.32 2.24 2.14 2.81 ...
-$ LeafC.perdrymass.mean : num 6.14 6.25 6.31 6.1 6.12 ...
-$ Leaf.delta.15N.mean : num 2.68 2.44 2.32 2.78 2.72 ...
-$ Stem.cond.dens.mean : num 5.44 4.82 5.44 6.41 6.69 ...
-$ Seed.num.rep.unit.mean : num 4.57 4.34 1.48 9.87 8.48 ...
-$ Wood.vessel.length.mean: num 6.32 6.17 6.18 6.13 6.16 ...
-$ Seed.length.mean : num 1.18 1.45 1.78 0.16 1.21 ...
-$ Disp.unit.leng.mean : num 1.173 2.119 2.064 0.039 1.287 ...
-$ SLA.sd : num 0.1105 0.1754 0.3161 0.0694 0.0596 ...
-$ PlantHeight.sd : num 0.185 0.347 1.296 0.274 0.137 ...
-$ SeedMass.sd : num 0.288 1.016 0.234 0.201 0.317 ...
-$ LDMC.sd : num 0.038 0.0469 0.0449 0.0433 0.0321 ...
-$ StemDens.sd : num 0.0629 0.0582 0.0511 0.0787 0.0635 ...
-$ LeafArea.sd : num 0.876 0.235 1.89 0.135 0.196 ...
-$ LeafN.sd : num 0.0857 0.0902 0.1229 0.0493 0.0443 ...
-$ LeafP.sd : num 0.2173 0.0974 0.142 0.0356 0.0611 ...
-$ LeafNperArea.sd : num 0.0983 0.0916 0.065 0.0425 0.053 ...
-$ Leaffreshmass.sd : num 0.454 0.189 0.776 0.116 0.144 ...
-$ LeafNPratio.sd : num 0.1921 0.0522 0.1184 0.032 0.0477 ...
-$ LeafC.perdrymass.sd : num 0.00741 0.00767 0.00706 0.00634 0.00469 ...
-$ Leaf.delta.15N.sd : num 0.0363 0.0376 0.0382 0.0296 0.0198 ...
-$ Stem.cond.dens.sd : num 0.0996 0.3641 0.2232 0.1064 0.1064 ...
-$ Seed.num.rep.unit.sd : num 0.319 0.51 0.684 0.295 0.408 ...
-$ Wood.vessel.length.sd : num 0.0435 0.0518 0.0673 0.0392 0.0291 ...
-$ Seed.length.sd : num 0.0508 0.0744 0.0992 0.0817 0.0651 ...
-$ Disp.unit.leng.sd : num 0.0721 0.1367 0.1409 0.0423 0.0754 ...
-'
-
-index6 <- match(TRY.all.mean.sd.2$StandSpeciesName,backbone.splot.try3$names.sPlot.TRY)
-any(is.na(index6)) #T
-TRY.all.mean.sd.2$StandSpeciesName[is.na(index6)]
-# 3958 species
-TRY.all.mean.sd.2$species <- backbone.splot.try3$name.short.correct[index6]
-dim(TRY.all.mean.sd.2) # 40790 39
-
-
-index1 <- match(DT$PlotObservationID, splot.header$PlotObservationID)
-any(is.na(index1)) #F
-length(index1) # 24241941
-length(unique(index1[!is.na(index1)])) #1117940
-
-index2 <- match(splot.header$PlotObservationID,DT$PlotObservationID)
-any(is.na(index2)) #F
-length(index2) # 1117940
-length(unique(index2[!is.na(index2)])) # 1117940
-dim(splot.header)
-# 1117940 3
-identical(as.character(splot.header$PlotObservationID),as.character(DT$PlotObservationID[index2]))
-# T
-
-any(is.na(DT$species)) #T
-# that should not be the case!
-length(DT$species[!is.na(DT$species)]) #24221565
-length(DT$species) #24241941
-24241941-24221565 # 20376 NA names!!!
-# it gives nonsense to match them with traits
-index7 <- match(DT$species,TRY.all.mean.sd.2$species)
-length(index7) #24241941
-length(index7[!is.na(index7)]) #21040927, with previous TRY version that was: 19841429
-24241941 - 21040927 # 3201014 entries have no traits
-(24241941 - 3201014)/24241941*100
-# which are 86.79555% of plots that have traits,
-# previously: 18.15247 %
-
-# reduce DT, and splot.header
-#splot.header <- splot.header[!is.na(DT$species[index2]),]
-#dim(splot.header) #1117476 3
-DT <- DT[!is.na(DT$species),]
-index7 <- match(DT$species,TRY.all.mean.sd.2$species)
-length(index7) #24221565
-
-
-### CWM ###
-mean(TRY.all.mean.sd.2$SLA.mean) # 2.611657
-min(TRY.all.mean.sd.2$SLA.mean) # -1.297031
-max(TRY.all.mean.sd.2$SLA.mean) # 5.893029
-# example
-DT$trait <- NA
-DT$trait <- TRY.all.mean.sd.2$SLA.mean[index7]
-length(DT$trait[!is.na(DT$trait)]) # 21020551
-length(DT$trait[is.na(DT$trait)]) # 3201014
-str(DT)
-
-colnames(TRY.all.mean.sd.2)
-CWM2 <- DT[,list(CWM.SLA = weighted.mean(trait,Cover,na.rm = T)),by=PlotObservationID]
-# I have checked that we do not need to prefilter only those entries that have a trait value
-# this works fine
-str(CWM2)
-dim(CWM2) #1117898
-which(colnames(TRY.all.mean.sd.2)=="SLA.mean") # 3
-which(colnames(TRY.all.mean.sd.2)=="Disp.unit.leng.mean") # 20
-CWM <- array(NA,c(dim(CWM2)[1],18),dimnames=list(CWM2$PlotObservationID,
- colnames(TRY.all.mean.sd.2)[3:20]))
-rm(CWM2)
-for (i in 1:18){
- DT$trait <- NA
- DT$trait <- TRY.all.mean.sd.2[index7,i+2]
- CWM[,i] <- DT[,list(CWM.trait= weighted.mean(trait,Relative.cover,na.rm = T)),by=PlotObservationID]$CWM.trait
-
-}
-str(CWM) #num [1:1117898, 1:18]
-CWM[1:20,]
-tail(CWM)
-write.csv(CWM,file = "CWM.csv", row.names = T)
-
-CWM <- read.csv("C://Daten//iDiv2//splot2//CWM.csv")
-str(CWM) #1117898 obs. of 19 variables
-dimnames(CWM)[[1]] <- CWM$X
-CWM <- CWM[,-1]
-tail(CWM)
-
-##### MAP ####
-index3 <- match(dimnames(CWM)[[1]],splot.header$PlotObservationID)
-length(index3) #1117898
-any(is.na(index3)) # F
-#length(index3[is.na(index3)]) #422
-splot.header2 <- splot.header[index3,]
-'CWM1$Longitude <- splot.header$Longitude[index10]
-CWM1$Latitude <- splot.header$Latitude[index10]
-'
-str(CWM) # num [1:1117898, 1:18]
-CWM <- CWM[!(is.na(splot.header2$Longitude)|is.na(splot.header2$Latitude)|is.na(CWM[,"SLA.mean"])),]
-str(CWM) # num [1:1111307, 1:18]
-1111307/1117898*100 # 99.41041% of plots have traits, Long and Lat
-any(is.na(CWM)) #F
-# match header again
-index4 <- match(dimnames(CWM)[[1]],splot.header$PlotObservationID)
-length(index4) #1111307
-any(is.na(index4)) # F
-splot.header2 <- splot.header[index4,]
-# now in the same sequence as CWM
-
-library(raster)
-CRSlonlat <- CRS("+proj=longlat +datum=WGS84")
-coords <- cbind(splot.header2$Longitude, splot.header2$Latitude)
-coords <- SpatialPoints(coords, proj4string=CRSlonlat)
-map.CWM <- SpatialPointsDataFrame(coords, as.data.frame(CWM), proj4string=CRSlonlat)
-str(map.CWM)
-
-library(maptools)
-data(wrld_simpl)
-plot(wrld_simpl)
-points(coordinates(map.CWM), cex=0.1, col="red", pch=16)
-# sPlot2withSLA.png
-max(coordinates(map.CWM)[,1]) # 179.5901
-max(coordinates(map.CWM)[,2]) # 80.14912
-min(coordinates(map.CWM)[,1]) # -162.7414
-min(coordinates(map.CWM)[,2]) # -64.78
-
-
-r <- raster(nrows=180, ncols=360, xmn=-180, xmx=180, ymn=-90, ymx=90)
-str(r)
-raster.map.SLA <- rasterize(map.CWM,r,map.CWM$SLA.mean, run="mean")
-str(raster.map.SLA)
-plot(wrld_simpl)
-plot(raster.map.SLA, add=T)
-# SLAgridded1.png
-
-library(rworldmap)
-raster.map.SLA.grid <- as(raster.map.SLA, 'SpatialGridDataFrame')
-library(RColorBrewer)
-# Set a color ramp
-rf <- brewer.pal(8, 'Spectral')
-mapParams <- mapGriddedData(raster.map.SLA.grid,addLegend=F,
- numCats = 8, colourPalette=rf)
-str(mapParams)
-mapParams$cutVector
-# addMapLegend(plottedData=raster.map.SLA.grid,legendShrink=0.5)
-#do.call(addMapLegend, c( mapParams, legendLabels="all", legendWidth=0.5,
-# legendIntervals="data", legendMar = 2 ) )
-do.call(addMapLegend, c( mapParams, legendLabels="all", legendWidth=0.3,
- legendIntervals="page", legendMar = 4, legendShrink=1.7,horizontal=F))
-#round(as.numeric(exp(mapParams$cutVector)),1)
-mtext(round(as.numeric(exp(mapParams$cutVector)),1),side=4,line=-2,
- at=seq(-150,170, by=39), las=2)
-#mtext("ln(SLA)",side=1,line=0, at=230, las=1)
-#mtext("ln(SLA)",side=4,line=-2,at=-170, las=2)
-# SLAgridded2.png
-
-############################################################
-# Regression with Bioclim
-############################################################
-'BIO1 = Annual Mean Temperature
-BIO2 = Mean Diurnal Range (Mean of monthly (max temp - min temp))
-BIO3 = Isothermality (BIO2/BIO7) (* 100)
-BIO4 = Temperature Seasonality (standard deviation *100)
-BIO5 = Max Temperature of Warmest Month
-BIO6 = Min Temperature of Coldest Month
-BIO7 = Temperature Annual Range (BIO5-BIO6)
-BIO8 = Mean Temperature of Wettest Quarter
-BIO9 = Mean Temperature of Driest Quarter
-BIO10 = Mean Temperature of Warmest Quarter
-BIO11 = Mean Temperature of Coldest Quarter
-BIO12 = Annual Precipitation
-BIO13 = Precipitation of Wettest Month
-BIO14 = Precipitation of Driest Month
-BIO15 = Precipitation Seasonality (Coefficient of Variation)
-BIO16 = Precipitation of Wettest Quarter
-BIO17 = Precipitation of Driest Quarter
-BIO18 = Precipitation of Warmest Quarter
-BIO19 = Precipitation of Coldest Quarter'
-
-index11 <- match(splot.header2$PlotObservationID,splot.bioclim$PLOT_ID)
-length(index11) # 1111307
-length(index11[!is.na(index11)]) # 1111307
-
-#index12 <- match(CWM2$PlotObservationID, splot.header$PlotObservationID)
-#any(is.na(index12)) #F
-str(splot.bioclim$BIO_01)
-plot(CWM[,"SLA.mean"]~splot.bioclim$BIO_01[index11],
- cex=0.05, xlab="BIO_1", ylab="CWM ln(SLA)",cex.lab=1.5, cex.axis=1.5)
-model1 <- lm(CWM[,"SLA.mean"]~splot.bioclim$BIO_01[index11])
-model2 <- lm(CWM[,"SLA.mean"]~splot.bioclim$BIO_01[index11]+I(splot.bioclim$BIO_01[index11]^2))
-summary(model1)
-'Coefficients:
- Estimate Std. Error t value Pr(>|t|)
-(Intercept) 2.945e+00 9.167e-04 3212.9 <2e-16 ***
-splot.bioclim$BIO_01[index11] -1.204e-03 9.067e-06 -132.8 <2e-16 ***
-Multiple R-squared: 0.01562, Adjusted R-squared: 0.01562
-F-statistic: 1.764e+04 on 1 and 1111305 DF, p-value: < 2.2e-16
-'
-summary(model2)
-' Estimate Std. Error t value Pr(>|t|)
-(Intercept) 2.806e+00 1.113e-03 2520.3 <2e-16 ***
-splot.bioclim$BIO_01[index11] 2.188e-03 1.857e-05 117.8 <2e-16 ***
-I(splot.bioclim$BIO_01[index11]^2) -1.657e-05 7.985e-08 -207.6 <2e-16 ***
-Multiple R-squared: 0.05241
-'
-abline(model1, lwd=3, col="blue")
-model2$coefficients
-x.new <- seq(-100,300,1)
-y.new <- model2$coefficients[1]+model2$coefficients[2]*x.new+model2$coefficients[3]*x.new^2
-lines(x.new,y.new, lwd=3, col="blue")
-# SLA_BIO1.png
-
-# devide in quantiles
-quantile(splot.bioclim$BIO_01[index11],probs = seq(0, 1, 1/3))
-' 0% 33.33333% 66.66667% 100%
- -227 79 99 303 '
-# -> most points (i.e. one third) are in the temperate zone between 7.9 and 9.9°C
-x1 <- splot.bioclim$BIO_01[index11][splot.bioclim$BIO_01[index11]<=79]
-model1a <- lm(CWM[,"SLA.mean"][splot.bioclim$BIO_01[index11]<=79]~splot.bioclim$BIO_01[index11][splot.bioclim$BIO_01[index11]<=79])
-summary(model1a)
-'Coefficients:
- Estimate Std. Error t value Pr(>|t|)
-(Intercept) 2.697e+00 1.244e-03 2167.2 <2e-16 ***
-splot.bioclim$BIO_01[index11][splot.bioclim$BIO_01[index11] <= 79] 2.466e-03 2.117e-05 116.5 <2e-16 ***
----
-Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
-
-Residual standard error: 0.4027 on 371873 degrees of freedom
-Multiple R-squared: 0.03519, Adjusted R-squared: 0.03519
-F-statistic: 1.356e+04 on 1 and 371873 DF, p-value: < 2.2e-16'
-model1b <- lm(CWM[,"SLA.mean"][splot.bioclim$BIO_01[index11]<=99]~splot.bioclim$BIO_01[index11][splot.bioclim$BIO_01[index11]<=99])
-summary(model1b)
-' Estimate Std. Error t value Pr(>|t|)
-(Intercept) 2.692e+00 1.161e-03 2319.8 <2e-16 ***
-splot.bioclim$BIO_01[index11][splot.bioclim$BIO_01[index11] <= 99] 2.777e-03 1.527e-05 181.9 <2e-16 ***
----
-Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
-
-Residual standard error: 0.3954 on 750010 degrees of freedom
-Multiple R-squared: 0.04224, Adjusted R-squared: 0.04224
-F-statistic: 3.308e+04 on 1 and 750010 DF, p-value: < 2.2e-16'
-model1c <- lm(CWM[,"SLA.mean"][splot.bioclim$BIO_01[index11]>99]~splot.bioclim$BIO_01[index11][splot.bioclim$BIO_01[index11]>99])
-summary(model1c)
-' Estimate Std. Error t value Pr(>|t|)
-(Intercept) 3.236e+00 2.587e-03 1251.1 <2e-16 ***
-splot.bioclim$BIO_01[index11][splot.bioclim$BIO_01[index11] > 99] -3.774e-03 1.855e-05 -203.4 <2e-16 ***
----
-Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
-
-Residual standard error: 0.443 on 361293 degrees of freedom
-Multiple R-squared: 0.1028, Adjusted R-squared: 0.1028
-F-statistic: 4.138e+04 on 1 and 361293 DF, p-value: < 2.2e-16
-'
-
-plot(CWM[,"SLA.mean"]~splot.bioclim$BIO_01[index11],
- cex=0.05, xlab="BIO_1", ylab="CWM ln(SLA)",cex.lab=1.5, cex.axis=1.5)
-x.new <- seq(-100,79,1)
-y.new <- model1a$coefficients[1]+model1a$coefficients[2]*x.new
-lines(x.new,y.new, lwd=3, col="blue")
-x.new <- seq(79,99,1)
-y.new <- model1b$coefficients[1]+model1b$coefficients[2]*x.new
-lines(x.new,y.new, lwd=3, col="green")
-x.new <- seq(99,300,1)
-y.new <- model1c$coefficients[1]+model1c$coefficients[2]*x.new
-lines(x.new,y.new, lwd=3, col="red")
-
-### calculate a GAM
-library(gam)
-x <- splot.bioclim$BIO_01[index11]
-model3 <- gam(CWM[,"SLA.mean"]~lo(x, span=0.5))
-summary(model3)
-' Null Deviance: 207269.5 on 1111306 degrees of freedom
-Residual Deviance: 187452.7 on 1111302 degrees of freedom
-AIC: 1175898
-Anova for Parametric Effects
- Df Sum Sq Mean Sq F value Pr(>F)
-lo(x, span = 0.5) 1 3238 3238.0 19197 < 2.2e-16 ***
-Residuals 1111302 187453 0.2
-Anova for Nonparametric Effects
- Npar Df Npar F Pr(F)
-(Intercept)
-lo(x, span = 0.5) 2.7 37028 < 2.2e-16 ***
-'
-
-x.new <- seq(-100,300,1)
-y.new <- predict(model3, data.frame(x=x.new), se = F)
-# se=T does not work
-plot(CWM[,"SLA.mean"]~splot.bioclim$BIO_01[index11],
- cex=0.05, xlab="BIO_1", ylab="CWM ln(SLA)",cex.lab=1.5, cex.axis=1.5)
-lines(x.new,y.new, lwd=3, col="blue")
-
-### N to P ratio
-model4 <- lm(CWM[,"LeafNPratio.mean"]~splot.bioclim$BIO_01[index11])
-summary(model4)
-' Estimate Std. Error t value Pr(>|t|)
-(Intercept) 2.290e+00 4.334e-04 5284.0 <2e-16 ***
-splot.bioclim$BIO_01[index11] 1.531e-03 4.287e-06 357.2 <2e-16 ***
----
-Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
-
-Residual standard error: 0.2026 on 1111305 degrees of freedom
-Multiple R-squared: 0.103, Adjusted R-squared: 0.103
-'
-plot(CWM[,"LeafNPratio.mean"]~splot.bioclim$BIO_01[index11],
- cex=0.05, xlab="BIO_01", ylab="CWM ln(Leaf N to P ratio)",cex.lab=1.5, cex.axis=1.5)
-abline(model4, lwd=3, col="blue")
-
-### turn the regression around
-# does not make sense
-model5 <- lm(splot.bioclim$BIO_01[index11]~CWM[,"SLA.mean"])
-summary(model5)
-' Estimate Std. Error t value Pr(>|t|)
-(Intercept) 127.41676 0.28026 454.6 <2e-16 ***
-CWM[, "SLA.mean"] -12.97407 0.09769 -132.8 <2e-16 ***
-Multiple R-squared: 0.01562, Adjusted R-squared: 0.01562 '
-
-model6 <- lm(splot.bioclim$BIO_01[index11]~CWM[,"LeafNPratio.mean"])
-summary(model6)
-'Coefficients:
- Estimate Std. Error t value Pr(>|t|)
-(Intercept) -72.7479 0.4591 -158.5 <2e-16 ***
-CWM[, "LeafNPratio.mean"] 67.2535 0.1883 357.2 <2e-16 ***
----
-Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
-
-Residual standard error: 42.46 on 1111305 degrees of freedom
-Multiple R-squared: 0.103, Adjusted R-squared: 0.103 '
-
-model7 <- lm(splot.bioclim$BIO_01[index11]~CWM[,"SLA.mean"]+CWM[,"LeafNPratio.mean"])
-summary(model7)
-' Estimate Std. Error t value Pr(>|t|)
-(Intercept) -41.59378 0.55839 -74.49 <2e-16 ***
-CWM[, "SLA.mean"] -9.08978 0.09355 -97.17 <2e-16 ***
-CWM[, "LeafNPratio.mean"] 65.04103 0.18885 344.40 <2e-16 ***
----
-Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
-
-Residual standard error: 42.28 on 1111304 degrees of freedom
-Multiple R-squared: 0.1106, Adjusted R-squared: 0.1106 '
-x <- splot.bioclim$BIO_01[index11]
-model8 <- lm(x~CWM[,"LeafNPratio.mean"]*CWM[,"SLA.mean"])
-summary(model8)
-'Coefficients:
- Estimate Std. Error t value Pr(>|t|)
-(Intercept) -432.2500 2.2645 -190.9 <2e-16 ***
-CWM[, "LeafNPratio.mean"] 227.1681 0.9304 244.2 <2e-16 ***
-CWM[, "SLA.mean"] 135.8572 0.8202 165.6 <2e-16 ***
-CWM[, "LeafNPratio.mean"]:CWM[, "SLA.mean"] -60.2280 0.3386 -177.9 <2e-16 ***
----
-Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
-
-Residual standard error: 41.69 on 1111303 degrees of freedom
-Multiple R-squared: 0.1352, Adjusted R-squared: 0.1352 '
-
-library(scatterplot3d) # für 3D-Plots
-s3d1 <- scatterplot3d(CWM[,"SLA.mean"],CWM[,"LeafNPratio.mean"], x,
- xlab="CWM ln(Leaf N to P ratio)", ylab="CWM ln(SLA)", zlab="BIO_1",
- pch=1,cex.symbols=0.05,highlight.3d=T,type="p")
-#fit <- lm(AnzF ~ Prop.n+Prop.ba,data=data6[data6$Sp=="Qupe",])
-s3d1$plane3d(model8)
-# does not work
-### producing the plane with points
-min(CWM[,"SLA.mean"])
-max(CWM[,"SLA.mean"])
-x2 <- sort(rep(seq(0.1, 5, by=0.1),35))
-min(CWM[,"LeafNPratio.mean"])
-max(CWM[,"LeafNPratio.mean"])
-y2 <- rep(seq(1, 4.4, by=0.1),50)
-func <- function(x,y) { -432.2500 + 227.1681 * x + 135.8572*y -60.2280*x*y}
-z2 <- func(x2,y2)
-s3d1 <- scatterplot3d(CWM[,"SLA.mean"],CWM[,"LeafNPratio.mean"], x,
- xlab="CWM ln(Leaf N to P ratio)", ylab="CWM ln(SLA)", zlab="BIO_1",
- pch=1,cex.symbols=0.05,highlight.3d=T,type="p",angle=70,scale.y=0.5)
-s3d1$points3d(x2,y2,z2,col="blue", type="p", pch=20, cex=0.5)
-#
-x <- splot.bioclim$BIO_01[index11]
-model7 <- lm(CWM[,"SLA.mean"]~x*CWM[,"LeafNPratio.mean"])
-summary(model7)
-' Estimate Std. Error t value Pr(>|t|)
-(Intercept) 1.948e+00 9.377e-03 207.8 <2e-16 ***
-x 1.278e-02 7.985e-05 160.0 <2e-16 ***
-CWM[, "LeafNPratio.mean"] 3.916e-01 3.856e-03 101.6 <2e-16 ***
-x:CWM[, "LeafNPratio.mean"] -5.472e-03 3.166e-05 -172.8 <2e-16 ***
----
-Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
-
-Residual standard error: 0.4213 on 1111303 degrees of freedom
-Multiple R-squared: 0.04842, Adjusted R-squared: 0.04842 '
-min(x)
-max(x)
-x2 <- sort(rep(seq(-100, 290, by=10),35))
-min(CWM[,"LeafNPratio.mean"])
-max(CWM[,"LeafNPratio.mean"])
-y2 <- rep(seq(1, 4.4, by=0.1),40)
-func <- function(x,y) { 1.948e+00 + 1.278e-02 * x + 3.916e-01*y -5.472e-03*x*y}
-z2 <- func(x2,y2)
-s3d1 <- scatterplot3d(x,CWM[,"LeafNPratio.mean"], CWM[,"SLA.mean"],
- xlab="BIO_1", ylab="CWM ln(Leaf N to P ratio)", zlab="CWM ln(SLA)",
- pch=1,cex.symbols=0.05,highlight.3d=T,type="p",angle=50,scale.y=1)
-s3d1$points3d(x2,y2,z2,col="blue", type="p", pch=20, cex=0.5)
-# final plot
-# some BIO_01 and LeafNPratio combinations do not exist
-?cut
-x.cut <- cut(x, breaks=seq(-110, 290, by=10))
-y.cut <- cut(CWM[,"LeafNPratio.mean"],breaks=seq(0.9, 4.4, by=0.1))
-z.cut <- table(x.cut,y.cut)
-head(z.cut)
-dim(z.cut) #40 35
-z2.cut <- table(x2,y2)
-head(z2.cut)
-dim(z2.cut) #40 35
-
-#library(reshape2)
-z.cut <- melt(z.cut)
-str(z.cut) # 1400 obs. of 3 variables:
-head(z.cut)
-z2.cut <- melt(z2.cut)
-str(z2.cut) #1400 obs. of 3 variables:
-head(z2.cut)
-z2.cut$value <- 0
-z2.cut$value <- z.cut$value
-z2.cut$xy <- paste(z2.cut$x2,z2.cut$y2,sep="_")
-#z3.cut <- acast(z2.cut, x2~y2)
-#head(z3.cut)
-xy <- paste(x2,y2,sep="_")
-index20 <- match(xy,z2.cut$xy)
-values.present <- vector(mode="logical", length=length(x2))
-values.present <- z2.cut$value[index20] > 0
-table(values.present)
-'z.cut$value[z.cut$x.cut==x.cut[1] & z.cut$y.cut==y.cut[1]]
-values.present <- function(z.cut,x,y){
- z.cut$value[z.cut$x.cut==x & z.cut$y.cut==y]
-}
-x1 <- data.frame(x,x.cut, y=CWM[,"LeafNPratio.mean"],y.cut)
-head(x1)
-z1 <- apply(x1,1,FUN=values.present(z.cut, x1$x.cut, x1$y.cut))
-values.present(z.cut, x1$x.cut[1], x1$y.cut[1])
-'
-
-s3d1 <- scatterplot3d(x,CWM[,"LeafNPratio.mean"], CWM[,"SLA.mean"],
- xlab="BIO_1", ylab="CWM ln(Leaf N to P ratio)", zlab="CWM ln(SLA)",
- pch=1,cex.symbols=0.05,highlight.3d=T,type="p",angle=-10,scale.y=0.8)
-s3d1$points3d(x2[values.present],y2[values.present],z2[values.present],col="blue", type="p", pch=20, cex=0.5)
-
-
-#### Collect all CWM regressions ####
-names(splot.bioclim) # from 4 (BIO_01 to 49 T_MAR
-CWM.bioclim <- array(NA,c(46,18,5),dimnames=list(names(splot.bioclim)[4:49],dimnames(CWM)[[2]],
- c("p_lin","r2_lin","p_qua1", "p_qua2","r2_qua")))
-for (i in 1:46){
- for (j in 1:18){
- model1 <- lm(CWM[,j]~splot.bioclim[index11,i+3])
- CWM.bioclim[i,j,1] <- summary(model1)$coefficients[2,4]
- CWM.bioclim[i,j,2] <- summary(model1)$r.squared
- model2 <- lm(CWM[,j]~splot.bioclim[index11,i+3]+I(splot.bioclim[index11,i+3]^2))
- CWM.bioclim[i,j,3] <- summary(model2)$coefficients[2,4]
- CWM.bioclim[i,j,4] <- summary(model2)$coefficients[3,4]
- CWM.bioclim[i,j,5] <- summary(model2)$r.squared
- }
-}
-str(CWM.bioclim) #num [1:46, 1:18, 1:5]
-write.csv(data.frame(names(splot.bioclim)[4:49],CWM.bioclim[,,2]),file = "CWM_bioclim_r2_lin.csv", row.names = FALSE)
-write.csv(data.frame(names(splot.bioclim)[4:49],CWM.bioclim[,,5]),file = "CWM_bioclim_r2_qua.csv", row.names = FALSE)
-
-
-which(CWM.bioclim[,1,2]==max(CWM.bioclim[,1,2])) # BIO_02
-which(CWM.bioclim[,1,5]==max(CWM.bioclim[,1,5])) # BIO_02
-which(CWM.bioclim[,,2]==max(CWM.bioclim[,,2]),arr.ind = TRUE)
-# row col
-#BIO_02 2 9
-CWM.bioclim[2,9,2] #0.1297723
-which(CWM.bioclim[,,5]==max(CWM.bioclim[,,5]),arr.ind = TRUE)
-# row col
-#BIO_02 2 9
-CWM.bioclim[2,9,5] # 0.1563736
-names(splot.bioclim)[2+3] #"BIO_02"
-dimnames(CWM)[[2]][9] #"LeafNperArea.mean"
-plot(CWM[,"LeafNperArea.mean"]~splot.bioclim$BIO_02[index11],
- cex=0.05, xlab="BIO_02", ylab="CWM ln(Leaf N per Area)",cex.lab=1.5, cex.axis=1.5)
-model1 <- lm(CWM[,9]~splot.bioclim[index11,2+3])
-abline(model1, lwd=3, col="blue")
-# MeanNperArea_BIO_02.png
-png(filename="MeanNperArea_BIO_02.png", width=800, height=800, units = "px", bg="white")
- plot(CWM[,"LeafNperArea.mean"]~splot.bioclim$BIO_02[index11],
- cex=0.05, xlab="BIO_02", ylab="CWM ln(Leaf N per Area)",cex.lab=1.5, cex.axis=1.5)
- #model1 <- lm(CWM[,9]~splot.bioclim[index11,2+3])
- abline(model1, lwd=3, col="blue")
-dev.off() # creates a file
-png(filename="SLA_BIO_02.png", width=800, height=800, units = "px", bg="white")
- plot(CWM[,"SLA.mean"]~splot.bioclim$BIO_02[index11],
- cex=0.05, xlab="BIO_02", ylab="CWM ln(SLA)",cex.lab=1.5, cex.axis=1.5)
- model1 <- lm(CWM[,1]~splot.bioclim[index11,2+3])
-abline(model1, lwd=3, col="blue")
-dev.off() # creates a file
-
-#### PCA of all CWM and all Bioclim variables ####
-library(vegan)
-pca1 <- rda(CWM,scale=T)
-plot(pca1)
-str(summary(pca1))
-#text(summary(pca1)$species[,c(1,2)],dimnames(summary(pca1)$species)[[1]])
-traitnames <- dimnames(summary(pca1)$species)[[1]]
-traitnames <- substr(traitnames,1,nchar(traitnames)-5)
-traitcoord <- summary(pca1)$species[,c(1,2)]
-which(traitnames=="SeedMass") #3
-traitcoord[3,2] <- traitcoord[3,2]-0.5
-ordiplot(pca1, type="n", ylim=c(-10,10))
-text(traitcoord,traitnames, col="red", cex=1)
-points(plotcoord*30, cex=0.05)
-#env1 <- envfit(pca1,splot.bioclim[index11,c(4:49)],permutations = 99)
-#with permutations = 999:
-#Error: cannot allocate vector of size 4.1 Gb
-plot(env1)
-png(filename="PCA1.png", width=1000, height=1000, units = "px", bg="white")
- ordiplot(pca1, type="n")
- text(traitcoord,traitnames)
- plot(env1)
-dev.off() # creates a file
-'rm(env1)
-rm(backbone.splot.try3)
-rm(model1)
-rm(model2)
-rm(pca1)
-rm(splot.header)
-rm(splot.header2)
-rm(traitdist)
-rm(coords)
-rm(TRY.all.mean.sd.2)
-rm(index7)
-rm(raster.map.SLA)
-rm(trait.vec)
-rm(FD)
-rm(splot.bioclim)
-'
-gc()
-
-### produce a PCA with kernel density
-# From Díaz et al. 2016, Nature
-library(ks)
-plotcoord <- summary(pca1)$sites[,c(1,2)]
-str(plotcoord) #num [1:1111307, 1:2]
-H <- Hpi(x=plotcoord) # optimal bandwidth estimation
-# does not work
-est<- kde(x=plotcoord, H=H, compute.cont=TRUE) # kernel density estimation
-
-library(raster)
-library(fBasics)
-# Color palette
-#plotcoord <- plotcoord*30
-min(plotcoord[,1]) #-0.2639705
-max(plotcoord[,1]) # 0.2415964
-min(plotcoord[,2]) #-0.2595053
-max(plotcoord[,2]) # 0.3560549
-x.cut <- cut(plotcoord[,1], breaks=seq(min(plotcoord[,1]), max(plotcoord[,1]), by=0.01))
-y.cut <- cut(plotcoord[,2], breaks=seq(min(plotcoord[,2]), max(plotcoord[,2]), by=0.01))
-z.cut <- table(x.cut,y.cut)
-head(z.cut)
-dim(z.cut) #50 61
-z.cut <- melt(z.cut)
-str(z.cut) # 3050 obs. of 3 variables:
-head(z.cut)
-z.cut$xy <- paste(z.cut$x.cut,z.cut$y.cut,sep="_")
-str(x.cut)
-length(x.cut) # 1111307
-xy <- paste(x.cut,y.cut,sep="_")
-index30 <- match(xy,z.cut$xy)
-length(index30)
-any(is.na(index30))
-head(log10(z.cut$value[index30]))
-values.present <- vector(mode="numeric", length=length(xy))
-min(log10(z.cut$value[index30]), na.rm=T) #0
-max(log10(z.cut$value[index30]), na.rm=T) # 3.9878
-values.cut <- cut(log10(z.cut$value[index30]), breaks=seq(min(log10(z.cut$value[index30]), na.rm=T),
- 4, by=0.4))
-#values.present <- z.cut$value[index30]
-ordiplot(pca1, type="n")
-#points(plotcoord*30, cex=0.05, col=rainbowPalette(n=10,
-# s = 1, v = 1, start = 3/6, end = 5/6, alpha = 1)[values.cut])
-points(plotcoord*30, cex=0.05, col=rev(greyPalette(n=10, start = 0, end = 0.8,
- gamma = 2.2, alpha = NULL))[values.cut])
-text(traitcoord,traitnames, col="red", cex=1)
-
-
-ext.factor <- 30
-r.pca <- raster(nrows=300, ncols=300, xmn=min(plotcoord[,1])*ext.factor, xmx=max(plotcoord[,2])*ext.factor,
- ymn=min(plotcoord[,1])*ext.factor, ymx=max(plotcoord[,2]*ext.factor))
-# attention, xmn and ymn as well as xmx and ymx have to be the same, otherwise the raster is shifted
-r.pca.raster <- rasterize(plotcoord*ext.factor, r.pca, fun="count")
-str(r.pca.raster)
-plot(log10(r.pca.raster), asp=0, col=seqPalette(n=10, name=c("Blues")),
- main="", cex.main=0.7, xlab="", ylab="", add=T)
-plot(log10(r.pca.raster), asp=0, col=rev(greyPalette(n=10, start = 0, end = 0.8,
- gamma = 2.2, alpha = NULL)),
- main="", cex.main=0.7, xlab="", ylab="", add=T)
-plot(log10(r.pca.raster), asp=0, col=rainbowPalette(n=10, s = 1, v = 1, start = 3/6, end = 5/6, alpha = 1),
- main="", cex.main=0.7, xlab="", ylab="", add=T, legend=F)
-
-'## graph for the PCA ##
-ordiplot(pca1, type="n", xlim=c(-10,13), ylim=c(-10,10))
-rasterImage(as.raster(r.pca.raster), -10, -10, 20, 10, interpolate = F,
- col=rainbowPalette(n=10, s = 1, v = 1, start = 3/6, end = 5/6, alpha = 1))
-plot(log10(r.pca.raster), asp=0, col=rainbowPalette(n=10, s = 1, v = 1, start = 3/6, end = 5/6, alpha = 1),
- main="", cex.main=0.7, xlab="", ylab="", add=T, legend=F)
-text(traitcoord,traitnames, col="red", cex=1)
-'
-
-
-rda1 <- rda(CWM~.,splot.bioclim[index11,c(4:49)],scale=T)
-plot(rda1)
-summary(rda1)
-rda1
-' Inertia Proportion Rank
-Total 18.0000 1.0000
-Constrained 2.4811 0.1378 15
-Unconstrained 15.5189 0.8622 15
-Inertia is correlations
-Eigenvalues for constrained axes:
- RDA1 RDA2 RDA3 RDA4 RDA5 RDA6 RDA7 RDA8 RDA9 RDA10 RDA11 RDA12 RDA13 RDA14 RDA15
-1.2230 0.6020 0.3275 0.1972 0.0653 0.0226 0.0156 0.0107 0.0059 0.0050 0.0031 0.0012 0.0010 0.0007 0.0003
-
-Eigenvalues for unconstrained axes:
-PC1 PC2 PC3 PC4 PC5 PC6 PC7 PC8 PC9 PC10 PC11 PC12 PC13 PC14 PC15
-4.594 3.453 1.486 1.426 1.095 0.798 0.648 0.603 0.483 0.343 0.261 0.122 0.084 0.074 0.047
-'
-#anova(rda1, perm=99)
-#anova(rda1, by="terms", permutations = how(nperm=9))
-library(parallel)
-detectCores() #32
-# Set up the cluster
-#clusterType <- if(length(find.package("snow", quiet = TRUE))) "SOCK" else "PSOCK"
-#clust <- try(makeCluster(getOption("cl.cores", 2), type = clusterType))
-cl <- makeCluster(10)
-#registerDoParallel(cl)
-options(na.action = "na.pass")
-#clusterExport(cl, "results")
-
-#print(system.time(anova(rda1, by="terms", perm=99, parallel = getOption("mc.cores"))))
-print(system.time(out1 <- anova(rda1, by="margin", permutations =99, parallel = cl)))
-'user system elapsed
-5304.802 1311.088 14522.410
-'
-out1
-
-traitnames <- dimnames(summary(rda1)$species)[[1]]
-traitnames <- substr(traitnames,1,nchar(traitnames)-5)
-traitcoord <- summary(rda1)$species[,c(1,2)]
-bioclimcoord <- summary(rda1)$biplot[,c(1,2)]
-which(traitnames=="LeafNperArea") #9
-traitcoord[9,1] <- traitcoord[9,1]+0.5
-which(traitnames=="StemDens") #5
-traitcoord[5,2] <- traitcoord[5,2]+0.3
-which(traitnames=="Disp.unit.leng") #18
-traitcoord[18,2] <- traitcoord[18,2]-0.3
-ordiplot(rda1, type="n", xlim=c(-7,7))
-arrows(rep(0,45),rep(0,45),bioclimcoord[,1]*8,bioclimcoord[,2]*8, col="cadetblue4", length=0.1)
-text(bioclimcoord*8,dimnames(bioclimcoord)[[1]], col="cadetblue4", cex=0.8)
-text(traitcoord,traitnames, col="red")
-# copied as first RDA graph into ppt
-
-#rda1 <- rda(CWM ~ ., splot.bioclim[index11,c(4:49)],scale=T)
-
-### Stepwise RDA ###
-rda2 <- rda(CWM ~ 1,splot.bioclim[index11,c(4:49)],scale=T)
-#rda3 <- step(rda2, scope = formula(rda1), test = "perm")
-rda3 <- ordiR2step(rda2, scope=formula(rda1), direction = "both", Pin = 0.05, R2scope = F, permutations = how(nperm = 99), trace = TRUE)
-#CWM ~ GDD5 + BIO_02 + BIO_12 + T_MAY + P_MAY + T_AUG ....
-#rm(rda1)
-
-# now rebuild the RDA sequentially
-pca1 <- rda(CWM,scale=T)
-pca1
-5.527+ 3.906+ 1.953+ 1.507+ 1.375+ 0.886+ 0.726+ 0.623+ 0.499 +0.360+ 0.295+ 0.131+ 0.087+ 0.077+ 0.048
-(5.527+ 3.906)/18 # 52.40556%
-
-traitnames <- dimnames(summary(pca1)$species)[[1]]
-traitnames <- substr(traitnames,1,nchar(traitnames)-5)
-traitcoord <- summary(pca1)$species[,c(1,2)]
-which(traitnames=="SeedMass") #3
-traitcoord[3,2] <- traitcoord[3,2]-0.5
-ordiplot(pca1, type="n")
-text(traitcoord,traitnames, col="red")
-# copied as second RDA graph into ppt
-
-
-rda4 <- rda(CWM ~ GDD5,splot.bioclim[index11,c(4:49)],scale=T)
-rda4
-traitnames <- dimnames(summary(rda4)$species)[[1]]
-traitnames <- substr(traitnames,1,nchar(traitnames)-5)
-traitcoord <- summary(rda4)$species[,c(1,2)]
-plotcoord <- summary(rda4)$sites[,c(1,2)]
-str(plotcoord) #num [1:1111307, 1:2]
-bioclimcoord <- summary(rda4)$biplot[,c(1,2)]
-
-min(plotcoord[,1]) #-0.5340595
-max(plotcoord[,1]) # 0.7760822
-min(plotcoord[,2]) #-0.2688765
-max(plotcoord[,2]) # 0.2326464
-x.cut <- cut(plotcoord[,1], breaks=seq(min(plotcoord[,1]), max(plotcoord[,1]), by=0.01))
-y.cut <- cut(plotcoord[,2], breaks=seq(min(plotcoord[,2]), max(plotcoord[,2]), by=0.01))
-z.cut <- table(x.cut,y.cut)
-head(z.cut)
-dim(z.cut) #131 50
-z.cut <- melt(z.cut)
-str(z.cut) # 6550 obs. of 3 variables:
-head(z.cut)
-z.cut$xy <- paste(z.cut$x.cut,z.cut$y.cut,sep="_")
-str(x.cut)
-length(x.cut) # 1111307
-xy <- paste(x.cut,y.cut,sep="_")
-index30 <- match(xy,z.cut$xy)
-length(index30)
-any(is.na(index30))
-head(log10(z.cut$value[index30]))
-values.present <- vector(mode="numeric", length=length(xy))
-min(log10(z.cut$value[index30]), na.rm=T) #0
-max(log10(z.cut$value[index30]), na.rm=T) # 3.956168
-values.cut <- cut(log10(z.cut$value[index30]), breaks=seq(min(log10(z.cut$value[index30]), na.rm=T),
- 4, by=0.4))
-
-
-ordiplot(rda4, type="n")
-points(plotcoord*30, cex=0.05, col=rev(greyPalette(n=10, start = 0, end = 0.8,
- gamma = 2.2, alpha = NULL))[values.cut])
-
-arrows(rep(0,1),rep(0,1),bioclimcoord[1]*8,bioclimcoord[2]*8, col="blue", length=0.1)
-text(bioclimcoord[1]*8+0.8,bioclimcoord[2]*8-0.2,"GDD5", col="blue", cex=1)
-text(traitcoord,traitnames, col="red")
-# copied as third RDA graph into ppt
-
-rda5 <- rda(CWM ~ GDD5+BIO_02,splot.bioclim[index11,c(4:49)],scale=T)
-rda5
-traitnames <- dimnames(summary(rda5)$species)[[1]]
-traitnames <- substr(traitnames,1,nchar(traitnames)-5)
-traitcoord <- summary(rda5)$species[,c(1,2)]
-bioclimcoord <- summary(rda5)$biplot[,c(1,2)]
-plotcoord <- summary(rda5)$sites[,c(1,2)]
-str(plotcoord) #num [1:1111307, 1:2]
-x.cut <- cut(plotcoord[,1], breaks=seq(min(plotcoord[,1]), max(plotcoord[,1]), by=0.01))
-y.cut <- cut(plotcoord[,2], breaks=seq(min(plotcoord[,2]), max(plotcoord[,2]), by=0.01))
-z.cut <- table(x.cut,y.cut)
-head(z.cut)
-dim(z.cut) #131 50
-z.cut <- melt(z.cut)
-str(z.cut) # 6550 obs. of 3 variables:
-head(z.cut)
-z.cut$xy <- paste(z.cut$x.cut,z.cut$y.cut,sep="_")
-str(x.cut)
-length(x.cut) # 1111307
-xy <- paste(x.cut,y.cut,sep="_")
-index30 <- match(xy,z.cut$xy)
-length(index30)
-any(is.na(index30))
-head(log10(z.cut$value[index30]))
-values.present <- vector(mode="numeric", length=length(xy))
-min(log10(z.cut$value[index30]), na.rm=T) #0
-max(log10(z.cut$value[index30]), na.rm=T) # 3.399154
-values.cut <- cut(log10(z.cut$value[index30]), breaks=seq(min(log10(z.cut$value[index30]), na.rm=T),
- 4, by=0.4))
-
-
-
-ordiplot(rda5, type="n", xlim=c(-10,10), ylim=c(-10,10))
-points(plotcoord*30, cex=0.05, col=rev(greyPalette(n=10, start = 0, end = 0.8,
- gamma = 2.2, alpha = NULL))[values.cut])
-arrows(rep(0,2),rep(0,2),bioclimcoord[,1]*8,bioclimcoord[,2]*8, col="blue", length=0.1)
-text(bioclimcoord*8,dimnames(bioclimcoord)[[1]], col="blue", cex=1)
-text(traitcoord,traitnames, col="red")
-# copied as fourth RDA graph into ppt
-
-
-plot(CWM[,"LeafNPratio.mean"]~splot.bioclim$GDD5[index11],
- cex=0.05, xlab="GDD5", ylab="CWM ln(leaf N to P ratio)",cex.lab=1.5, cex.axis=1.5)
-model1 <- lm(CWM[,"LeafNPratio.mean"]~splot.bioclim$GDD5[index11])
-summary(model1)
-' Estimate Std. Error t value Pr(>|t|)
-(Intercept) 2.297e+00 3.891e-04 5905.2 <2e-16 ***
-splot.bioclim$GDD5[index11] 6.349e-06 1.636e-08 388.1 <2e-16 ***
-Multiple R-squared: 0.1193, Adjusted R-squared: 0.1193
-'
-abline(model1, lwd=3, col="blue")
-# LeafNPRatio_GDD5
-
-plot(CWM[,"SLA.mean"]~jitter(splot.bioclim$BIO_02[index11]),
- cex=0.05, xlab="BIO_02", ylab="CWM ln(SLA)",cex.lab=1.5, cex.axis=1.5)
-model1 <- lm(CWM[,"SLA.mean"]~splot.bioclim$BIO_02[index11])
-summary(model1)
-' Estimate Std. Error t value Pr(>|t|)
-(Intercept) 3.399e+00 1.448e-03 2348.2 <2e-16 ***
-splot.bioclim$BIO_02[index11] -6.253e-03 1.551e-05 -403.3 <2e-16 ***
-
-Residual standard error: 0.4034 on 1111305 degrees of freedom
-Multiple R-squared: 0.1277, Adjusted R-squared: 0.1277
-'
-abline(model1, lwd=3, col="blue")
-# SLA_BIO_02
-
-plot(CWM[,"LeafNperArea.mean"]~jitter(splot.bioclim$BIO_02[index11]),
- cex=0.05, xlab="BIO_02", ylab="CWM ln(Leaf N per area)",cex.lab=1.5, cex.axis=1.5)
-model1 <- lm(CWM[,"LeafNperArea.mean"]~splot.bioclim$BIO_02[index11])
-summary(model1)
-' Estimate Std. Error t value Pr(>|t|)
-(Intercept) -5.989e-02 8.475e-04 -70.67 <2e-16 ***
-splot.bioclim$BIO_02[index11] 3.696e-03 9.079e-06 407.09 <2e-16 ***
----
-Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
-
-Residual standard error: 0.2362 on 1111305 degrees of freedom
-Multiple R-squared: 0.1298, Adjusted R-squared: 0.1298
-'
-abline(model1, lwd=3, col="blue")
-# LeafNper_BIO_02
-
-##### FD calculations ####
-j <- 3
-nam <- DT$species[DT$PlotObservationID==j]
-abu <- DT$Relative.cover[DT$PlotObservationID==j]
-trait <- DT$trait[DT$PlotObservationID==j]
-abundance.weighted=T
-
-FD.fun <- function(trait, abu){
- res <- as.double(NA)
- if (length(trait[!is.na(trait)])>0){
- res <- 0
- #nam <- as.character(nam[!is.na(trait)])
- abu <- as.data.frame(abu[!is.na(trait)])
- #rownames(abu) <- nam
- trait <- as.data.frame(trait[!is.na(trait)])
- #rownames(trait) <- nam
- dis <- dist(trait,method="euclidean")
- if (sum(dis)>1){
- sqrt.dis <- sqrt(dis)
- # taking the square root is required to obtain the same result from the divc function
- # as divc takes the square of distances values, as suggested by
- # Champely & Chessel (2002, Env. Ecol. Stat. 9: 167-177)
- #abu <- as.data.frame(abu)
- res <- unlist(divc(abu, sqrt.dis))
- #Rao's Q: sumi sumj (dist ij prop i prop j)'
- }
- }
- res
-}
-
-mpd.fun <- function(nam, trait, abu,abundance.weighted){
- #res <- numeric(1) # carries result from the function
- res <- as.double(NA)
- if (length(trait[!is.na(trait)])>0){
- # dis <- as.matrix(dist(t(t1),method="euclidean"))
- abu <- t(abu)
- colnames(abu) <- nam
- trait <- t(trait)
- colnames(trait) <- nam
- #dis <- as.matrix(vegdist(t(t1),method="euclidean", na.rm=T))
- dis <- as.matrix(dist(t(trait),method="euclidean"))
- res <- mpd(abu, dis, abundance.weighted=abundance.weighted)
- }
- res
-}
-
-# from mpd source code
-'N <- dim(samp)[1]
-mpd <- numeric(N)
-for (i in 1:N) {
- sppInSample <- names(samp[i, samp[i, ] > 0])
- if (length(sppInSample) > 1) {
- sample.dis <- dis[sppInSample, sppInSample]
- if (abundance.weighted) {
- sample.weights <- t(as.matrix(samp[i, sppInSample,
- drop = FALSE])) %*% as.matrix(samp[i, sppInSample,
- drop = FALSE])
- mpd[i] <- weighted.mean(sample.dis, sample.weights)
- }
- else {
- mpd[i] <- mean(sample.dis[lower.tri(sample.dis)])
- }
- }
- else {
- mpd[i] <- NA
- }
-}
-mpd
-'
-# from Oliver, who used picante
-# does not accept missing values in the function
-library(picante)
-# define function for mean pairwise distance
-'mpd.fun.dt <- function(abu, dis, nam, abundance.weighted){
- abu <- t(abu)
- colnames(abu) <- nam
- res <- mpd(abu, dis, abundance.weighted=abundance.weighted)
-}
-# define funtion of mean nearest neighbor distance
-mntd.fun.dt <- function(abu, dis, nam, abundance.weighted){
- abu <- t(abu)
- colnames(abu) <- nam
- res <- mntd(abu, dis, abundance.weighted=abundance.weighted)
- res
-}
-species.list <- unique(DT$species)
-any(is.na(species.list)) #T, there is one NA
-length(species.list) #60909
-species.list <- species.list[!is.na(species.list)]
-length(species.list) #60908
-index8 <- match(species.list,TRY.all.mean.sd.2$species)
-str(index8) # 60908
-length(index8[!is.na(index8)]) #21001
-#species.list2 <- TRY.all.mean.sd.2$species[index8]
-#species.list2 <- species.list2[!is.na(species.list2)]
-#length(species.list2) # 21001
-TRY.all.mean.sd.2.na <- TRY.all.mean.sd.2[index8,]
-str(TRY.all.mean.sd.2.na)
-TRY.all.mean.sd.2.na$species2 <- species.list
-rm(index9)
-'
-FD <- array(NA,c(dim(CWM)[1],18),dimnames=list(dimnames(CWM)[[1]],
- colnames(TRY.all.mean.sd.2)[3:20]))
-str(FD) #logi [1:1117898, 1:18]
-'for (i in 1:18){
- # create distance matrix
- #trait.vec <- TRY.all.mean.sd.2.na[,i+2]
- #names(trait.vec) <- TRY.all.mean.sd.2.na$species2
- #trait.dist.mat <- as.matrix(dist(trait.vec))
- # too large for my notebook
- #mpd[,i] <- DT[,list(mpd.trait= mpd.fun.dt(Relative.cover, trait.dist.mat, species, abundance.weighted=TRUE)),by=PlotObservationID]$mpd.trait
- print(i)
- DT$trait <- NA
- DT$trait <- TRY.all.mean.sd.2[index7,i+2]
- FD[,i] <- DT[,list(mpd.trait= mpd.fun(species, trait, Relative.cover, abundance.weighted=TRUE)),by=PlotObservationID]$mpd.trait
-}
-write.csv(FD,file = "FD2.csv", row.names = T)
-'
-for (i in 1:18){
- # create distance matrix
- #trait.vec <- TRY.all.mean.sd.2.na[,i+2]
- #names(trait.vec) <- TRY.all.mean.sd.2.na$species2
- #trait.dist.mat <- as.matrix(dist(trait.vec))
- # too large for my notebook
- #mpd[,i] <- DT[,list(mpd.trait= mpd.fun.dt(Relative.cover, trait.dist.mat, species, abundance.weighted=TRUE)),by=PlotObservationID]$mpd.trait
- print(i)
- DT$trait <- NA
- DT$trait <- TRY.all.mean.sd.2[index7,i+2]
- # recalculate relative cover for this trait
- FD[,i] <- DT[,list(FD.trait= FD.fun(trait, Relative.cover)),by=PlotObservationID]$FD.trait
-}
-write.csv(FD,file = "FD3.csv", row.names = T)
-
-FD <- read.csv("C://Daten//iDiv2//splot2//FD3.csv")
-str(FD)
-dimnames(FD)[[1]] <- FD$X
-FD <- FD[,-1]
-tail(FD)
-
-#### match FD to splot.header ####
-index12 <- match(dimnames(FD)[[1]],splot.header$PlotObservationID)
-length(index12) #1117898
-any(is.na(index12)) # F
-splot.header3 <- splot.header[index12,]
-str(FD) # num [1:1117898, 1:18]
-FD <- FD[!(is.na(splot.header3$Longitude)|is.na(splot.header3$Latitude)|is.na(FD[,"SLA.mean"])),]
-str(FD) # num [1:1111307 , 1:18]
-1111307 /1117898*100 # 99.41041 % of plots have traits, Long and Lat
-any(is.na(FD)) #F
-# match header again
-index13 <- match(dimnames(FD)[[1]],splot.header$PlotObservationID)
-length(index13) #1111307
-any(is.na(index13)) # F
-splot.header3 <- splot.header[index13,]
-# now in the same sequence as FD
-
-
-#### Collect all FD regressions ####
-names(splot.bioclim) # from 4 BIO_01 to 49 T_MAR
-index14 <- match(splot.header3$PlotObservationID,splot.bioclim$PLOT_ID)
-length(index14) # 1111307
-length(index14[!is.na(index14)]) # 1111307
-
-FD.bioclim <- array(NA,c(46,18,5),dimnames=list(names(splot.bioclim)[4:49],dimnames(FD)[[2]],
- c("p_lin","r2_lin","p_qua1", "p_qua2","r2_qua")))
-for (i in 1:46){
- for (j in 1:18){
- model1 <- lm(FD[,j]~splot.bioclim[index14,i+3])
- FD.bioclim[i,j,1] <- summary(model1)$coefficients[2,4]
- FD.bioclim[i,j,2] <- summary(model1)$r.squared
- model2 <- lm(FD[,j]~splot.bioclim[index14,i+3]+I(splot.bioclim[index14,i+3]^2))
- FD.bioclim[i,j,3] <- summary(model2)$coefficients[2,4]
- FD.bioclim[i,j,4] <- summary(model2)$coefficients[3,4]
- FD.bioclim[i,j,5] <- summary(model2)$r.squared
- }
-}
-str(FD.bioclim) #num [1:46, 1:18, 1:5]
-write.csv(data.frame(names(splot.bioclim)[4:49],FD.bioclim[,,2]),file = "FD2_bioclim_r2_lin.csv", row.names = FALSE)
-write.csv(data.frame(names(splot.bioclim)[4:49],FD.bioclim[,,5]),file = "FD2_bioclim_r2_qua.csv", row.names = FALSE)
-
-library(vegan)
-FD.pca1 <- rda(FD,scale=T)
-FD.pca1
-'Eigenvalues for unconstrained axes:
-PC1 PC2 PC3 PC4 PC5 PC6 PC7 PC8
-8.597 1.876 1.378 0.886 0.790 0.780 0.583 0.551
-(Showed only 8 of all 18 unconstrained eigenvalues)
-'
-plot(FD.pca1)
-summary(FD.pca1)
-#text(summary(FD.pca1)$species[,c(1,2)],dimnames(summary(FD.pca1)$species)[[1]])
-traitnames <- dimnames(summary(FD.pca1)$species)[[1]]
-traitnames <- substr(traitnames,1,nchar(traitnames)-5)
-traitcoord <- summary(FD.pca1)$species[,c(1,2)]
-which(traitnames=="Seed.length") #17
-traitcoord[17,2] <- traitcoord[17,2]+0.3
-which(traitnames=="LeafC.perdrymass") #12
-traitcoord[12,2] <- traitcoord[12,2]-0.4
-which(traitnames=="Leaf.delta.15N") #13
-traitcoord[13,2] <- traitcoord[13,2]+0.4
-which(traitnames=="LeafNPratio") #11
-traitcoord[11,2] <- traitcoord[11,2]-0.3
-
-plotcoord <- summary(FD.pca1)$sites[,c(1,2)]
-str(plotcoord) #num [1:1111307, 1:2]
-x.cut <- cut(plotcoord[,1], breaks=seq(min(plotcoord[,1]), max(plotcoord[,1]), by=0.01))
-y.cut <- cut(plotcoord[,2], breaks=seq(min(plotcoord[,2]), max(plotcoord[,2]), by=0.01))
-z.cut <- table(x.cut,y.cut)
-head(z.cut)
-dim(z.cut) #131 50
-z.cut <- melt(z.cut)
-str(z.cut) # 6550 obs. of 3 variables:
-head(z.cut)
-z.cut$xy <- paste(z.cut$x.cut,z.cut$y.cut,sep="_")
-str(x.cut)
-length(x.cut) # 1111307
-xy <- paste(x.cut,y.cut,sep="_")
-index30 <- match(xy,z.cut$xy)
-length(index30)
-any(is.na(index30))
-head(log10(z.cut$value[index30]))
-values.present <- vector(mode="numeric", length=length(xy))
-min(log10(z.cut$value[index30]), na.rm=T) #0
-max(log10(z.cut$value[index30]), na.rm=T) # 3.399154
-values.cut <- cut(log10(z.cut$value[index30]), breaks=seq(min(log10(z.cut$value[index30]), na.rm=T),
- 4, by=0.4))
-
-ordiplot(FD.pca1, type="n")
-#text(traitcoord,traitnames)
-points(plotcoord*30, cex=0.05, col=rev(greyPalette(n=10, start = 0, end = 0.8,
- gamma = 2.2, alpha = NULL))[values.cut])
-
-#arrows(rep(0,45),rep(0,45),bioclimcoord[,1]*8,bioclimcoord[,2]*8, col="cadetblue4", length=0.1)
-text(traitcoord,traitnames, col="red")
-# copied as first RDA graph into ppt
-
-env1 <- envfit(FD.pca1,splot.bioclim[index14,c(4:49)],permutations = 0)
-env1
-#with permutations = 999:
-#Error: cannot allocate vector of size 4.1 Gb
-plot(env1)
-png(filename="PCA1.png", width=1000, height=1000, units = "px", bg="white")
-ordiplot(pca1, type="n")
-text(traitcoord,traitnames)
-plot(env1)
-dev.off() # creates a file
-
-### Stepwise RDA ###
-### RDA ###
-FD.rda1 <- rda(FD~.,splot.bioclim[index14,c(4:49)],scale=T)
-FD.rda1
-FD.rda2 <- rda(FD ~ 1,splot.bioclim[index14,c(4:49)],scale=T)
-FD.rda3 <- ordistep(FD.rda2, scope=formula(FD.rda1), direction = "both", Pin = 0.1,
- permutations = how(nperm = 10), trace = TRUE)
-'Start: FD ~ 1
-
-Df AIC F Pr(>F)
-+ BIO_14 1 3191679 20602.8 0.09091 .
-+ BIO_17 1 3191898 20379.5 0.09091 .
-+ BIO_12 1 3197070 15125.2 0.09091 .
-+ P_MAY 1 3198071 14110.8 0.09091 .
-+ P_APR 1 3198194 13986.4 0.09091 .
-+ BIO_18 1 3198767 13406.3 0.09091 .
-+ P_MAR 1 3199448 12717.1 0.09091 .
-+ P_FEB 1 3199451 12714.5 0.09091 .
-+ P_JUN 1 3199643 12519.9 0.09091 .
-+ BIO_05 1 3199804 12357.7 0.09091 .
-+ BIO_02 1 3200062 12096.5 0.09091 .
-+ P_NOV 1 3200602 11551.1 0.09091 .
-+ P_DEC 1 3200823 11328.1 0.09091 .
-+ P_JAN 1 3201335 10810.2 0.09091 .
-+ BIO_10 1 3202194 9943.1 0.09091 .
-+ T_JUL 1 3202226 9911.4 0.09091 .
-+ T_AUG 1 3202685 9448.0 0.09091 .
-+ GDD5 1 3202914 9217.7 0.09091 .
-+ T_JUN 1 3202934 9196.7 0.09091 .
-+ BIO_16 1 3203459 8668.3 0.09091 .
-+ BIO_19 1 3203685 8440.1 0.09091 .
-+ P_AUG 1 3204079 8043.8 0.09091 .
-+ BIO_13 1 3204383 7737.1 0.09091 .
-+ T_MAY 1 3204450 7669.7 0.09091 .
-+ P_OCT 1 3204755 7362.4 0.09091 .
-+ T_SEP 1 3204986 7130.0 0.09091 .
-+ P_SEP 1 3205174 6941.4 0.09091 .
-+ BIO_15 1 3205378 6735.8 0.09091 .
-+ P_JUL 1 3205600 6512.4 0.09091 .
-+ T_APR 1 3205746 6365.1 0.09091 .
-+ BIO_01 1 3205841 6270.3 0.09091 .
-+ BIO_03 1 3205882 6228.8 0.09091 .
-+ T_OCT 1 3206048 6062.0 0.09091 .
-+ AR 1 3206562 5545.3 0.09091 .
-+ BIO_07 1 3207063 5042.2 0.09091 .
-+ T_NOV 1 3207547 4555.5 0.09091 .
-+ T_MAR 1 3207659 4443.2 0.09091 .
-+ T_FEB 1 3208177 3923.3 0.09091 .
-+ T_DEC 1 3208260 3840.3 0.09091 .
-+ BIO_11 1 3208334 3765.7 0.09091 .
-+ T_JAN 1 3208465 3634.3 0.09091 .
-+ BIO_08 1 3208844 3254.4 0.09091 .
-+ BIO_06 1 3209374 2723.0 0.09091 .
-+ BIO_09 1 3209543 2553.3 0.09091 .
-+ PET 1 3209642 2454.5 0.09091 .
-+ BIO_04 1 3210074 2021.3 0.09091 .
----
-Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
-
-Step: FD ~ BIO_14
-
-Df AIC F Pr(>F)
-- BIO_14 1 3212091 20603 0.09091 .
----
-Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
-
-Df AIC F Pr(>F)
-+ BIO_02 1 3184224 7481.80 0.09091 .
-+ BIO_03 1 3184435 7269.37 0.09091 .
-+ GDD5 1 3187255 4435.29 0.09091 .
-+ BIO_05 1 3187772 3915.70 0.09091 .
-+ T_JAN 1 3188026 3661.00 0.09091 .
-+ BIO_01 1 3188040 3647.18 0.09091 .
-+ T_FEB 1 3188048 3638.72 0.09091 .
-+ T_MAR 1 3188075 3612.29 0.09091 .
-+ T_APR 1 3188082 3605.14 0.09091 .
-+ BIO_11 1 3188088 3598.96 0.09091 .
-+ T_DEC 1 3188123 3563.70 0.09091 .
-+ T_NOV 1 3188208 3478.25 0.09091 .
-+ T_OCT 1 3188225 3461.84 0.09091 .
-+ T_MAY 1 3188357 3329.42 0.09091 .
-+ BIO_06 1 3188491 3194.91 0.09091 .
-+ BIO_10 1 3188609 3076.07 0.09091 .
-+ T_SEP 1 3188617 3068.49 0.09091 .
-+ BIO_13 1 3188664 3021.74 0.09091 .
-+ T_JUN 1 3188669 3016.26 0.09091 .
-+ BIO_16 1 3188709 2976.11 0.09091 .
-+ P_SEP 1 3188724 2961.09 0.09091 .
-+ BIO_07 1 3188728 2956.84 0.09091 .
-+ BIO_12 1 3188782 2902.48 0.09091 .
-+ T_AUG 1 3188788 2897.30 0.09091 .
-+ T_JUL 1 3188879 2805.50 0.09091 .
-+ P_AUG 1 3189181 2503.15 0.09091 .
-+ BIO_18 1 3189256 2427.84 0.09091 .
-+ BIO_04 1 3189277 2406.93 0.09091 .
-+ P_JUL 1 3189398 2285.76 0.09091 .
-+ BIO_15 1 3189605 2078.03 0.09091 .
-+ P_APR 1 3189650 2032.61 0.09091 .
-+ P_JUN 1 3189703 1979.86 0.09091 .
-+ BIO_09 1 3189912 1770.71 0.09091 .
-+ P_MAR 1 3189965 1717.19 0.09091 .
-+ PET 1 3189984 1698.34 0.09091 .
-+ AR 1 3190003 1678.99 0.09091 .
-+ P_MAY 1 3190005 1677.55 0.09091 .
-+ P_OCT 1 3190063 1619.65 0.09091 .
-+ P_FEB 1 3190490 1192.28 0.09091 .
-+ BIO_17 1 3190497 1184.39 0.09091 .
-+ BIO_19 1 3190555 1126.59 0.09091 .
-+ P_NOV 1 3190570 1111.68 0.09091 .
-+ BIO_08 1 3190582 1100.05 0.09091 .
-+ P_JAN 1 3190893 787.98 0.09091 .
-+ P_DEC 1 3190955 726.52 0.09091 .
----
-Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
-
-Step: FD ~ BIO_14 + BIO_02
-
-Df AIC F Pr(>F)
-- BIO_02 1 3191679 7481.8 0.09091 .
-- BIO_14 1 3200062 15953.2 0.09091 .
----
-Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1'
-
-FD.rda4 <- rda(FD ~ BIO_14,splot.bioclim[index14,c(4:49)],scale=T)
-FD.rda4
-' Inertia Proportion Rank
-Total 18.0000 1.0000
-Constrained 0.3276 0.0182 1
-Unconstrained 17.6724 0.9818 18
-Inertia is correlations
-
-Eigenvalues for constrained axes:
-RDA1
-0.3276
-
-Eigenvalues for unconstrained axes:
-PC1 PC2 PC3 PC4 PC5 PC6 PC7 PC8
-8.307 1.850 1.378 0.884 0.789 0.780 0.578 0.551
-'
-traitnames <- dimnames(summary(FD.rda4)$species)[[1]]
-traitnames <- substr(traitnames,1,nchar(traitnames)-5)
-traitcoord <- summary(FD.rda4)$species[,c(1,2)]
-bioclimcoord <- summary(FD.rda4)$biplot[,c(1,2)]
-plotcoord <- summary(FD.rda4)$sites[,c(1,2)]
-str(plotcoord) #num [1:1111307, 1:2]
-x.cut <- cut(plotcoord[,1], breaks=seq(min(plotcoord[,1]), max(plotcoord[,1]), by=0.01))
-y.cut <- cut(plotcoord[,2], breaks=seq(min(plotcoord[,2]), max(plotcoord[,2]), by=0.01))
-z.cut <- table(x.cut,y.cut)
-head(z.cut)
-dim(z.cut) #131 50
-z.cut <- melt(z.cut)
-str(z.cut) # 6550 obs. of 3 variables:
-head(z.cut)
-z.cut$xy <- paste(z.cut$x.cut,z.cut$y.cut,sep="_")
-str(x.cut)
-length(x.cut) # 1111307
-xy <- paste(x.cut,y.cut,sep="_")
-index30 <- match(xy,z.cut$xy)
-length(index30)
-any(is.na(index30))
-head(log10(z.cut$value[index30]))
-values.present <- vector(mode="numeric", length=length(xy))
-min(log10(z.cut$value[index30]), na.rm=T) #0
-max(log10(z.cut$value[index30]), na.rm=T) # 3.399154
-values.cut <- cut(log10(z.cut$value[index30]), breaks=seq(min(log10(z.cut$value[index30]), na.rm=T),
- 4, by=0.4))
-
-ordiplot(FD.rda4, type="n", xlim=c(-10,10),ylim=c(-10,12))
-points(plotcoord*15, cex=0.05, col=rev(greyPalette(n=10, start = 0, end = 0.8,
- gamma = 2.2, alpha = NULL))[values.cut])
-arrows(rep(0,1),rep(0,1),bioclimcoord[1]*8,bioclimcoord[2]*8, col="blue", length=0.1)
-text(bioclimcoord[1]*8,bioclimcoord[2]*8+0.5,"BIO_14", col="blue", cex=1)
-text(traitcoord,traitnames, col="red")
-# copied as third RDA graph into ppt
-
-FD.rda5 <- rda(FD ~ BIO_14+BIO_02,splot.bioclim[index14,c(4:49)],scale=T)
-FD.rda5
-'
-Inertia Proportion Rank
-Total 18.00000 1.00000
-Constrained 0.44582 0.02477 2
-Unconstrained 17.55418 0.97523 18
-Inertia is correlations
-
-Eigenvalues for constrained axes:
-RDA1 RDA2
-0.3641 0.0817
-
-Eigenvalues for unconstrained axes:
-PC1 PC2 PC3 PC4 PC5 PC6 PC7 PC8
-8.267 1.816 1.377 0.880 0.788 0.780 0.568 0.530 '
-traitnames <- dimnames(summary(FD.rda5)$species)[[1]]
-traitnames <- substr(traitnames,1,nchar(traitnames)-5)
-traitcoord <- summary(FD.rda5)$species[,c(1,2)]
-bioclimcoord <- summary(FD.rda5)$biplot[,c(1,2)]
-
-plotcoord <- summary(FD.rda5)$sites[,c(1,2)]
-str(plotcoord) #num [1:1111307, 1:2]
-x.cut <- cut(plotcoord[,1], breaks=seq(min(plotcoord[,1]), max(plotcoord[,1]), by=0.01))
-y.cut <- cut(plotcoord[,2], breaks=seq(min(plotcoord[,2]), max(plotcoord[,2]), by=0.01))
-z.cut <- table(x.cut,y.cut)
-head(z.cut)
-dim(z.cut) #131 50
-z.cut <- melt(z.cut)
-str(z.cut) # 6550 obs. of 3 variables:
-head(z.cut)
-z.cut$xy <- paste(z.cut$x.cut,z.cut$y.cut,sep="_")
-str(x.cut)
-length(x.cut) # 1111307
-xy <- paste(x.cut,y.cut,sep="_")
-index30 <- match(xy,z.cut$xy)
-length(index30)
-any(is.na(index30))
-head(log10(z.cut$value[index30]))
-values.present <- vector(mode="numeric", length=length(xy))
-min(log10(z.cut$value[index30]), na.rm=T) #0
-max(log10(z.cut$value[index30]), na.rm=T) # 3.399154
-values.cut <- cut(log10(z.cut$value[index30]), breaks=seq(min(log10(z.cut$value[index30]), na.rm=T),
- 4, by=0.4))
-
-ordiplot(FD.rda5, type="n", xlim=c(-10,10),ylim=c(-10,12))
-points(plotcoord*15, cex=0.05, col=rev(greyPalette(n=10, start = 0, end = 0.8,
- gamma = 2.2, alpha = NULL))[values.cut])
-arrows(rep(0,2),rep(0,2),bioclimcoord[,1]*4,bioclimcoord[,2]*4, col="blue", length=0.1)
-text(bioclimcoord*4,dimnames(bioclimcoord)[[1]], col="blue", cex=1)
-text(traitcoord,traitnames, col="red")
-# copied as fourth RDA graph into ppt
-
-plot(FD[,"Seed.length.mean"]~splot.bioclim$BIO_14[index14],
- cex=0.05, xlab="BIO_14", ylab="FD ln(seed length)",cex.lab=1.5, cex.axis=1.5)
-model1 <- lm(FD[,"Seed.length.mean"]~splot.bioclim$BIO_14[index14])
-summary(model1)
-' Estimate Std. Error t value Pr(>|t|)
-(Intercept) 2.584e-01 3.193e-04 809.1 <2e-16 ***
-splot.bioclim$BIO_14[index14] 1.407e-04 6.813e-07 206.5 <2e-16 ***
-Multiple R-squared: 0.03697, Adjusted R-squared: 0.03697
-'
-abline(model1, lwd=3, col="blue")
-#
-
-plot(FD[,"Seed.num.rep.unit.mean"]~splot.bioclim$BIO_02[index14],
- cex=0.05, xlab="BIO_02", ylab="FD ln(Seed number of reproductive unit)",cex.lab=1.5, cex.axis=1.5)
-model1 <- lm(FD[,"Seed.num.rep.unit.mean"]~splot.bioclim$BIO_02[index14])
-summary(model1)
-' Estimate Std. Error t value Pr(>|t|)
-(Intercept) 1.251e+00 1.799e-03 695.3 <2e-16 ***
-splot.bioclim$BIO_02[index14] -3.751e-03 1.927e-05 -194.6 <2e-16 ***
-Multiple R-squared: 0.03296, Adjusted R-squared: 0.03296 '
-abline(model1, lwd=3, col="blue")
-
-
-
-
-
-
-
-
-
-
-
-
-# mpd presence-absence
-system.time(
- mpd.sla.no.abu <- DT[,list(mpd.sla = mntd.fun.dt(Cover.., SLA.dist.mat, StandSpeciesName, abundance.weighted=FALSE)),by=PlotObservationID]
-)
-
-
-###########
-levels(splot.header2$COUNTRY_CO)
-levels(splot.header$Country)
-CWM2$Country <- splot.header$Country[index12]
-data(countryExData)
-str(countryExData)
-sPDF <- joinCountryData2Map(splot.header, joinCode = "NAME",nameJoinColumn = "Country")
-?mapCountryData
-#mapCountryData(sPDF, nameColumnToPlot="BIODIVERSITY" )
-
-#### Climatic space, from Jonathon Lenoir ####
-library(fBasics)
-# Color palette
-
-# Import bio1 (MAT) of the world
-bio <- read.table(paste("/home/helge/sPlot/Bioclim/", "BIO.txt", sep = ""),sep=";", header=T)
-str(bio) #9033453 obs. of 22 variables
-r <- raster(nrows=33, ncols=33, xmn=min(bio$BIO_01), xmx=max(bio$BIO_01), ymn=min(bio$BIO_12), ymx=max(bio$BIO_12))
-mattar <- data.frame(bio$BIO_01,bio$BIO_12)
-str(mattar)
-names(mattar) <- c("bio1","bio12")
-mattar_r <- rasterize(mattar, r, fun="count")
-str(mattar_r)
-plot(log10(mattar_r), asp=0, col=rev(divPalette(n=100, name=c("RdYlBu"))),
- main="N 2.5 arc-minutes cells per bin (log10 scale)", cex.main=0.7,
- xlab="MAT (°C*10)", ylab="TAR (mm)")
-# Bio1_Bio12_world.png
-
-#plot(mattar_r/max(getValues(mattar_r), na.rm=TRUE), asp=0, col=rev(divPalette(n=100, name=c("RdYlBu"))), main="Relative proportion of 10 arc-minutes cells per bin", xlab="MAR (°C*10)", ylab="TAR (mm)")
-# Plot the sampling effort of sPlot across the MAT-TAR climatic space
-str(splot.bioclim)
-# Make it a spatial data frame
-CRSlonlat <- CRS("+proj=longlat +datum=WGS84")
-splot.bioclim.coords <- cbind(splot.bioclim$X, splot.bioclim$Y)
-splot.bioclim.coords <- SpatialPoints(splot.bioclim.coords, proj4string=CRSlonlat)
-splot.bioclim2 <- data.frame(splot.bioclim$BIO_01,splot.bioclim$BIO_12)
-names(splot.bioclim2) <- c("bio1","bio12")
-#splot.bioclim2 <- SpatialPointsDataFrame(splot.bioclim.coords, splot.bioclim2, proj4string=CRSlonlat)
-#str(splot.bioclim2)
-#sPlot_mattar <- extract(bio1, coordinates(sPlot))
-#sPlot_mattar <- as.data.frame(sPlot_mattar)
-#names(sPlot_mattar) <- c("bio1")
-#sPlot_mattar$bio12 <- extract(bio12, coordinates(sPlot))
-sPlot_mattar_r <- rasterize(splot.bioclim2, r, fun="count")
-str(sPlot_mattar_r)
-plot(mattar_r, asp=0, col="grey", legend=FALSE, main="N vegetation plots per bin (log10 scale)",
- xlab="MAT (°C)", ylab="TAR (mm)")
-# Bio1_Bio12_world_grey.png
-plot(log10(sPlot_mattar_r), asp=0, col=rev(divPalette(n=100, name=c("RdYlBu"))),
- add=TRUE)
-# Bio1_Bio12_sPlot_SLA.png
-
-# Build the PC1-PC2 climatic space
-load("/home/helge/sPlot/Bioclim/pca_data.RData")
-str(pca.data)
-r <- raster(nrows=33, ncols=33, xmn=min(pca.data$PC1), xmx=max(pca.data$PC1), ymn=min(pca.data$PC2), ymx=max(pca.data$PC2))
-pca_r <- rasterize(pca.data[, c("PC1", "PC2")], r, fun="count")
-plot(log10(pca_r), asp=0, col=rev(divPalette(n=100, name=c("RdYlBu"))),
- main="N 2.5 arc-minutes cells per bin (log10 scale)",
- xlab="PC1 (cold-to-warm)", ylab="PC2 (dry-to-wet)")
-plot(pca_r/max(getValues(pca_r), na.rm=TRUE), asp=0,
- col=rev(divPalette(n=100, name=c("RdYlBu"))),
- main="Relative proportion of 2.5 arc-minutes cells per bin",
- xlab="PC1 (cold-to-warm)", ylab="PC2 (dry-to-wet)")
-# Plot the sampling effort of sPlot across the PC1-PC2 climatic space
-index13 <- match(splot.bioclim$RAST_ID,pca.data$RAST_ID)
-sPlot_pca <- pca.data[index13, c("PC1", "PC2")]
-sPlot_pca_r <- rasterize(sPlot_pca, r, fun="count")
-plot(pca_r, asp=0, col="grey", legend=FALSE,
- main="N vegetation plots per bin (log10 scale)",
- xlab="PC1 (cold-to-warm)", ylab="PC2 (dry-to-wet)")
-plot(log10(sPlot_pca_r), asp=0, col=rev(divPalette(n=100, name=c("RdYlBu"))),
- add=TRUE)
-
-#### For Marten Winter ####
-# species x country Matrix
-species.country <- table(DT$species,splot.header$Country[index30])
-str(species.country) # int [1:60908, 1:130]
-names(dimnames(species.country))
-names(dimnames(species.country)) <- c("species","country")
-
-library(reshape2)
-?melt.table
-species.country2 <- melt(species.country)
-str(species.country2) # 188470 obs. of 3 variables:
-names(species.country2)
-names(species.country2)[3] <- "p_a"
-table(species.country2$p_a)
-
-species.country2$p_a[species.country2$p_a>0] <- 1
-species.country2 <- species.country2[species.country2$p_a==1,]
-str(species.country2) # 188470 obs. of 3 variables:
-species.country2 <- species.country2[,-3]
-str(species.country2) # 188470
-write.csv(species.country2,file =
- paste("/home/oliver/shared/", "species.country2.csv", sep = ""), row.names = FALSE)
-levels(species.country2$country) # 130
-length(species.country2$country[species.country2$country=="Algeria"])
-# 97
-sum(species.country[,2])
-# 799, this is the amount of records for DZA = Algeria
-sum(species.country[,126])
-# 1467200, sum of records for United States
-
-TaxaNat <- read.csv(paste("/home/oliver/shared/Glonaf/", "TaxaNat.csv", sep = ""))
-str(TaxaNat) #178779 obs. of 2
-
-GlonafRegion <- read.csv(paste("/home/oliver/shared/Glonaf/", "GlonafRegion.csv", sep = ""))
-str(GlonafRegion) #843 obs. of 10
-trim.trailing <- function (x) sub("\\s+$", "", x)
-# returns string w/o leading or trailing whitespace
-GlonafRegion$Country_sPlot <- as.factor(trim.trailing(as.character(GlonafRegion$Country_sPlot)))
-index14 <- match(levels(species.country2$country),levels(GlonafRegion$Country_sPlot))
-index14
-levels(GlonafRegion$Country)
-
-data.frame(levels(species.country2$country),levels(GlonafRegion$Country_sPlot)[index14])
-## all correct
-
-index15 <- match(levels(GlonafRegion$Country_sPlot),GlonafRegion$Country_sPlot)
-index15
-index16 <- match(species.country2$country,levels(species.country2$country))
-any(is.na(index16)) # F
-length(levels(species.country2$species)) #60908
-'species.country2$species.country <- paste(species.country2$species,
-levels(species.country2$country)[index16],sep="_")
-species.country2$species.country <- paste(species.country2$species,
-levels(GlonafRegion$Country_sPlot)[index14][index16],sep="_")
-species.country2$species.country <- paste(species.country2$species,
-GlonafRegion$Country_sPlot[index15][index14][index16],sep="_")
-'
-species.country2$species.country <- paste(species.country2$species,
- GlonafRegion$ISO.country[index15][index14][index16],sep="_")
-str(species.country2) # 188470 obs. of 3 variables:
-head(species.country2)
-species.country2[c(160000:160100),]
-## fits!
-
-index17 <- match(TaxaNat$Region,GlonafRegion$ID.region)
-any(is.na(index17)) # F
-
-TaxaNat$species.country <-paste(TaxaNat$Taxon2,GlonafRegion$ISO.country[index17], sep="_")
-
-index18 <- match(species.country2$species.country,TaxaNat$species.country)
-length(index18) #188470
-length(index18[!is.na(index18)]) # 9283 out of 188470 species x country occurrences
-# are alien occurrences
-9283/188470*100 # 4.925452 % species x country occurrences
-species.country2$alien <- F
-species.country2$alien[!is.na(index18)] <- T
-species.country2$ISO.country <- GlonafRegion$ISO.country[index15][index14][index16]
-str(species.country2)
-head(species.country2)
-length(species.country2$ISO.country[species.country2$ISO.country=="DZA"])
-# 12939, this cannot be
-#species.country2$ISO.country[species.country2$ISO.country=="DZA"]
-sum(species.country[,2])
-# 799, this is the amount of records for DZA = Algeria
-length(species.country2$ISO.country[species.country2$ISO.country=="DZA"
- & !is.na(species.country2$ISO.country)])
-# 97!
-
-#table(species.country2$ISO.country)
-#species.country3 <- acast(species.country2, species~ISO.country, length)
-length(species.country2$alien[species.country2$alien])
-# 9283 alien records in sPlot
-length(species.country2$alien[species.country2$alien & species.country2$ISO.country=="DZA"])
-# 3 alien records in DZA = Algeria
-length(species.country[species.country[,2]>0,2])
-# 97, however there are only 97 records in total
-length(species.country2$alien[species.country2$alien & species.country2$ISO.country=="USA"])
-# 1911 alien species in USA
-length(species.country[species.country[,126]>0,2])
-# 9933 species in total in USA
-alien.richness <- as.matrix(table(species.country2$ISO.country,species.country2$alien))[,2]
-
-species.country3 <- species.country
-species.country3[species.country3>0] <- 1
-country.richness <- colSums(species.country3)
-country.richness
-names(country.richness)
-data.frame(names(country.richness),GlonafRegion$ISO.country[index15][index14])
-country.richness2 <- data.frame(GlonafRegion$ISO.country[index15][index14],as.numeric(country.richness))
-names(country.richness2) <- c("ISO.country","richness")
-country.richness2$ISO.country
-country.richness2
-
-index19 <- match(country.richness2$ISO.country,names(alien.richness))
-index19
-alien.richness
-country.richness2$alien.richness <- alien.richness[index19]
-country.richness2
-country.richness2$alien.prop <- country.richness2$alien.richness/country.richness2$richness
-country.richness2$alien.prop
-
-sPDF <- joinCountryData2Map(country.richness2,joinCode = "ISO3",
- nameJoinColumn = "ISO.country")
-'112 codes from your data successfully matched countries in the map
-18 codes from your data failed to match with a country code in the map
-132 codes from the map werent represented in your data
-'
-mapParams <- mapCountryData(sPDF, nameColumnToPlot="richness",addLegend=FALSE)
-do.call(addMapLegend, c(mapParams, legendWidth=0.5, legendMar = 2,legendLabels="all"))
-mapParams <- mapCountryData(sPDF, nameColumnToPlot="richness",
- catMethod=seq(0,10000,by=1000),addLegend=FALSE)
-do.call(addMapLegend, c(mapParams, legendWidth=0.5, legendMar = 2,legendLabels="all"))
-
-mapParams <- mapCountryData(sPDF, nameColumnToPlot="alien.richness",
- catMethod=seq(0,2000,by=200),addLegend=FALSE)
-do.call(addMapLegend, c(mapParams, legendWidth=0.5, legendMar = 2,legendLabels="all"))
-
-max(country.richness2$alien.prop, na.rm=T) #0.2823779
-mapParams <- mapCountryData(sPDF, nameColumnToPlot="alien.prop",addLegend=FALSE)
-mapParams <- mapCountryData(sPDF, nameColumnToPlot="alien.prop",
- catMethod=seq(0,0.3,by=0.02),addLegend=FALSE)
-do.call(addMapLegend, c(mapParams, legendWidth=0.5, legendMar = 2,legendLabels="all"))
-
-
-
-
-##################### OLD ###################
-library(spdep)
-library(maptools)
-library(hexbin)
-library(RColorBrewer)
-# Use a lighter version of sPlot focusing on mean annual temperature MAT (BIO1)
-# and total annual rainfall TAR (BIO12)
-sPlot <- splot.header[, c("LONGITUDE", "LATITUDE", "BIO_1", "BIO_12")]
-sPlot$LONGITUDE <- as.numeric(as.character(sPlot$LONGITUDE))
-sPlot$LATITUDE <- as.numeric(as.character(sPlot$LATITUDE))
-sPlot$BIO_1 <- as.numeric(sPlot$BIO_1)/10
-sPlot$BIO_12 <- as.numeric(sPlot$BIO_12)
-str(sPlot)
-# Make it a spatial data frame
-CRSlonlat <- CRS("+proj=longlat +datum=WGS84")
-coords <- cbind(sPlot$LONGITUDE, sPlot$LATITUDE)
-coords <- SpatialPoints(coords, proj4string=CRSlonlat)
-sPlot <- SpatialPointsDataFrame(coords, sPlot, proj4string=CRSlonlat)
-str(sPlot)
-# Plot the data on the world map to check the coordinates and projection
-data(wrld_simpl)
-plot(wrld_simpl)
-points(coordinates(sPlot), cex=0.1, col="red", pch=16)
-# Set a color ramp
-rf <- colorRampPalette(rev(brewer.pal(11, 'Spectral')))
-# Plot the plot distribution within the MAT-TAR bidimentionnal climatic space
-hexbinplot(BIO_12~BIO_1, data=sPlot@data, colramp=rf, trans=log, inv=exp,
- border="white", xlab="MAT", ylab="TAR", xbins=30)
-# Display only bins with a minimum of 10 plots and a maximum of 1000 plots
-hexbinplot(BIO_12~BIO_1, data=sPlot@data, colramp=rf, trans=log, inv=exp,
- border="white", xlab="MAT", ylab="TAR", xbins=30, mincnt=10, maxcnt=1000)
-#
-
-
-
-
-
-
-
-
-################ OLD ###################
-
-###
-load("C:\\Daten\\iDiv\\sPlot\\pca_data.Rdata")
-str(pca.data)
-any(is.na(pca.data$PC1)) #F
-
-load("C:\\Daten\\iDiv\\sPlot\\sPlot_header_pca.RData")
-str(splot.header2)
-any(is.na(pca.data$PC1)) #F
-
-library(raster)
-r <-
- splot_pca_r <- rasterize(splot_pca, r, run="count")
-plot(splot_pca_r, asp=0, col=rev(divPalette(n=100, name="RdY1Bu"))))
-plot(splot_r, asp=0, col=rev(divPalette(n=100, name="RdY1Bu"))))
-
-index <- which(getValues(splot_pca_r)>0)
-which(getValues(pca_r)>0)
-poinra
-getValues(splot_pca_r)[index] # 136 values
-?sample
-splo_pca$binid <- cellFromXY(pca_r, splot_pca)
-
-
-###########################
-load("C:\\Daten\\iDiv\\sPlot\\splot_pca2.Rdata")
-str(splot_pca)
-splot_pca <- data.table(splot_pca)
-setkey(splot_pca,binid)
-inter <- splot_pca[,list(count=length(PCA1)), by=binid]
-splot_pca2 <- splot_pca[,list(PCA1, PCA2,count=length(PCA1)), by=binid]
-
-str(splot_pca2)
-hist(splot_pca2$binid)
-inter2 <- splot_pca2[count>100,list(count=sum(PCA1)), by=binid]
-str(inter2)
-[,list(total.cover=sum(Cover..)), by=PlotObservationID]
-<- DT[,list(total.cover=sum(Cover..)), by=PlotObservationID]
-
-sample(splot_pca$PCA1,100)
-splot_pca2 <- splot_pca[table(splot_pca$binid)>100,sample(PCA1,100),by=binid]
-str(splot_pca2)
-
-resample <- function(PCA1_,PCA2_,binid_){
- if (length(PCA1)>100) {
-
- }
- else {
- 1
- }
- ]index10 <- match(species,labels(d))
- d3 <- as.dist(dist[index10,index10])
- if (sum(d3)==0){
- NA
- }
- divc(as.data.frame(rel.cover),d3)$diversity
-}
-splot_pca2 <- splot_pca[length(PCA1)<=100,res_ID := 1]
-head(splot_pca2,100)
-table(splot_pca2$binid,splot_pca2$count, exclude=NULL)
-table(splot_pca$binid)
-splot_pca2 <- splot_pca[,res_ID := 0]
-splot_pca2 <- splot_pca[table(splot_pca2$binid)<=100,res_ID := 1]
-splot.selection <- splot_pca[,list(res = resample(PCA1,PCA2,binid)),by=binid]
-#length(splot_pca[binid>100])
-hist(splot_pca$binid)
-splot.selection <- splot_pca[,list(length(binid)),by=binid]
-splot.selection <- splot_pca[,list(count=length(PCA1)),by=binid]
-splot.selection <- splot_pca[,,by=binid]
-hist(splot.selection$count)
-str(splot.selection)
-str(splot_pca2)
-splot.selection <- splot_pca[,list(nrow()),by=binid]
-splot.selection <- splot_pca[,count:=length(binid),by=binid]
-length(splot_pca,by=binid)
-########################################################################################
-
-### CWM ###
-#CWM1 <- DT[!is.na(plot.species.trait.value),list(CWM = weighted.mean(plot.species.trait.value,Relative.cover)),by=PlotObservationID]
-CWM1 <- DT[,list(CWM = weighted.mean(plot.species.trait.value,Relative.cover)),by=PlotObservationID]
-# see http://stackoverflow.com/questions/3367190/aggregate-and-weighted-mean-in-r
-CWM1[1:100,]
-hist(CWM1$CWM)
-write.csv((x, file = paste(path.sPlot, "CWM_SLA_all.csv", sep = ""), row.names = FALSE))
-
-### FD ###
-library(ade4) # f?r divc (Rao's Q)
-# this does not work: ~50 GB RAM required for 40k x 40k species matrix:
-'ktab1 <- ktab.list.df(list(as.data.frame(species.trait.value[!is.na(species.trait.value)])))
-d <- dist.ktab(ktab1, type=c("Q"), option = "scaledBYrange")
-#labels(d)
-#str(d)
-#attr(d,"Labels")
-attr(d,"Labels") <- dimnames(species.trait.value)[[1]][!is.na(species.trait.value)]
-labels(d)
-any(is.na(d)) # F
-dist <- as.matrix(d)
-
-#tree.trait.weighted <- t(t(veg.KHG3.prop) * all.traits.KHG[,j])
-#dist.sqrt <- sqrt(dist)
-d.sqrt <- sqrt(d)
-labels(d.sqrt)
-dist <- as.matrix(d.sqrt)
-'
-'Rao <- function(rel.cover,species){
-index10 <- match(species,labels(d))
-d3 <- as.dist(dist[index10,index10])
-if (sum(d3)==0){
-NA
-}
-divc(as.data.frame(rel.cover),d3)$diversity
-}
-FD.all <- DT[!is.na(plot.species.trait.value),list(Rao.Q = Rao(Relative.cover,Matched.concept)),by=PlotObservationID]
-'
-#species.trait.value.scaled <- species.trait.value/(max(species.trait.value)-min(species.trait.value))
-species.trait.value.scaled <- TRY.mean.sd.splot$SLA.mean/(max(TRY.mean.sd.splot$SLA.mean)-min(TRY.mean.sd.splot$SLA.mean))
-
-# scale by range, which saves time for distance calculations
-
-library(reshape2)
-
-package.size = 1000
-#number.plots.total <- length(levels(DT[!is.na(plot.species.trait.value)]$PlotObservationID)) # 705272
-number.plots.total <- length(levels(DT$PlotObservationID)) # 701611 vs. 705272
-FD.all <- array(0,number.plots.total, dimnames=list(levels(DT$PlotObservationID)))
-number.packages <-ceiling(number.plots.total/package.size) #702
-pb <- winProgressBar(title = "progress bar", min = 0, max =number.packages, width = 300)
-for (k in 1:number.packages) {
- setWinProgressBar(pb, k, title=paste("Row ",k,"of ", number.packages, " is processed"))
- ifelse (k==number.packages, package<-number.plots.total%%package.size, package <- package.size)
- #index3 <- match(DT[!is.na(plot.species.trait.value)]$PlotObservationID,levels(DT[!is.na(plot.species.trait.value)]$PlotObservationID)[seq((k-1)*package.size+1,k*package.size,1)])
- index3 <- match(DT$PlotObservationID,levels(DT$PlotObservationID)[seq((k-1)*package.size+1,k*package.size,1)])
- #DT2 <- DT[!is.na(plot.species.trait.value) & !is.na(index3)]
- #DT2 <- DT2[!is.na(plot.species.trait.value)]
- #length(index3[!is.na(index3)])
- DT2 <- DT[!is.na(index3)]
- #str(DT2)
- #veg.matrix1 <- acast(DT2, Matched.concept ~ PlotObservationID, value.var="Relative.cover", fill=0)
- veg.matrix1 <- acast(DT2, names.stand ~ PlotObservationID, value.var="Relative.cover", fun.aggregate=sum, fill=0)
- # fun.aggregate = sum -> sums up multiple entries per species!
- # this can happen with the unified species nomenclature!
- #length(unique(dimnames(veg.matrix1)[[2]]))# 1000
- #length(unique(DT2$PlotObservationID)) # 1000
- #veg.matrix1[veg.matrix1 > 0] <- 1
- #any(colSums(veg.matrix1)!=1) #F
- index4 <- match(dimnames(veg.matrix1)[[1]],species.list)
- species.trait.value.scaled2 <- species.trait.value.scaled[index4]
- species.list2 <- species.list[index4]
- ktab1 <- ktab.list.df(list(as.data.frame(species.trait.value.scaled2)))
- d <- dist.ktab(ktab1, type=c("Q"), option = "noscale")
- #labels(d)
- #str(d)
- #attr(d,"Labels")
- attr(d,"Labels") <- species.list2
- #labels(d)
- #any(is.na(d)) # F
- #dist <- as.matrix(d)
- #as.matrix(d)[1:8,1:8]
- d <- sqrt(d)
- #dist <- as.matrix(d)
- #rm(d)
- if (k<number.packages){
- FD.all[seq((k-1)*package.size+1,k*package.size,1)] <-
- divc(as.data.frame(veg.matrix1),d)$diversity
- #DT2[,list(Rao.Q = Rao(Relative.cover,Matched.concept)),by=PlotObservationID]$Rao.Q
- any(colSums(veg.matrix1)==0) # F o.K.
- any(colSums(veg.matrix1)!=1) # T ups....
- hist(colSums(veg.matrix1))
- any(rowSums(veg.matrix1)==0) # F o.K.
- } else {
- FD.all[seq((k-1)*package.size+1,number.plots.total,1)] <-
- divc(as.data.frame(veg.matrix1),d)$diversity
- #DT2[,list(Rao.Q = Rao(Relative.cover,Matched.concept)),by=PlotObservationID]
- }
-}
-close(pb)
-#FD.all
-write.csv(FD.all, file = paste(path.sPlot, "FDall.csv", sep = ""), row.names = FALSE)
-