########################################################### #### Calculation of activity and restriction times #### ########################################################### # Script by Octavio Jiménez-Robles. octavio.jimenez.robles (at) gmail.com # Used in the analyses of Jiménez-Robles & De la Riva. 2019. # Lizards in the mist: Thermal niches constrained by habitat and microclimates in the Andes of southern Bolivia. # Journal of Biogeography. #--------------------------------------------------------------------------------------------------------------- # 1. Quantifying Activity and Restriction Times #### # Reading Operative Temperatures TeSam2012 <-read.table("input_data/TeSam2012.csv", sep= ";", header=TRUE) # Data from the campaign 2012-2013 TeSam2013 <-read.table("input_data/TeSam2013.csv", sep= ";", header=TRUE) # Data from the campaign 2013-2014 # Physiological limits species.fisio <- read.table("output_data/SpeciesPhysiology.csv", sep=";", header=T) # Created in "Lizards in the mist. Script 1.Summary of thermal thresholds of activity.R" # Adjustment of heating and cooling rates TeCorrections <- read.table("output_data/TeCorrections.csv", sep=";", header=T) # Created in "Lizards in the mist. Script 1.Summary of thermal thresholds of activity.R" species.list <- list("puna", "ornatus", "sp", "orientalis") Transectslist <- c("1A", "1B", "1C", "1D", "1E", "1F", "1G", "1H", "1I", "1J", "2A", "2B", "2C", "2D", "2E", "2F", "2G", "2H", "2I", "2J", "3A", "3B", "3C", "3D", "3E", "3F", "3G", "3H", "3I", "3J") source("functions/degreeMinutes.R") for(l in 1:length(species.list)){ Vmax <- species.fisio$TVmax[l] TEmerge <- species.fisio$TEmerge[l] VmaxP <- species.fisio$TVmaxP[l] TEmergeP <- species.fisio$TEmergeP[l] if(l==4){ TEmergeP <- TEmerge} # This is because we did not have a representative number of Tb data from active pregnant females of Liolaemus orientalis. So we decided to use the same lower thermal threshold than the rest of individuals of the species. intercept <- TeCorrections$MixedInt[l] slope <- TeCorrections$MixedSlp[l] Activity12 <- data.frame(Days = levels(as.factor(TeSam2012$JULDAY)), matrix(ncol=30, NA)) Activity13 <- data.frame(Days = levels(as.factor(TeSam2013$JULDAY)), matrix(ncol=30, NA)) names(Activity12)[-1] <-Transectslist names(Activity13) <-names(Activity12) Restriction12 <- Restriction12P <- Restriction12cali <- Restriction12caliP <- Activity12P <- Activity12cali <- Activity12caliP <- Activity12 Restriction13 <- Restriction13P <- Restriction13cali <- Restriction13caliP <- Activity13P <- Activity13cali <- Activity13caliP <- Activity13 for (i in 1:length(Transectslist)){ Transect <- Transectslist[[i]] TransectLETTER<-c("HA", "HB", "HC", "HD", "HE", "HF", "HG", "HH", "HI", "HJ", "IA", "IB", "IC", "ID", "IE", "IF", "IG", "IH", "II", "IJ", "LA", "LB", "LC", "LD", "LE", "LF", "LG", "LH", "LI", "LJ")[[i]] TeTransect12 <- TeSam2012[, c(1,2,grep(TransectLETTER, substring(colnames(TeSam2012),1,2)))] TeTransect13 <- TeSam2013[, c(1,2,grep(TransectLETTER, substring(colnames(TeSam2013),1,2)))] if(i>0&i<6 | i==9| i>10&i<16 | i>20&i<26){ for (j in 1:length(levels(as.factor(TeSam2012$JULDAY)))){ Fragment12 <- subset(TeTransect12, JULDAY==levels(as.factor(TeSam2012$JULDAY))[[j]]) Fragment12 <- Fragment12[, !apply(is.na(Fragment12), 2, all)]# remove columns that are all NA # If the sum of any of the rows between 6am and 7pm (ORDTimes = 360 & 1140) # it does not equal to 0 because there are missing values in all the sensors for that time frame if(any(rowSums(Fragment12[181:570,-1:-2], na.rm=T)!=0)){ DM <- degreeMinutes(timeseries = Fragment12[,-1:-2], timeresolution = 2, upperY = Vmax, lowerY = TEmerge) DM2 <- degreeMinutes8(timeseries = Fragment12[,-1:-2], timeresolution = 2, upperY = Vmax, lowerY = TEmerge, Int = intercept, Slope = slope) DMP <- degreeMinutes(timeseries = Fragment12[,-1:-2], timeresolution = 2, upperY = VmaxP, lowerY = TEmergeP) DMP2 <- degreeMinutes8(timeseries = Fragment12[,-1:-2], timeresolution = 2, upperY = VmaxP, lowerY = TEmergeP, Int = intercept, Slope = slope) Activity12[j, i+1] <- DM[[2]][1] Activity12cali[j, i+1] <- DM2[[2]][1] Restriction12[j, i+1] <- DM[[2]][2] Restriction12cali[j, i+1] <- DM2[[2]][2] Activity12P[j, i+1] <- DMP[[2]][1] Activity12caliP[j, i+1] <- DMP2[[2]][1] Restriction12P[j, i+1] <- DMP[[2]][2] Restriction12caliP[j, i+1] <- DMP2[[2]][2] } } # end of for j loop } for (k in 1:length(levels(as.factor(TeSam2013$JULDAY)))){ Fragment13 <- subset(TeTransect13, JULDAY==levels(as.factor(TeSam2013$JULDAY))[[k]]) Fragment13 <- Fragment13[, !apply(is.na(Fragment13), 2, all)]# remove columns that are all NA # If the sum of any of the rows between 6am and 7pm (ORDTimes = 360 & 1140) # it does not equal to 0 because there are missing values in all the sensors for that time frame if(any(rowSums(Fragment13[61:190,-1:-2], na.rm=T)!=0)){ # If there are data... DM <- degreeMinutes(timeseries = Fragment13[,-1:-2], timeresolution = 6, upperY = Vmax, lowerY = TEmerge) DM2 <- degreeMinutes8(timeseries = Fragment13[,-1:-2], timeresolution = 6, upperY = Vmax, lowerY = TEmerge, Int = intercept, Slope = slope) DMP <- degreeMinutes(timeseries = Fragment13[,-1:-2], timeresolution = 6, upperY = VmaxP, lowerY = TEmergeP) DMP2 <- degreeMinutes8(timeseries = Fragment13[,-1:-2], timeresolution = 6, upperY = VmaxP, lowerY = TEmergeP, Int = intercept, Slope = slope) Activity13[k, i+1] <- DM[[2]][1] Activity13cali[k, i+1] <- DM2[[2]][1] Restriction13[k, i+1] <- DM[[2]][2] Restriction13cali[k, i+1] <- DM2[[2]][2] Activity13P[k, i+1] <- DMP[[2]][1] Activity13caliP[k, i+1] <- DMP2[[2]][1] Restriction13P[k, i+1] <- DMP[[2]][2] Restriction13caliP[k, i+1] <- DMP2[[2]][2] } } # end of for k loop cat(Transect, ' ') } # end of for i loop assign(paste("Activity12", substring(species.list[[l]],1,3), sep = "."), Activity12) assign(paste("Activity12cali", substring(species.list[[l]],1,3), sep = "."), Activity12cali) assign(paste("Activity12P", substring(species.list[[l]],1,3), sep = "."), Activity12P) assign(paste("Activity12caliP", substring(species.list[[l]],1,3), sep = "."), Activity12caliP) assign(paste("Activity13", substring(species.list[[l]],1,3), sep = "."), Activity13) assign(paste("Activity13cali", substring(species.list[[l]],1,3), sep = "."), Activity13cali) assign(paste("Activity13P", substring(species.list[[l]],1,3), sep = "."), Activity13P) assign(paste("Activity13caliP", substring(species.list[[l]],1,3), sep = "."), Activity13caliP) assign(paste("Restriction12", substring(species.list[[l]],1,3), sep = "."), Restriction12) assign(paste("Restriction12cali", substring(species.list[[l]],1,3), sep = "."), Restriction12cali) assign(paste("Restriction12P", substring(species.list[[l]],1,3), sep = "."), Restriction12P) assign(paste("Restriction12caliP", substring(species.list[[l]],1,3), sep = "."), Restriction12caliP) assign(paste("Restriction13", substring(species.list[[l]],1,3), sep = "."), Restriction13) assign(paste("Restriction13cali", substring(species.list[[l]],1,3), sep = "."), Restriction13cali) assign(paste("Restriction13P", substring(species.list[[l]],1,3), sep = "."), Restriction13P) assign(paste("Restriction13caliP", substring(species.list[[l]],1,3), sep = "."), Restriction13caliP) cat(species.list[[l]], "\n") rm(Fragment12, Fragment13, i, j, k, TEmerge, Vmax, TEmergeP, VmaxP, DM, DM2, DMP, DMP2) } # end of for l loop # It is recommended to save these results after such long time of calculations # save(Restriction12.ori, Restriction12.orn, Restriction12.pun, Restriction12.sp, # Restriction12cali.ori, Restriction12cali.orn, Restriction12cali.pun, Restriction12cali.sp, # Restriction12caliP.ori, Restriction12caliP.orn, Restriction12caliP.pun, Restriction12caliP.sp, # Restriction12P.ori, Restriction12P.orn, Restriction12P.pun, Restriction12P.sp, # Restriction13.ori, Restriction13.orn, Restriction13.pun, Restriction13.sp, # Restriction13cali.ori, Restriction13cali.orn, Restriction13cali.pun, Restriction13cali.sp, Restriction13caliP, # Restriction13caliP.ori, Restriction13caliP.orn, Restriction13caliP.pun, Restriction13caliP.sp, # Restriction13P.ori, Restriction13P.orn, Restriction13P.pun, Restriction13P.sp, # Activity12.ori, Activity12.orn, Activity12.pun, Activity12.sp, # Activity12cali.ori, Activity12cali.orn, Activity12cali.pun, Activity12cali.sp, # Activity12caliP.ori, Activity12caliP.orn, Activity12caliP.pun, Activity12caliP.sp, # Activity12P.ori, Activity12P.orn, Activity12P.pun, Activity12P.sp, # Activity13.ori, Activity13.orn,Activity13.pun, Activity13.sp, # Activity13cali.ori, Activity13cali.orn, Activity13cali.pun, Activity13cali.sp, # Activity13caliP.ori, Activity13caliP.orn, Activity13caliP.pun, Activity13caliP.sp, # Activity13P.ori, Activity13P.orn, Activity13P.pun, Activity13P.sp, # file = "TemporaryBackup.RData") # # load("TemporaryBackup.RData") # We will use only the data in which we had more operative temperatures records at the same time for most of the transects, # which was during 31 days in the last fieldwork campaign (Julian Days 349 to 379: 15 December 2013 to 14 January 2014) Act.pun <- Activity13.pun[46:76,]; Act.punP <- Activity13P.pun[46:76,]; Act.puncali <- Activity13cali.pun[46:76,]; Act.puncaliP <- Activity13caliP.pun[46:76,] Res.pun <- Restriction13.pun[46:76,]; Res.punP <- Restriction13P.pun[46:76,]; Res.puncali <- Restriction13cali.pun[46:76,]; Res.puncaliP <- Restriction13caliP.pun[46:76,] Act.orn <- Activity13.orn[46:76,]; Act.ornP <- Activity13P.orn[46:76,]; Act.orncali <- Activity13cali.orn[46:76,]; Act.orncaliP <- Activity13caliP.orn[46:76,] Res.orn <- Restriction13.orn[46:76,]; Res.ornP <- Restriction13P.orn[46:76,]; Res.orncali <- Restriction13cali.orn[46:76,]; Res.orncaliP <- Restriction13caliP.orn[46:76,] Act.sp <- Activity13.sp[46:76,]; Act.spP <- Activity13P.sp[46:76,]; Act.spcali <- Activity13cali.sp[46:76,]; Act.spcaliP <- Activity13caliP.sp[46:76,] Res.sp <- Restriction13.sp[46:76,]; Res.spP <- Restriction13P.sp[46:76,]; Res.spcali <- Restriction13cali.sp[46:76,]; Res.spcaliP <- Restriction13caliP.sp[46:76,] Act.ori <- Activity13.ori[46:76,]; Act.oriP <- Activity13P.ori[46:76,]; Act.oricali <- Activity13cali.ori[46:76,]; Act.oricaliP <- Activity13caliP.ori[46:76,] Res.ori <- Restriction13.ori[46:76,]; Res.oriP <- Restriction13P.ori[46:76,]; Res.oricali <- Restriction13cali.ori[46:76,]; Res.oricaliP <- Restriction13caliP.ori[46:76,] # 2. Extrapolation of some missing values #### # In the selected period (15 December 2013 to 14 January 2014) # we are still missing the data of transects 1C, 1E and 2D # (due to unfortunate losses of the dataloggers or damage by weather conditions) # Because we have data from the 2012 campaign for those, # we can extrapolate the hours of activity and restriction for those three transects in the selected period. # We do that by linear models using the data of the closest transects in space, # and verifying they give a high explained deviance (R2>0.8): # After trying several combinations of predictors we decided that the best predictors where the data of: # · transects 1B and 1D for the data of transect 1C (R2 = 0.93) # · transects 1I and 2E for the data of transect 1E (R2 = 0.84 - 0.86) # · transects 2A, 2B and 2C for the data of transect 2D (R2 = 0.93 - 0.95) # Data for L. puna in transect 1C modelo_1C <- lm(`1C`~`1B`+`1D`, data=Activity12.pun); summary(modelo_1C) # 1C pre1C <- predict(modelo_1C, newdata=Activity13.pun[46:76,], type="response") pre1C[!is.na(pre1C) & pre1C<0] <- 0 Act.pun[,4] <- pre1C modelo_1C <- lm(`1C`~`1B`+`1D`, data=Activity12P.pun); summary(modelo_1C) # 1C P pre1C <- predict(modelo_1C, newdata=Activity13P.pun[46:76,], type="response") pre1C[!is.na(pre1C) & pre1C<0] <- 0 Act.punP[,4] <- pre1C modelo_1C <- lm(`1C`~`1B`+`1D`, data=Activity12cali.pun); summary(modelo_1C) # 1C cali pre1C <- predict(modelo_1C, newdata=Activity13cali.pun[46:76,], type="response") pre1C[!is.na(pre1C) & pre1C<0] <- 0 Act.puncali[,4] <- pre1C modelo_1C <- lm(`1C`~`1B`+`1D`, data=Activity12caliP.pun); summary(modelo_1C) # 1C caliP pre1C <- predict(modelo_1C, newdata=Activity13caliP.pun[46:76,], type="response") pre1C[!is.na(pre1C) & pre1C<0] <- 0 Act.puncaliP[,4] <- pre1C Restriction12.pun$`1C` # 1C Res Res.pun[,4] <- 0 Restriction12P.pun$`1C` # 1C Res P Res.punP[,4] <- 0 Restriction12cali.pun$`1C` # 1C Res cali Res.puncali[,4] <- 0 Restriction12caliP.pun$`1C` # 1C Res caliP Res.puncaliP[,4] <- 0 # Data for L. puna in transect 1E modelo_1E <- lm(`1E`~`1I`+`2E`, data=Activity12.pun); summary(modelo_1E) # 1E pre1E <- predict(modelo_1E, newdata=Activity13.pun[46:76,], type="response") pre1E[!is.na(pre1E) & pre1E<0] <- 0 Act.pun[,6] <- pre1E modelo_1E <- lm(`1E`~`1I`+`2E`, data=Activity12P.pun); summary(modelo_1E) # 1E P pre1E <- predict(modelo_1E, newdata=Activity13P.pun[46:76,], type="response") pre1E[!is.na(pre1E) & pre1E<0] <- 0 Act.punP[,6] <- pre1E modelo_1E <- lm(`1E`~`1I`+`2E`, data=Activity12cali.pun); summary(modelo_1E) # 1E cali pre1E <- predict(modelo_1E, newdata=Activity13cali.pun[46:76,], type="response") pre1E[!is.na(pre1E) & pre1E<0] <- 0 Act.puncali[,6] <- pre1E modelo_1E <- lm(`1E`~`1I`+`2E`, data=Activity12caliP.pun); summary(modelo_1E) # 1E caliP pre1E <- predict(modelo_1E, newdata=Activity13caliP.pun[46:76,], type="response") pre1E[!is.na(pre1E) & pre1E<0] <- 0 Act.puncaliP[,6] <- pre1E Restriction12.pun$`1E` # 1E Res Res.pun[,6] <- 0 Restriction12P.pun$`1E` # 1E Res P Res.punP[,6] <- 0 Restriction12cali.pun$`1E` # 1E Res cali Res.puncali[,6] <- 0 Restriction12caliP.pun$`1E` # 1E Res caliP Res.puncaliP[,6] <- 0 # Data for L. puna in transect 2D modelo_2D <- lm(`2D`~`2A`+`2B`+`2C`, data=Activity12.pun); summary(modelo_2D) # 2D pre2D <- predict(modelo_2D, newdata=Activity13.pun[46:76,], type="response") pre2D[!is.na(pre2D) & pre2D<0] <- 0 Act.pun[,15] <- pre2D modelo_2D <- lm(`2D`~`2A`+`2B`+`2C`, data=Activity12P.pun); summary(modelo_2D) # 2D P pre2D <- predict(modelo_2D, newdata=Activity13P.pun[46:76,], type="response") pre2D[!is.na(pre2D) & pre2D<0] <- 0 Act.punP[,15] <- pre2D modelo_2D <- lm(`2D`~`2A`+`2B`+`2C`, data=Activity12cali.pun); summary(modelo_2D) # 2D cali pre2D <- predict(modelo_2D, newdata=Activity13cali.pun[46:76,], type="response") pre2D[!is.na(pre2D) & pre2D<0] <- 0 Act.puncali[,15] <- pre2D modelo_2D <- lm(`2D`~`2A`+`2B`+`2C`, data=Activity12caliP.pun); summary(modelo_2D) # 2D caliP pre2D <- predict(modelo_2D, newdata=Activity13caliP.pun[46:76,], type="response") pre2D[!is.na(pre2D) & pre2D<0] <- 0 Act.puncaliP[,15] <- pre2D Restriction12.pun$`2D` # 2D Res Res.pun[,15] <- 0 Restriction12P.pun$`2D` # 2D Res P Res.punP[,15] <- 0 Restriction12cali.pun$`2D` # 2D Res cali Res.puncali[,15] <- 0 Restriction12caliP.pun$`2D` # 2D Res caliP Res.puncaliP[,15] <- 0 # Data for L. ornatus in transect 1C modelo_1C <- lm(`1C`~`1B`+`1D`, data=Activity12.orn); summary(modelo_1C) # 1C pre1C <- predict(modelo_1C, newdata=Activity13.orn[46:76,], type="response") pre1C[!is.na(pre1C) & pre1C<0] <- 0 Act.orn[,4] <- pre1C modelo_1C <- lm(`1C`~`1B`+`1D`, data=Activity12P.orn); summary(modelo_1C) # 1C P pre1C <- predict(modelo_1C, newdata=Activity13P.orn[46:76,], type="response") pre1C[!is.na(pre1C) & pre1C<0] <- 0 Act.ornP[,4] <- pre1C modelo_1C <- lm(`1C`~`1B`+`1D`, data=Activity12cali.orn); summary(modelo_1C) # 1C cali pre1C <- predict(modelo_1C, newdata=Activity13cali.orn[46:76,], type="response") pre1C[!is.na(pre1C) & pre1C<0] <- 0 Act.orncali[,4] <- pre1C modelo_1C <- lm(`1C`~`1B`+`1D`, data=Activity12caliP.orn); summary(modelo_1C) # 1C caliP pre1C <- predict(modelo_1C, newdata=Activity13caliP.orn[46:76,], type="response") pre1C[!is.na(pre1C) & pre1C<0] <- 0 Act.orncaliP[,4] <- pre1C Restriction12.orn$`1C` # 1C Res Res.orn[,4] <- 0 Restriction12P.orn$`1C` # 1C Res P Res.ornP[,4] <- 0 Restriction12cali.orn$`1C` # 1C Res cali Res.orncali[,4] <- 0 Restriction12caliP.orn$`1C` # 1C Res caliP Res.orncaliP[,4] <- 0 # Data for L. ornatus in transect 1E modelo_1E <- lm(`1E`~`1I`+`2E`, data=Activity12.orn); summary(modelo_1E) # 1E pre1E <- predict(modelo_1E, newdata=Activity13.orn[46:76,], type="response") pre1E[!is.na(pre1E) & pre1E<0] <- 0 Act.orn[,6] <- pre1E modelo_1E <- lm(`1E`~`1I`+`2E`, data=Activity12P.orn); summary(modelo_1E) # 1E P pre1E <- predict(modelo_1E, newdata=Activity13P.orn[46:76,], type="response") pre1E[!is.na(pre1E) & pre1E<0] <- 0 Act.ornP[,6] <- pre1E modelo_1E <- lm(`1E`~`1I`+`2E`, data=Activity12cali.orn); summary(modelo_1E) # 1E cali pre1E <- predict(modelo_1E, newdata=Activity13cali.orn[46:76,], type="response") pre1E[!is.na(pre1E) & pre1E<0] <- 0 Act.orncali[,6] <- pre1E modelo_1E <- lm(`1E`~`1I`+`2E`, data=Activity12caliP.orn); summary(modelo_1E) # 1E caliP pre1E <- predict(modelo_1E, newdata=Activity13caliP.orn[46:76,], type="response") pre1E[!is.na(pre1E) & pre1E<0] <- 0 Act.orncaliP[,6] <- pre1E Restriction12.orn$`1E` # 1E Res Res.orn[,6] <- 0 Restriction12P.orn$`1E` # 1E Res P Res.ornP[,6] <- 0 Restriction12cali.orn$`1E` # 1E Res cali Res.orncali[,6] <- 0 Restriction12caliP.orn$`1E` # 1E Res caliP Res.orncaliP[,6] <- 0 # Data for L. ornatus in transect 2D modelo_2D <- lm(`2D`~`2A`+`2B`+`2C`, data=Activity12.orn); summary(modelo_2D) # 2D pre2D <- predict(modelo_2D, newdata=Activity13.orn[46:76,], type="response") pre2D[!is.na(pre2D) & pre2D<0] <- 0 Act.orn[,15] <- pre2D modelo_2D <- lm(`2D`~`2A`+`2B`+`2C`, data=Activity12P.orn); summary(modelo_2D) # 2D P pre2D <- predict(modelo_2D, newdata=Activity13P.orn[46:76,], type="response") pre2D[!is.na(pre2D) & pre2D<0] <- 0 Act.ornP[,15] <- pre2D modelo_2D <- lm(`2D`~`2A`+`2B`+`2C`, data=Activity12cali.orn); summary(modelo_2D) # 2D cali pre2D <- predict(modelo_2D, newdata=Activity13cali.orn[46:76,], type="response") pre2D[!is.na(pre2D) & pre2D<0] <- 0 Act.orncali[,15] <- pre2D modelo_2D <- lm(`2D`~`2A`+`2B`+`2C`, data=Activity12caliP.orn); summary(modelo_2D) # 2D caliP pre2D <- predict(modelo_2D, newdata=Activity13caliP.orn[46:76,], type="response") pre2D[!is.na(pre2D) & pre2D<0] <- 0 Act.orncaliP[,15] <- pre2D Restriction12.orn$`2D` # 2D Res Res.orn[,15] <- 0 Restriction12P.orn$`2D` # 2D Res P Res.ornP[,15] <- 0 Restriction12cali.orn$`2D` # 2D Res cali Res.orncali[,15] <- 0 Restriction12caliP.orn$`2D` # 2D Res caliP Res.orncaliP[,15] <- 0 # Data for Liolaemus sp. in transect 1C modelo_1C <- lm(`1C`~`1B`+`1D`, data=Activity12.sp); summary(modelo_1C) # 1C pre1C <- predict(modelo_1C, newdata=Activity13.sp[46:76,], type="response") pre1C[!is.na(pre1C) & pre1C<0] <- 0 Act.sp[,4] <- pre1C modelo_1C <- lm(`1C`~`1B`+`1D`, data=Activity12P.sp); summary(modelo_1C) # 1C P pre1C <- predict(modelo_1C, newdata=Activity13P.sp[46:76,], type="response") pre1C[!is.na(pre1C) & pre1C<0] <- 0 Act.spP[,4] <- pre1C modelo_1C <- lm(`1C`~`1B`+`1D`, data=Activity12cali.sp); summary(modelo_1C) # 1C cali pre1C <- predict(modelo_1C, newdata=Activity13cali.sp[46:76,], type="response") pre1C[!is.na(pre1C) & pre1C<0] <- 0 Act.spcali[,4] <- pre1C modelo_1C <- lm(`1C`~`1B`+`1D`, data=Activity12caliP.sp); summary(modelo_1C) # 1C caliP pre1C <- predict(modelo_1C, newdata=Activity13caliP.sp[46:76,], type="response") pre1C[!is.na(pre1C) & pre1C<0] <- 0 Act.spcaliP[,4] <- pre1C Restriction12.sp$`1C` # 1C Res Res.sp[,4] <- 0 Restriction12P.sp$`1C` # 1C Res P Res.spP[,4] <- 0 Restriction12cali.sp$`1C` # 1C Res cali Res.spcali[,4] <- 0 Restriction12caliP.sp$`1C` # 1C Res caliP Res.spcaliP[,4] <- 0 # Data for Liolaemus sp. in transect 1E modelo_1E <- lm(`1E`~`1I`+`2E`, data=Activity12.sp); summary(modelo_1E) # 1E pre1E <- predict(modelo_1E, newdata=Activity13.sp[46:76,], type="response") pre1E[!is.na(pre1E) & pre1E<0] <- 0 Act.sp[,6] <- pre1E modelo_1E <- lm(`1E`~`1I`+`2E`, data=Activity12P.sp); summary(modelo_1E) # 1E P pre1E <- predict(modelo_1E, newdata=Activity13P.sp[46:76,], type="response") pre1E[!is.na(pre1E) & pre1E<0] <- 0 Act.spP[,6] <- pre1E modelo_1E <- lm(`1E`~`1I`+`2E`, data=Activity12cali.sp); summary(modelo_1E) # 1E cali pre1E <- predict(modelo_1E, newdata=Activity13cali.sp[46:76,], type="response") pre1E[!is.na(pre1E) & pre1E<0] <- 0 Act.spcali[,6] <- pre1E modelo_1E <- lm(`1E`~`1I`+`2E`, data=Activity12caliP.sp); summary(modelo_1E) # 1E caliP pre1E <- predict(modelo_1E, newdata=Activity13caliP.sp[46:76,], type="response") pre1E[!is.na(pre1E) & pre1E<0] <- 0 Act.spcaliP[,6] <- pre1E Restriction12.sp$`1E` # 1E Res Res.sp[,6] <- 0 Restriction12P.sp$`1E` # 1E Res P Res.spP[,6] <- 0 Restriction12cali.sp$`1E` # 1E Res cali Res.spcali[,6] <- 0 Restriction12caliP.sp$`1E` # 1E Res caliP Res.spcaliP[,6] <- 0 # Data for Liolaemus sp. in transect 2D modelo_2D <- lm(`2D`~`2A`+`2B`+`2C`, data=Activity12.sp); summary(modelo_2D) # 2D pre2D <- predict(modelo_2D, newdata=Activity13.sp[46:76,], type="response") pre2D[!is.na(pre2D) & pre2D<0] <- 0 Act.sp[,15] <- pre2D modelo_2D <- lm(`2D`~`2A`+`2B`+`2C`, data=Activity12P.sp); summary(modelo_2D) # 2D P pre2D <- predict(modelo_2D, newdata=Activity13P.sp[46:76,], type="response") pre2D[!is.na(pre2D) & pre2D<0] <- 0 Act.spP[,15] <- pre2D modelo_2D <- lm(`2D`~`2A`+`2B`+`2C`, data=Activity12cali.sp); summary(modelo_2D) # 2D cali pre2D <- predict(modelo_2D, newdata=Activity13cali.sp[46:76,], type="response") pre2D[!is.na(pre2D) & pre2D<0] <- 0 Act.spcali[,15] <- pre2D modelo_2D <- lm(`2D`~`2A`+`2B`+`2C`, data=Activity12caliP.sp); summary(modelo_2D) # 2D caliP pre2D <- predict(modelo_2D, newdata=Activity13caliP.sp[46:76,], type="response") pre2D[!is.na(pre2D) & pre2D<0] <- 0 Act.spcaliP[,15] <- pre2D Restriction12.sp$`2D` # 2D Res Res.sp[,15] <- 0 Restriction12P.sp$`2D` # 2D Res P Res.spP[,15] <- 0 Restriction12cali.sp$`2D` # 2D Res cali Res.spcali[,15] <- 0 Restriction12caliP.sp$`2D` # 2D Res caliP Res.spcaliP[,15] <- 0 # Data for L. orientalis in transect 1C modelo_1C <- lm(`1C`~`1B`+`1D`, data=Activity12.ori); summary(modelo_1C) # 1C pre1C <- predict(modelo_1C, newdata=Activity13.ori[46:76,], type="response") pre1C[!is.na(pre1C) & pre1C<0] <- 0 Act.ori[,4] <- pre1C modelo_1C <- lm(`1C`~`1B`+`1D`, data=Activity12P.ori); summary(modelo_1C) # 1C P pre1C <- predict(modelo_1C, newdata=Activity13P.ori[46:76,], type="response") pre1C[!is.na(pre1C) & pre1C<0] <- 0 Act.oriP[,4] <- pre1C modelo_1C <- lm(`1C`~`1B`+`1D`, data=Activity12cali.ori); summary(modelo_1C) # 1C cali pre1C <- predict(modelo_1C, newdata=Activity13cali.ori[46:76,], type="response") pre1C[!is.na(pre1C) & pre1C<0] <- 0 Act.oricali[,4] <- pre1C modelo_1C <- lm(`1C`~`1B`+`1D`, data=Activity12caliP.ori); summary(modelo_1C) # 1C caliP pre1C <- predict(modelo_1C, newdata=Activity13caliP.ori[46:76,], type="response") pre1C[!is.na(pre1C) & pre1C<0] <- 0 Act.oricaliP[,4] <- pre1C Restriction12.ori$`1C` # 1C Res Res.ori[,4] <- 0 Restriction12P.ori$`1C` # 1C Res P Res.oriP[,4] <- 0 Restriction12cali.ori$`1C` # 1C Res cali Res.oricali[,4] <- 0 Restriction12caliP.ori$`1C` # 1C Res caliP Res.oricaliP[,4] <- 0 # Data for L. orientalis in transect 1E modelo_1E <- lm(`1E`~`1I`+`2E`, data=Activity12.ori); summary(modelo_1E) # 1E pre1E <- predict(modelo_1E, newdata=Activity13.ori[46:76,], type="response") pre1E[!is.na(pre1E) & pre1E<0] <- 0 Act.ori[,6] <- pre1E modelo_1E <- lm(`1E`~`1I`+`2E`, data=Activity12P.ori); summary(modelo_1E) # 1E P pre1E <- predict(modelo_1E, newdata=Activity13P.ori[46:76,], type="response") pre1E[!is.na(pre1E) & pre1E<0] <- 0 Act.oriP[,6] <- pre1E modelo_1E <- lm(`1E`~`1I`+`2E`, data=Activity12cali.ori); summary(modelo_1E) # 1E cali pre1E <- predict(modelo_1E, newdata=Activity13cali.ori[46:76,], type="response") pre1E[!is.na(pre1E) & pre1E<0] <- 0 Act.oricali[,6] <- pre1E modelo_1E <- lm(`1E`~`1I`+`2E`, data=Activity12caliP.ori); summary(modelo_1E) # 1E caliP pre1E <- predict(modelo_1E, newdata=Activity13caliP.ori[46:76,], type="response") pre1E[!is.na(pre1E) & pre1E<0] <- 0 Act.oricaliP[,6] <- pre1E Restriction12.ori$`1E` # 1E Res Res.ori[,6] <- 0 Restriction12P.ori$`1E` # 1E Res P Res.oriP[,6] <- 0 Restriction12cali.ori$`1E` # 1E Res cali Res.oricali[,6] <- 0 Restriction12caliP.ori$`1E` # 1E Res caliP Res.oricaliP[,6] <- 0 # Data for L. orientalis in transect 2D modelo_2D <- lm(`2D`~`2A`+`2B`+`2C`, data=Activity12.ori); summary(modelo_2D) # 2D pre2D <- predict(modelo_2D, newdata=Activity13.ori[46:76,], type="response") pre2D[!is.na(pre2D) & pre2D<0] <- 0 Act.ori[,15] <- pre2D modelo_2D <- lm(`2D`~`2A`+`2B`+`2C`, data=Activity12P.ori); summary(modelo_2D) # 2D P pre2D <- predict(modelo_2D, newdata=Activity13P.ori[46:76,], type="response") pre2D[!is.na(pre2D) & pre2D<0] <- 0 Act.oriP[,15] <- pre2D modelo_2D <- lm(`2D`~`2A`+`2B`+`2C`, data=Activity12cali.ori); summary(modelo_2D) # 2D cali pre2D <- predict(modelo_2D, newdata=Activity13cali.ori[46:76,], type="response") pre2D[!is.na(pre2D) & pre2D<0] <- 0 Act.oricali[,15] <- pre2D modelo_2D <- lm(`2D`~`2A`+`2B`+`2C`, data=Activity12caliP.ori); summary(modelo_2D) # 2D caliP pre2D <- predict(modelo_2D, newdata=Activity13caliP.ori[46:76,], type="response") pre2D[!is.na(pre2D) & pre2D<0] <- 0 Act.oricaliP[,15] <- pre2D Restriction12.ori$`2D` # 2D Res Res.ori[,15] <- 0 Restriction12P.ori$`2D` # 2D Res P Res.oriP[,15] <- 0 Restriction12cali.ori$`2D` # 2D Res cali Res.oricali[,15] <- 0 Restriction12caliP.ori$`2D` # 2D Res caliP Res.oricaliP[,15] <- 0 # 3. Plot of activity hours per transect #### par(mfrow=c(1,1), mar=c(5,5,.5,.5)) boxplot(Act.pun[,-1]/60, at=seq(1, 150, by=5), col="blue4", border=rgb(0,0,0,.5), xaxt='n', ylim=c(0,11), ylab='Average potential activity (h)', xlab='Transects') abline(v=seq(7.5, 150, by=10), lwd=30, col=rgb(0,0,0,.25)) abline(v=c(50,100), lwd=2) boxplot(Act.pun[,-1]/60, at=seq(1, 150, by=5), col="blue4", border=rgb(0,0,0,.5), xaxt='n', yaxt='n', add=T) # Solo datos del periodo de solapamiento boxplot(Act.orn[,-1]/60, at=seq(2, 150, by=5), col="orange", border=rgb(0,0,0,.5), xaxt='n', yaxt='n', add=T) boxplot(Act.sp[,-1]/60, at=seq(3, 150, by=5), col="yellow", border=rgb(0,0,0,.5), xaxt='n', yaxt='n', add=T) boxplot(Act.ori[,-1]/60, at=seq(4, 150, by=5), col="forestgreen", border=rgb(0,0,0,.5), xaxt='n', yaxt='n', add=T) axis(at=seq(2.5, 150, by=5), labels=Transectslist, side=1) axis(1, at=c(25,75,125), pos=-1.2, labels=c("High (4170–4300 m)", "Intermediate (3970–4010 m )", "Low (3740–3780 m m)"), tck=0, lwd=0) points(y=colMeans(Act.pun[,-1], na.rm=T)/60, x=seq(1, 150, by=5), pch="-", col='red') points(y=colMeans(Act.orn[,-1], na.rm=T)/60, x=seq(2, 150, by=5), pch="-", col='red') points(y=colMeans(Act.sp[,-1], na.rm=T)/60, x=seq(3, 150, by=5), pch="-", col='red') points(y=colMeans(Act.ori[,-1], na.rm=T)/60, x=seq(4, 150, by=5), pch="-", col='red') par(mfrow=c(1,1), mar=c(5,5,.5,.5)) boxplot(Act.punP[,-1]/60, at=seq(1, 150, by=5), col="blue4", border=rgb(0,0,0,.5), xaxt='n', ylim=c(0,11), ylab='Average potential activity (h)', xlab='Transects') abline(v=seq(7.5, 150, by=10), lwd=30, col=rgb(0,0,0,.25)) abline(v=c(50,100), lwd=2) boxplot(Act.punP[,-1]/60, at=seq(1, 150, by=5), col="blue4", border=rgb(0,0,0,.5), xaxt='n', yaxt='n', add=T) # Solo datos del periodo de solapamiento boxplot(Act.ornP[,-1]/60, at=seq(2, 150, by=5), col="orange", border=rgb(0,0,0,.5), xaxt='n', yaxt='n', add=T) boxplot(Act.spP[,-1]/60, at=seq(3, 150, by=5), col="yellow", border=rgb(0,0,0,.5), xaxt='n', yaxt='n', add=T) boxplot(Act.oriP[,-1]/60, at=seq(4, 150, by=5), col="forestgreen", border=rgb(0,0,0,.5), xaxt='n', yaxt='n', add=T) axis(at=seq(2.5, 150, by=5), labels=Transectslist, side=1) axis(1, at=c(25,75,125), pos=-1.2, labels=c("High (4170–4300 m)", "Intermediate (3970–4010 m )", "Low (3740–3780 m m)"), tck=0, lwd=0) points(y=colMeans(Act.punP[,-1], na.rm=T)/60, x=seq(1, 150, by=5), pch="-", col='red') points(y=colMeans(Act.ornP[,-1], na.rm=T)/60, x=seq(2, 150, by=5), pch="-", col='red') points(y=colMeans(Act.spP[,-1], na.rm=T)/60, x=seq(3, 150, by=5), pch="-", col='red') points(y=colMeans(Act.oriP[,-1], na.rm=T)/60, x=seq(4, 150, by=5), pch="-", col='red') par(mfrow=c(1,1), mar=c(5,5,.5,.5)) boxplot(Act.puncali[,-1]/60, at=seq(1, 150, by=5), col="blue4", border=rgb(0,0,0,.5), xaxt='n', ylim=c(0,11), ylab='Average potential activity (h)', xlab='Transects') abline(v=seq(7.5, 150, by=10), lwd=30, col=rgb(0,0,0,.25)) abline(v=c(50,100), lwd=2) boxplot(Act.puncali[,-1]/60, at=seq(1, 150, by=5), col="blue4", border=rgb(0,0,0,.5), xaxt='n', yaxt='n', add=T) # Solo datos del periodo de solapamiento boxplot(Act.orncali[,-1]/60, at=seq(2, 150, by=5), col="orange", border=rgb(0,0,0,.5), xaxt='n', yaxt='n', add=T) boxplot(Act.spcali[,-1]/60, at=seq(3, 150, by=5), col="yellow", border=rgb(0,0,0,.5), xaxt='n', yaxt='n', add=T) boxplot(Act.oricali[,-1]/60, at=seq(4, 150, by=5), col="forestgreen", border=rgb(0,0,0,.5), xaxt='n', yaxt='n', add=T) axis(at=seq(2.5, 150, by=5), labels=Transectslist, side=1) axis(1, at=c(25,75,125), pos=-1.2, labels=c("High (4170–4300 m)", "Intermediate (3970–4010 m )", "Low (3740–3780 m m)"), tck=0, lwd=0) points(y=colMeans(Act.puncali[,-1], na.rm=T)/60, x=seq(1, 150, by=5), pch="-", col='red') points(y=colMeans(Act.orncali[,-1], na.rm=T)/60, x=seq(2, 150, by=5), pch="-", col='red') points(y=colMeans(Act.spcali[,-1], na.rm=T)/60, x=seq(3, 150, by=5), pch="-", col='red') points(y=colMeans(Act.oricali[,-1], na.rm=T)/60, x=seq(4, 150, by=5), pch="-", col='red') # 4. Saving the results #### save(Activity12.pun, Activity13.pun, Activity12.orn, Activity13.orn, Activity12.sp, Activity13.sp, Activity12.ori, Activity13.ori, file= "output_data/SummaryActivities_Sama.RData") save(Activity12P.pun, Activity13P.pun, Activity12P.orn, Activity13P.orn, Activity12P.sp, Activity13P.sp, Activity12P.ori, Activity13P.ori, file= "output_data/SummaryPActivities_Sama.RData") save(Activity12cali.pun, Activity13cali.pun, Activity12cali.orn, Activity13cali.orn, Activity12cali.sp, Activity13cali.sp, Activity12cali.ori, Activity13cali.ori, file= "output_data/SummaryCalibratedActivities_Sama.RData") save(Activity12caliP.pun, Activity13caliP.pun, Activity12caliP.orn, Activity13caliP.orn, Activity12caliP.sp, Activity13caliP.sp, Activity12caliP.ori, Activity13caliP.ori, file= "output_data/SummaryCalibratedPActivities_Sama.RData") save(Restriction12.pun, Restriction13.pun, Restriction12.orn, Restriction13.orn, Restriction12.sp, Restriction13.sp, Restriction12.ori, Restriction13.ori, file= "output_data/SummaryRestrictions_Sama.RData") save(Restriction12P.pun, Restriction13P.pun, Restriction12P.orn, Restriction13P.orn, Restriction12P.sp, Restriction13P.sp, Restriction12P.ori, Restriction13P.ori, file= "output_data/SummaryPRestrictions_Sama.RData") save(Restriction12cali.pun, Restriction13cali.pun, Restriction12cali.orn, Restriction13cali.orn, Restriction12cali.sp, Restriction13cali.sp, Restriction12cali.ori, Restriction13cali.ori, file= "output_data/SummaryCalibratedRestrictions_Sama.RData") save(Restriction12caliP.pun, Restriction13caliP.pun, Restriction12caliP.orn, Restriction13caliP.orn, Restriction12caliP.sp, Restriction13caliP.sp, Restriction12caliP.ori, Restriction13caliP.ori, file= "output_data/SummaryCalibratedPRestrictions_Sama.RData") DegreeMinutesTable <- data.frame(Transect=Transectslist, Act.pun=(colMeans(Act.pun[,2:31], na.rm=T)), Act.orn=(colMeans(Act.orn[,2:31], na.rm=T)), Act.ori=(colMeans(Act.ori[,2:31], na.rm=T)), Act.sp=(colMeans(Act.sp[,2:31], na.rm=T)), Act.punP=(colMeans(Act.punP[,2:31], na.rm=T)), Act.ornP=(colMeans(Act.ornP[,2:31], na.rm=T)), Act.oriP=(colMeans(Act.oriP[,2:31], na.rm=T)), Act.spP=(colMeans(Act.spP[,2:31], na.rm=T)), Act.puncali=(colMeans(Act.puncali[,2:31], na.rm=T)), Act.orncali=(colMeans(Act.orncali[,2:31], na.rm=T)), Act.oricali=(colMeans(Act.oricali[,2:31], na.rm=T)), Act.spcali=(colMeans(Act.spcali[,2:31], na.rm=T)), Act.puncaliP=(colMeans(Act.puncaliP[,2:31], na.rm=T)), Act.orncaliP=(colMeans(Act.orncaliP[,2:31], na.rm=T)), Act.oricaliP=(colMeans(Act.oricaliP[,2:31], na.rm=T)), Act.spcaliP=(colMeans(Act.spcaliP[,2:31], na.rm=T)), Res.pun=(colMeans(Res.pun[,2:31], na.rm=T)), Res.orn=(colMeans(Res.orn[,2:31], na.rm=T)), Res.ori=(colMeans(Res.ori[,2:31], na.rm=T)), Res.sp=(colMeans(Res.sp[,2:31], na.rm=T)), Res.punP=(colMeans(Res.punP[,2:31], na.rm=T)), Res.ornP=(colMeans(Res.ornP[,2:31], na.rm=T)), Res.oriP=(colMeans(Res.oriP[,2:31], na.rm=T)), Res.spP=(colMeans(Res.spP[,2:31], na.rm=T)), Res.puncali=(colMeans(Res.puncali[,2:31], na.rm=T)), Res.orncali=(colMeans(Res.orncali[,2:31], na.rm=T)), Res.oricali=(colMeans(Res.oricali[,2:31], na.rm=T)), Res.spcali=(colMeans(Res.spcali[,2:31], na.rm=T)), Res.puncaliP=(colMeans(Res.puncaliP[,2:31], na.rm=T)), Res.orncaliP=(colMeans(Res.orncaliP[,2:31], na.rm=T)), Res.oricaliP=(colMeans(Res.oricaliP[,2:31], na.rm=T)), Res.spcaliP=(colMeans(Res.spcaliP[,2:31], na.rm=T))) write.table(DegreeMinutesTable, file="output_data/ActivityPerTransect_Sama.csv", sep=";")