finalize fish indicators for TCRMP and PRCRMP

indicator_processing/non_automated/PRCRMP_benthic.R

@@ -158,6 +158,9 @@ hist(ben$corspprich)
 # merge site reference info with database --------------------------
 
 ben2 <- merge(ben, met, by.x = "SITE.NAME", by.y = "Site.Name")
+dim(ben)
+dim(met)
+dim(ben2)
 table(ben2$SITE.NAME, ben2$YEAR)
 table(ben2$loc, ben2$YEAR)
 
@@ -208,7 +211,7 @@ mod <- summary(out1)$coef[1:length(yrs), 1]
 modse <- summary(out1)$coef[1:length(yrs), 2]
 
 par(mar = c(4, 4, 3, 1))
-plot(yrs, ind, ylim = c(min(ind) - 4, max(ind + 2)), main = varint)
+plot(yrs, ind, ylim = c(min(ind) * 0.6, max(ind) * 1.1), main = varint)
 lines(yrs, ind + indse, lty = 2, col = 1)
 lines(yrs, ind - indse, lty = 2, col = 1)
 

indicator_processing/non_automated/PRCRMP_fish.R

@@ -54,16 +54,16 @@ sitelis
 tail(fish$YEAR, 20)
 #fish <- fish[-which(is.na(fish$YEAR)), ]
 
-table(fish$LOCATION)
-table(fish$DEPTH.ZONE)
+table(fish$LOCATION, useNA = "always")
+table(fish$DEPTH.ZONE, useNA = "always")
 
-#fish$LOCATION[which(fish$LOCATION == "Mayagüez")] <- "Mayaguez"
+fish$LOCATION[which(fish$LOCATION == "Mayagüez")] <- "Mayaguez"
 fish$DEPTH.ZONE[grep("mediate", fish$DEPTH.ZONE)] <- "intermediate"
 fish$DEPTH.ZONE[grep("photic", tolower(fish$DEPTH.ZONE))] <- "mesophotic"
 fish$DEPTH.ZONE[grep("ery", fish$DEPTH.ZONE)] <- "very shal"
 fish$DEPTH.ZONE[grep("hallow", fish$DEPTH.ZONE)] <- "shallow"
 
-apply(fish[1:10], 2, table)
+apply(fish[1:10], 2, table, useNA = "always")
 
 num <- NA
 for (i in 1:ncol(fish)) { num[i] <- length(which(fish[, i] == ""))  } 
@@ -99,8 +99,9 @@ table(is.na(siz$YEAR))
 #siz <- siz[-which(is.na(siz$YEAR)), ]
 
 # check labels and correct ------------------------------------------
-table(siz$LOCATION)
-table(siz$DEPTH.ZONE)
+
+table(siz$LOCATION, useNA = "always")
+table(siz$DEPTH.ZONE, useNA = "always")
 siz$LOCATION[which(siz$LOCATION == "Mayagüez")] <- "Mayaguez"
 siz$DEPTH.ZONE[grep("mediate", siz$DEPTH.ZONE)] <- "intermediate"
 siz$DEPTH.ZONE[grep("photic", tolower(siz$DEPTH.ZONE))] <- "mesophotic"
@@ -144,10 +145,11 @@ splis2
 
 yrs <- sort(unique(fish$YEAR))
 
-apply(fish[1:10], 2, table)
+apply(fish[1:10], 2, table, useNA = "always")
 
 names(fish)
 names(fish)[which(names(fish) == "Abudefduf.saxatilis"): (ncol(fish))]
+table(fish$TRANSECT, useNA = "always")
 
 # run separately for fished and not fished -------------------------------
 
@@ -157,16 +159,18 @@ fish2 <- fish[which(names(fish) %in% splis2)]; grp <- "com"  # fished
 names(fish2)
 
 # calculate total fish density by site ---------------------------------
+
 fish$dens <- NA
 for (i in 1:nrow(fish2)) {  fish$dens[i] <- sum(as.numeric(fish2[i, ]))   }
 hist(fish$dens)
 table(fish$dens, useNA = "always")
 
 table(fish$SITE.NAME, fish$YEAR)
+image(table(fish$YEAR, fish$SITE.NAME))
 
 # analysis of average fish density ------------------------------------
 
-fish$var <- log(fish$dens+0.001)  
+fish$var <- log(fish$dens+0.1)  
 
 hist(fish$var)
 
@@ -205,59 +209,137 @@ lines(yrs, mod - modse, col = 2, lty = 2)
 
 cor(ind, mod)
 
-save(out1, file = "C:/Users/mandy.karnauskas/Desktop/Caribbean-ESR/indicator_data/PRCRMP/fish_density_PR.RData")
+# functions for converting lognormal mean and SE to normal space -------------
+lnorm.mean <- function(x1, x1e) {  exp(x1 + 0.5 * x1e^2)   }
+lnorm.se   <- function(x1, x1e) {  ((exp(x1e^2)-1)*exp(2 * x1 + x1e^2))^0.5  } 
+
+ind_norm   <- lnorm.mean(mod, modse)                # convert lognormal mean and SE to normal space
+indse_norm <- lnorm.se(mod, modse)
+
+plot(yrs, (ind_norm), ylim = c(0, 3))
+lines(yrs, (ind_norm + indse_norm), lty = 2, col = 1)
+lines(yrs, (ind_norm - indse_norm), lty = 2, col = 1)
 
-############ NOT FINALIZED BELOW 
-############ COPIED OVER FROM TRCRMP BUT NEED TO MODIFY FOR PR DATA
+points(yrs, tapply(fish$dens, fish$YEAR, mean, na.rm = T), col = 2)
+
+fin <- data.frame(cbind(yrs, ind_norm, indse_norm))
+fin
+
+save(fin, file = "C:/Users/mandy.karnauskas/Desktop/Caribbean-ESR/indicator_data/PRCRMP/fish_density_PR.RData")
+
+###########################  END  ###############################
+
+##############   CODE BELOW IS ABORTED   #######################
+# CODE FUNCTIONS BUT NOT SUFFICIENT DATA TO RUN ANALYSIS
 # size-based indicators -------------------------------------
 
-# set up empty matrix to put fish length data in ------------
-v1 <- rep(sitelis, length(yrs))
-v2 <- sort(rep(yrs, length(sitelis)))
-mat <- matrix(data = NA, nrow = length(v1), ncol = length(names(fish)[11:21]))
+head(siz)
+splis2
+
+# get extract columns for commercial species -----------------
+cols <- c()
+for (i in splis2) { 
+  co <- grep(i, names(siz))
+  if (length(co) > 0)  { cols <- c(cols, co) }
+  }
+
+cols <- sort(cols)
+names(siz)[cols]
+dim(siz)
+siz <- siz[, c(1:7, cols)]
+dim(siz)
+head(siz)
+
+yrs <- yrs[which(yrs >= min(siz$YEAR))]
+
+# reformat size data into different object --------------------
+
+#v1 <- rep(sitelis, length(yrs))
+#v2 <- sort(rep(yrs, length(sitelis)))
+v2 <- yrs
+v1 <- NA
+mat <- matrix(data = NA, nrow = length(v1), ncol = 12)
 mat <- data.frame(v1, v2, mat)
-names(mat) <- c("site", "yr", "X0.10", names(fish)[13:21], "total")
+names(mat) <- c("site", "yr", paste0("c", 1:12))
 head(mat)
 
-# summarize by site and year
-i <- 0
-for (j in sitelis) { 
+head(siz)
+names(siz)[grep(".c10", names(siz))] <- gsub(".c10", ".cX10", names(siz)[grep(".c10", names(siz))])
+names(siz)[grep(".c11", names(siz))] <- gsub(".c11", ".cX11", names(siz)[grep(".c11", names(siz))])
+names(siz)[grep(".c12", names(siz))] <- gsub(".c12", ".cX12", names(siz)[grep(".c12", names(siz))])
+names(siz)[grep(".c13", names(siz))] <- gsub(".c13", ".cX13", names(siz)[grep(".c13", names(siz))])
+names(siz)[grep(".c14", names(siz))] <- gsub(".c14", ".cX14", names(siz)[grep(".c14", names(siz))])
+names(siz)[grep(".c15", names(siz))] <- gsub(".c15", ".cX15", names(siz)[grep(".c15", names(siz))])
+names(siz)[grep(".c16", names(siz))] <- gsub(".c16", ".cX16", names(siz)[grep(".c16", names(siz))])
+names(siz)[grep(".c17", names(siz))] <- gsub(".c17", ".cX17", names(siz)[grep(".c17", names(siz))])
+names(siz)[grep(".c18", names(siz))] <- gsub(".c18", ".cX18", names(siz)[grep(".c18", names(siz))])
+names(siz)[grep(".c19", names(siz))] <- gsub(".c19", ".cX19", names(siz)[grep(".c19", names(siz))])
+names(siz)[grep(".c20", names(siz))] <- gsub(".c20", ".cX20", names(siz)[grep(".c20", names(siz))])
+names(siz)[grep(".c30", names(siz))] <- gsub(".c30", ".cX30", names(siz)[grep(".c30", names(siz))])
+
+names(siz)[grep(".c20", names(siz))]
+names(siz)[grep(".c30", names(siz))]
+
+# summarize by site and year ----------------------------------
+
+i <- 1
+# for (j in sitelis) { 
   for (y in yrs)     {  
-    f1 <- fish[which(fish$SampleYear == y & fish$Location == j), ]
-    su <- colSums(f1[11:21])
-    su2 <- su[1] + su[2]
+    f1 <- siz[which(siz$YEAR == y), ] # & siz$SITE.NAME == j), ]
+    if (nrow(f1) > 0) { 
+    mat[i, 3]  <- sum(f1[, grep(".c1", names(siz))], na.rm = T)
+    mat[i, 4]  <- sum(f1[, grep(".c2", names(siz))], na.rm = T)
+    mat[i, 5]  <- sum(f1[, grep(".c3", names(siz))], na.rm = T)
+    mat[i, 6]  <- sum(f1[, grep(".c4", names(siz))], na.rm = T)
+    mat[i, 7]  <- sum(f1[, grep(".c5", names(siz))], na.rm = T)
+    mat[i, 8]  <- sum(f1[, grep(".c6", names(siz))], na.rm = T)
+    mat[i, 9]  <- sum(f1[, grep(".c7", names(siz))], na.rm = T)
+    mat[i, 10] <- sum(f1[, grep(".c8", names(siz))], na.rm = T)
+    mat[i, 11] <- sum(f1[, grep(".c9", names(siz))], na.rm = T)
+    mat[i, 12] <- sum(f1[, grep(".cX10", names(siz))], na.rm = T)
+    mat[i, 13] <- sum(f1[, grep(".cX11", names(siz))], na.rm = T)
+    mat[i, 14] <- sum(f1[, grep(".cX12", names(siz))], na.rm = T)
+    }
     i <- i + 1
-    mat[i, 3:12] <- c(su2, su[3:length(su)])
-    mat[i, 13] <- sum(f1$Total)  }
+  #  }
 } 
-mat
-mat <- mat[-which(rowSums(mat[3:ncol(mat)]) == 0), ]
-mat
+
+head(mat)
+dim(mat)
+#mat <- mat[-which(rowSums(mat[3:ncol(mat)]) == 0), ]
+#mat <- mat[-which(is.na(rowSums(mat[3:ncol(mat)]))), ]
+dim(mat)
+which(rowSums(mat[3:ncol(mat)]) == 0)
+
+# mat is new data object with one site per year ------------
 
 table(mat$site, mat$yr)
 
-#for (i in 1:nrow(mat)) {  mat[i, 3:ncol(mat)] <- mat[i, 3:ncol(mat)] / sum(mat[i, 3:ncol(mat)]/ 100)  }
+plot((colSums(mat[3:14])), type = "h", main = "count data")
+plot(log(colSums(mat[3:14])), type = "h", main = "count data")
+
+barplot(as.matrix(log(mat[3:14] + 1)), beside = T, col = rainbow(length(yrs)), names.arg = names(mat)[3:14])
+
+#for (i in 1:nrow(mat)) {       # converting to proportions results in exactly the same estimates
+#  mat[i, 3:(ncol(mat)-1)] <- mat[i, 3:(ncol(mat)-1)] / sum(mat[i, 3:(ncol(mat)-1)]) * 100  }
 
-barplot(as.matrix(log(mat[3:12] + 1)), beside = T, col = rainbow(length(yrs)), names.arg = names(mat)[3:12])
+# look at size spectra and subset classes that are selected for ------------------
 
-plot((colSums(mat[3:12])), type = "h")
-plot(log(colSums(mat[3:12])), type = "h")
+classes <- 3:10  
 
-if (grp == "com")  {  classes <- 2:5  }
-if (grp != "com")  {  classes <- 1:5  }
+head(mat[, classes + 2], 30)
 
-head(mat[, classes + 2], 100)
+# calculate slope of the size spectrum for each site and year --------------------
 
 mat$slope <- NA
 mat$SE <- NA
 
-par(mfrow = c(8, 5))
-
+par(mfrow = c(4, 5), mex = 0.5)
 for (i in 1:nrow(mat)) { 
   ve <- as.numeric(log(mat[i, classes + 2]))
   plot(ve, type = "h", axes = F)
-  axis(1, at = 1:10, lab = names(mat)[3:12], las = 2)
-  mtext(side = 3, text = paste(mat$site[i], mat$yr[i]))
+  axis(1, at = 1:8, lab = names(mat[i, classes + 2]), las = 2)
+  mtext(side = 3, text = paste(mat$yr[i]))
 
   ve[which(ve == "-Inf")] <- NA
   out <- lm(ve ~ classes)
@@ -267,18 +349,15 @@ for (i in 1:nrow(mat)) {
     mat$SE[i] <- summary(out)$coef[2, 2]   }
 }
 
+# merge slope estimates with data file and view ---------------------
 
-mat <- merge(mat, met, by.x = "site", by.y = "Location")
+#mat <- merge(mat, met, by.x = "site", by.y = "Location")
 mat
 mat$yr <- as.factor(mat$yr)
 
 dev.off()
 
-plot(log(mat$total) ~ mat$site)
-plot(log(mat$total) ~ mat$yr)
-
-boxplot(mat$slope ~ mat$site)
-boxplot(mat$slope ~ mat$yr)
+plot(mat$yr, mat$slope)
 
 metric <- mat$slope
 

indicator_processing/non_automated/TRCRMP_benthic.R → indicator_processing/non_automated/TCRMP_benthic.R

@@ -45,6 +45,8 @@ cds$Meaning[grep("spp.", cds$Meaning)]
 cds$Cat2[grep("spp.", cds$Meaning)] <- "X"
 
 # create list of coral species ---------------------------
+
+cds$Meaning
 cds$Meaning[which(cds$Cat2 == "Coral")]
 coralspplis <- cds$Code[which(cds$Cat2 == "Coral")]
 coralspplis
@@ -79,6 +81,7 @@ ben <- ben[which(ben$Period == "Annual"), ]
 bensum <- bensum[which(bensum$Period == "Annual"), ]
 
 # create unique sample ID, for merging databases ------------------------
+
 ben$ID <- paste0(ben$SampleYear, ben$SampleMonth, ben$Location, ben$Transect)
 bensum$ID <- paste0(bensum$SampleYear, bensum$SampleMonth, bensum$Location, bensum$Transect)
 length(unique(ben$ID))
@@ -155,7 +158,7 @@ summary(out1)$coef[, 1]
 mod <- summary(out1)$coef[, 1]
 modse <- summary(out1)$coef[, 2]
 
-plot(2001:2021, ind, ylim = c(min(ind) - 2, max(ind + 2)))
+plot(2001:2021, ind, ylim = c(min(ind) - 2, max(ind + 2)), main = varint)
 lines(2001:2021, ind + indse)
 lines(2001:2021, ind - indse)
 points(2001:2021, mod, col = 2)
@@ -164,8 +167,8 @@ lines(2001:2021, mod - modse, col = 2, lty = 2)
 
 cor(ind, mod)
 
-if (varint == "sprich")  {  save(out1, file = "C:/Users/mandy.karnauskas/Desktop/Caribbean-ESR/indicator_data/PRCRMP/coralspprich_USVI.RData")  }
-if (varint == "percov")  {  save(out1, file = "C:/Users/mandy.karnauskas/Desktop/Caribbean-ESR/indicator_data/PRCRMP/percoralcov_USVI.RData")   }
+if (varint == "sprich")  {  save(out1, file = "C:/Users/mandy.karnauskas/Desktop/Caribbean-ESR/indicator_data/TCRMP/coralspprich_USVI.RData")  }
+if (varint == "percov")  {  save(out1, file = "C:/Users/mandy.karnauskas/Desktop/Caribbean-ESR/indicator_data/TCRMP/percoralcov_USVI.RData")   }
 }