# Script file for Elevations # R version 2.1.1 # see also OldWolfeErrorAnalysisSriptFile.txt and `Tightening the CLAMP' #################### #Preprocessing data #Read in the raw files Wolfe173ALL <- read.table("/Users/wag/Archive/Active/Ecophenetics and Comp. Index/CLAMP Database/Wolfe173ALL.txt", sep = '\t', header = TRUE, colClasses = 'character') JacobsALL <- read.table("/Users/wag/Archive/Active/Ecophenetics and Comp. Index/CLAMP Database/JacobsALL.txt", sep = '\t', header = TRUE, colClasses = 'character') GregoryALL <- read.table("/Users/wag/Archive/Active/Ecophenetics and Comp. Index/CLAMP Database/GregoryALL.txt", sep = '\t', header = TRUE, colClasses = 'character') KowalskiALL <- read.table("/Users/wag/Archive/Active/Ecophenetics and Comp. Index/CLAMP Database/KowalskiALL.txt", sep = '\t', header = TRUE, colClasses = 'character') # Abbreviations for the CLAMP variables CLAMP.vars.abbrev31 <- c('Lobd', 'Entr', 'TReg', 'TCls', 'TRnd', 'TAcu', 'TCmp', 'ZNan', 'ZLe1', 'ZLe2', 'ZMi1', 'ZMi2', 'ZMi3', 'ZMe1', 'ZMe2', 'ZMe3', 'AEmg', 'ARnd', 'AAcu', 'AAtn', 'BCor', 'BRnd', 'BAcu', 'Rlt1', 'Rb12', 'Rb23', 'Rb34', 'Rgt4', 'SObo', 'SElp', 'SOvt') # Make the CLAMP codes numeric in the Wolfe data set Wolfe173CLAMP <- Wolfe173ALL[,1:31] for(i in 1:31){ Wolfe173CLAMP[,i] <- as.numeric(Wolfe173CLAMP[,i]) } # Reorder and make the CLAMP codes numeric in the Jacobs data set JacobsCLAMP29 <- JacobsALL[,c(3, 2, 4, 8, 6, 7, 5, 24:30, 47, 10, 12, 11, 15, 13, 14, 19:23, 17, 18, 16)] for(i in 1:29){ JacobsCLAMP29[,i] <- 100 * as.numeric(JacobsCLAMP29[,i]) } JacobsCLAMP <- cbind(JacobsCLAMP29[,c(1:7)], rep(NA, 30), JacobsCLAMP29[,c(8:14)], rep(NA, 30), JacobsCLAMP29[,c(15:29)]) # Make the CLAMP codes numeric in the Gregory data set GregoryCLAMP <- GregoryALL[,20:50] for(i in 1:31){ GregoryCLAMP[,i] <- as.numeric(GregoryCLAMP[,i]) } # Make the CLAMP codes numeric in the Kowalski data set KowalskiCLAMP <- KowalskiALL[,1:31] for(i in 1:31){ KowalskiCLAMP[,i] <- 100 * as.numeric(KowalskiCLAMP[,i]) } # Make all the names the same abbreviations names (Wolfe173CLAMP) <- CLAMP.vars.abbrev31 names (JacobsCLAMP29) <- CLAMP.vars.abbrev31[c(1:7, 9:15, 17:31)] names (JacobsCLAMP) <- CLAMP.vars.abbrev31 names (GregoryCLAMP) <- CLAMP.vars.abbrev31 names (KowalskiCLAMP) <- CLAMP.vars.abbrev31 # Add row names to the Jacobs data row.names (JacobsCLAMP29) <- c('Mpanga Forest, UG', 'Combretaceous Woodland, UG', 'Arabuko-Sokoke Natural Reserve, KE', 'Mbololo Teita, KE', 'Mramagambo, UG', 'Bakundu, CAM', 'Tree Savanna, UG', 'Mt. Kenya, Meru, KE', 'Magamba Forest, TZ', 'Nidrowa, Okomu, NI', 'Budongo, UG', 'Lupa, TZ', 'Kibwezi, KE', 'Albizia Woodland, UG', 'Lengwe, MAL', 'Mwaniliana, Lowland, TZ', 'Usambara, Submontane, TZ', 'Kora, KE', 'St. Florais, CAR', 'Buganda Thicket, UG', ' Yangambi, CO', 'Madabu, CO', 'Kade, CO', 'Lobe, GAB', 'Kahuzi, CO', 'Mole, GH', 'Kaniama, CO', 'Wamba, CO', 'Manyara, TZ', 'Koitoboss, Elgon, KE') row.names (JacobsCLAMP) <- row.names (JacobsCLAMP29) # Put together the CLAMP and MAT data ALLCLAMP <- rbind (Wolfe173CLAMP, JacobsCLAMP, GregoryCLAMP, KowalskiCLAMP) ALLMAT <- as.numeric(c(Wolfe173ALL[,33], JacobsALL[,32], GregoryALL[,13], KowalskiALL[,36])) # Read in the supplementary data from Wolfe 1993 Wolfe1993 <- read.table("/Users/wag/Archive/Active/Ecophenetics and Comp. Index/CLAMP Database/Wolfe1993Data.txt", sep = '\t', header = TRUE, colClasses = 'character') # Create factor to indicate which study each flora came from Study <- as.factor(c(rep (1, 173), rep (2, 30), rep (3, 12), rep (4,30))) levels(Study) <- c('Wolfe173', 'Jacobs', 'Gregory', 'Kowalski') # Numbers of coded OTUs ALLn.coded <- as.numeric(c(Wolfe1993$n, JacobsALL[,43], GregoryALL[,1], KowalskiALL[,41])) # Determine total numbers (of OTUs) ALLn.total <- as.numeric(c(rep(NA, 173), JacobsALL[,42], rep(NA, 12), KowalskiALL[,40])) # Deal with Lat., Long., and elevation # N latitude and E longitude are positive # Wolfe LatDD <- as.numeric(Wolfe1993[,1]) + as.numeric(Wolfe1993[,2]) / 60 LongDD <- as.numeric(Wolfe1993[,4]) + as.numeric(Wolfe1993[,5]) / 60 WolfeSpatial <- data.frame(LatDeg = I(as.numeric(Wolfe1993[,1])), LatMin = I(as.numeric(Wolfe1993[,2])), LatSec = rep(NA, nrow(Wolfe1993)), LatDD = LatDD, NS = rep ('N', nrow(Wolfe1993)), LongDeg = I(as.numeric(Wolfe1993[,4])), LongMin = I(as.numeric(Wolfe1993[,5])), LongSec = rep(NA, nrow(Wolfe1993)), LongDD = LongDD, EW = rep ('W', nrow(Wolfe1993)), Elev = I(as.numeric(Wolfe1993$elev))) WolfeSpatial$NS[is.na(WolfeSpatial$LatDeg)] <- NA WolfeSpatial$EW[is.na(WolfeSpatial$LongDeg)] <- NA levels(WolfeSpatial$NS) <- list(N='N', S='S', E='E', W='W') levels(WolfeSpatial$EW) <- list(N='N', S='S', E='E', W='W') row.names(WolfeSpatial) <- row.names (Wolfe173ALL) #All the Wolfe data is N lat and W long, so just reverse the sign of the # decimal degree long....this may not work in general. WolfeSpatial$LongDD <- 0 - WolfeSpatial$LongDD # Fix a little typo in Jacobs JacobsALL$Longitude[4] <- '38deg15minE' # Jacobs JacobsLat <- t(as.data.frame(strsplit(as.character(JacobsALL$Latitude), split = '[a-z]{1,3}'))) JacobsLong <- t(as.data.frame(strsplit(JacobsALL$Longitude, split = '[a-z]{1,3}'))) LatDD <- as.numeric(JacobsLat[,1]) + as.numeric(JacobsLat[,2]) / 60 LongDD <- as.numeric(JacobsLong[,1]) + as.numeric(JacobsLong[,2]) / 60 JacobsSpatial <- data.frame(LatDeg = I(as.numeric(JacobsLat[,1])), LatMin = I(as.numeric(JacobsLat[,2])), LatSec = rep(NA, nrow(JacobsLat)), LatDD = LatDD, NS = JacobsLat[,3], LongDeg = I(as.numeric(JacobsLong[,1])), LongMin = I(as.numeric(JacobsLong[,2])), LongSec = rep(NA, nrow(JacobsLong)), LongDD = LongDD, EW = JacobsLong[,3], Elev = rep(NA, nrow(JacobsLong))) levels(JacobsSpatial$NS) <- list(N='N', S='S', E='E', W='W') levels(JacobsSpatial$EW) <- list(N='N', S='S', E='E', W='W') row.names(JacobsSpatial) <- row.names(JacobsALL) JacobsSpatial$LatDD[JacobsSpatial$NS == 'S'] <- 0 - JacobsSpatial$LatDD[JacobsSpatial$NS == 'S'] JacobsSpatial$LongDD[JacobsSpatial$EW == 'W'] <- 0 - JacobsSpatial$LongDD[JacobsSpatial$EW == 'W'] # Gregory GregorySpatial <- data.frame(LatDeg = rep(NA, nrow(GregoryALL)), LatMin = rep(NA, nrow(GregoryALL)), LatSec = rep(NA, nrow(GregoryALL)), LatDD = abs(as.numeric(GregoryALL$Lat.)), NS = rep('S', nrow(GregoryALL)), LongDeg = rep(NA, nrow(GregoryALL)), LongMin = rep(NA, nrow(GregoryALL)), LongSec = rep(NA, nrow(GregoryALL)), LongDD = abs(as.numeric(GregoryALL$Long.)), EW = rep('W', nrow(GregoryALL)), Elev = as.numeric(gsub('[^0-9]', '', GregoryALL$Z.m.))) levels(GregorySpatial$NS) <- list(N='N', S='S', E='E', W='W') levels(GregorySpatial$EW) <- list(N='N', S='S', E='E', W='W') row.names(GregorySpatial) <- row.names(GregoryALL) GregorySpatial$LatDD[GregorySpatial$NS == 'S'] <- 0 - GregorySpatial$LatDD[GregorySpatial$NS == 'S'] GregorySpatial$LongDD[GregorySpatial$EW == 'W'] <- 0 - GregorySpatial$LongDD[GregorySpatial$EW == 'W'] #Kowalski Pos <- t(as.data.frame(strsplit(KowalskiALL$Long., split = '[^0-9]'))) LatDD <- as.numeric(Pos[,1]) + as.numeric(Pos[,4]) / 60 LongDD <- as.numeric(Pos[,6]) + as.numeric(Pos[,9]) / 60 Elevs <- (as.numeric(KowalskiALL$Max..Elev...m.) + as.numeric(KowalskiALL$Min..Elev...m.)) / 2 KowalskiSpatial <- data.frame(LatDeg = as.numeric(Pos[,1]), LatMin = as.numeric(Pos[,2]), LatSec = rep(NA, nrow(KowalskiALL)), LatDD = LatDD, NS = substr(gsub('[^NSEW]', '', KowalskiALL$Long.), 1, 1), LongDeg = as.numeric(Pos[,4]), LongMin = as.numeric(Pos[,5]), LongSec = rep(NA, nrow(KowalskiALL)), LongDD = LongDD, EW = substr(gsub('[^NSEW]', '', KowalskiALL$Long.), 2, 2), Elev = Elevs) levels(KowalskiSpatial$NS) <- list(N='N', S='S', E='E', W='W') levels(KowalskiSpatial$EW) <- list(N='N', S='S', E='E', W='W') row.names(KowalskiSpatial) <- row.names(KowalskiALL) KowalskiSpatial$LatDD[KowalskiSpatial$NS == 'S'] <- 0 - KowalskiSpatial$LatDD[KowalskiSpatial$NS == 'S'] KowalskiSpatial$LongDD[KowalskiSpatial$EW == 'W'] <- 0 - KowalskiSpatial$LongDD[KowalskiSpatial$EW == 'W'] # Put the spatial data together ALLSpatial <- rbind (WolfeSpatial, JacobsSpatial, GregorySpatial, KowalskiSpatial) plot(ALLSpatial$LongDD, ALLSpatial$LatDD, col = as.numeric(Study)) #text(ALLSpatial$LongDD, ALLSpatial$LatDD, labels = 1:245) abline(h=0) abline(v=0) #################### # Some more climate variables ALLMAP <- as.numeric(c(rep(NA, 173), JacobsALL[,36], rep(NA, 12), KowalskiALL[,37])) ALLGSP <- as.numeric(c(Wolfe173ALL[,37], rep(NA, 30), GregoryALL[,15], rep(NA, 30))) GregoryALL[10,14] <- 17.4 ALLWMMT <- as.numeric(c(Wolfe173ALL[,34], JacobsALL[,33], GregoryALL[,14], rep(NA, 30))) ALLCMMT <- as.numeric(c(Wolfe173ALL[,35], JacobsALL[,34], rep(NA, 12), rep(NA, 30))) ALL <- cbind(n.coded = ALLn.coded, Study = Study, ALLSpatial, ALLCLAMP, MAT = ALLMAT, MAP = ALLMAP, GSP = ALLGSP, WMMT = ALLWMMT, CMMT = ALLCMMT) #################### # More sensible morphological variables #take the size-class bins and replace them with a decent continuous variable: lengths <- c(NA, 10, 20, 41, 77, 116, 160, NA, NA) widths <- c(NA, 3, 7, 15, 28, 45, 59, NA, NA) areas <- (2/3) * (lengths * widths) mean.area <- rep (0, 245) for (i in 1:245){ area.proportions <- c(ALLCLAMP$ZNan[i], ALLCLAMP$ZLe2[i], ALLCLAMP$ZLe2[i], ALLCLAMP$ZMi1[i], ALLCLAMP$ZMi2[i], ALLCLAMP$ZMi3[i], ALLCLAMP$ZMe1[i], ALLCLAMP$ZMe2[i], ALLCLAMP$ZMe3[i]) mean.area[i] <- mean(na.omit(areas * area.proportions)) } pairs(cbind(elev = ALLSpatial$Elev, mean.area, lat = ALLSpatial$LatDeg, entr = ALLCLAMP$Entr), col = as.numeric(Study)) summary(lm1 <- lm(ALLSpatial$Elev ~ mean.area)) summary(lm2 <- lm(ALLSpatial$Elev ~ mean.area + ALLCLAMP$Entr)) summary(lm3 <- lm(ALLSpatial$Elev ~ mean.area + ALLSpatial$LatDeg)) summary(lm4 <- lm(ALLSpatial$Elev ~ mean.area + ALLCLAMP$Entr + ALLSpatial$LatDeg)) summary(lm5 <- lm(ALLSpatial$Elev ~ mean.area * ALLCLAMP$Entr * ALLSpatial$LatDeg)) models <- list(lm1, lm2, lm3, lm4, lm5) nobs <- c(117, 117, 105, 105, 105) aics <- rep(0, 5) bics <- rep(0, 5) for(i in 1:5){ aics[i] <- AIC(models[[i]], k = 2) bics[i] <- AIC(models[[i]], k = log(nobs[i])) } data.frame(cbind(Model = 1:5, AIC = aics, BIC = bics)) #and the same for length-to-width-ratio: aspects <- 0.5, 1.5, 2.5, 3.5, 7 mean.aspect <- rep(0, 245) for (i in 1:245){ aspect.proportions <- c(ALLCLAMP$R<1:, ALLCLAMP$R1:2, ALLCLAMP$R2:3, ALLCLAMP$R3:4, ALLCLAMP$R:>4) mean.aspect[i] <- mean(na.omit(aspects * aspect.proportions)) } #names(KowalskiALL)[32:44] <- c('Lat.', 'Long.', 'Elev.min.m.', 'Elev.max.m.', 'MAT.C', 'MAPcm.', 'DryMos', 'MATDur', 'nspp', 'nscored', 'Area', 'TypeSampled', 'For.Type')