NB - Må kjøres på NINA sin Rstudio server
# Connect to server and get trim results data (MUST BE RUN IN RSTUDIO ON NINA SERVER)
# myconn <- DBI::dbConnect(odbc::odbc(),
# Driver = "FreeTDS",
# Server = "ninsql07.nina.no",
# Database = "TOVTaksering",
# UID = "TOVeLes",
# PWD = "gLuteusMax1mus",
# Port = 1433)
#
# # Load table TrimResults
# TrimResults = as_tibble(dbGetQuery(myconn, paste("SELECT * FROM TrimResults")))
#
# Species = as_tibble(dbGetQuery(myconn, paste("SELECT * FROM Art")))
#
# dbDisconnect(myconn)
#
# write.csv(Species, file = paste(getwd(), "/Species.csv", sep =""))
# write.csv(TrimResults, file = paste(getwd(), "/TrimResults.csv", sep =""))# Species table
Species = read.csv("P:/41201042_okologisk_tilstand_fastlandsnorge_2020_dataanaly/fjell2021/data/TRIM Fjellrype/Species.csv")
Species = Species %>% filter(ArtsID == 1540)
# Observation data
TrimResults = read.csv("P:/41201042_okologisk_tilstand_fastlandsnorge_2020_dataanaly/fjell2021/data/TRIM Fjellrype/TrimResults.csv")
TrimResults = TrimResults %>% filter(Species == 1540, Year >= 2010)
TrimResults = TrimResults %>% mutate(Count = ifelse(Count == -1, NA, Count))
TrimResults = TrimResults %>% dplyr::select(-X)
# Get data linking sites and regions
trimdata_orig_pass = read.csv("P:/41201042_okologisk_tilstand_fastlandsnorge_2020_dataanaly/Fuglindeks/data/ObsDat/ObsDat_Trim_pass_2020-09-29.csv")
trimdata_orig_line = read.csv("P:/41201042_okologisk_tilstand_fastlandsnorge_2020_dataanaly/Fuglindeks/data/ObsDat/ObsDat_Trim_line_2020-09-29.csv")
trimdata_orig = rbind(trimdata_orig_pass, trimdata_orig_line)
trimdata_orig = trimdata_orig %>% dplyr::select(site, region)
trimdata_orig = trimdata_orig %>% mutate(region = ifelse(region == "øst", "øst", ifelse(region == "sør", "sør", region))) %>% rename(Site = site)
trimdata_orig = trimdata_orig[!duplicated(trimdata_orig$Site), ]
fjellrype_info = Species %>% dplyr::select(ArtsID, Artsnavn_Norsk) %>% rename(Species = ArtsID)
fjellrype = left_join(TrimResults, fjellrype_info, by = "Species")
fjellrype = left_join(fjellrype, trimdata_orig, by = c("Site"))# Fjellrype - 1540
fjellNord = fjellrype %>% filter(region == "nord")
fjellMidt = fjellrype %>% filter(region == "midt")
fjellOst = fjellrype %>% filter(region == "øst")
fjellVest = fjellrype %>% filter(region == "vest")
fjellSor = fjellrype %>% filter(region == "sør")
# Run trim
fjellNordMod = trim(Count ~ Site + Year, data = fjellNord, model = 2, changepoints = "all", stepwise = FALSE, serialcor = TRUE, overdisp = TRUE)
fjellMidtMod = trim(Count ~ Site + Year, data = fjellMidt, model = 2, changepoints = "all", stepwise = FALSE, serialcor = TRUE, overdisp = TRUE)
fjellOstMod = trim(Count ~ Site + Year, data = fjellOst, model = 2, changepoints = "all", stepwise = FALSE, serialcor = TRUE, overdisp = TRUE)
fjellVestMod = trim(Count ~ Site + Year, data = fjellVest, model = 2, changepoints = "all", stepwise = FALSE, serialcor = TRUE, overdisp = TRUE)
fjellSorMod = trim(Count ~ Site + Year, data = fjellSor, model = 2, changepoints = "all", stepwise = FALSE, serialcor = TRUE, overdisp = TRUE)# Create the index table
fjellTab = data.frame(matrix(data = c("fjellrype", 1540), nrow = 1, ncol = 2, dimnames = list(c(), c("Species", "SpeciesID"))))
fjellTab = fjellTab[rep(1:nrow(fjellTab), rep(c(5), dim(fjellTab)[1])), ]
fjellTab1 = fjellTab %>% mutate(Rectype = 1, "2010" = 0, "2011" = 0, "2012" = 0, "2013" = 0, "2014" = 0, "2015" = 0, "2016" = 0, "2017" = 0, "2018" = 0, "2019" = 0, "2020" = 0, Region = rep(c("nord", "midt", "øst", "vest", "sør"), 1), Slope_add = 0, Slope_add_se = 0, Slope_mul = 0, Slope_mul_se = 0, p_trend = 0, Trend_class = 0)
fjellTab2 = fjellTab1 %>% mutate(Rectype = 2)
fjellTab3 = fjellTab1 %>% mutate(Rectype = 3)
fjellTab4 = fjellTab1 %>% mutate(Rectype = 4)
rypeNordIndex = index(fjellNordMod, which = "imputed", base = 1)
rypeNordTot = totals(fjellNordMod, which = "imputed", obs = TRUE)
rypeMidtIndex = index(fjellMidtMod, which = "imputed", base = 1)
rypeMidtTot = totals(fjellMidtMod, which = "imputed", obs = TRUE)
rypeOstIndex = index(fjellOstMod, which = "imputed", base = 1)
rypeOstTot = totals(fjellOstMod, which = "imputed", obs = TRUE)
rypeVestIndex = index(fjellVestMod, which = "imputed", base = 1)
rypeVestTot = totals(fjellVestMod, which = "imputed", obs = TRUE)
rypeSorIndex = index(fjellSorMod, which = "imputed", base = 1)
rypeSorTot = totals(fjellSorMod, which = "imputed", obs = TRUE)
fjellTab1[1, 4:14] = rypeNordIndex$imputed
fjellTab1[2, 4:14] = rypeMidtIndex$imputed
fjellTab1[3, 4:14] = rypeOstIndex$imputed
fjellTab1[4, 4:14] = rypeVestIndex$imputed
fjellTab1[5, 4:14] = rypeSorIndex$imputed
fjellTab2[1, 4:14] = rypeNordIndex$se_imp
fjellTab2[2, 4:14] = rypeMidtIndex$se_imp
fjellTab2[3, 4:14] = rypeOstIndex$se_imp
fjellTab2[4, 4:14] = rypeVestIndex$se_imp
fjellTab2[5, 4:14] = rypeSorIndex$se_imp
fjellTab3[1, 4:14] = rypeNordTot$imputed
fjellTab3[2, 4:14] = rypeMidtTot$imputed
fjellTab3[3, 4:14] = rypeOstTot$imputed
fjellTab3[4, 4:14] = rypeVestTot$imputed
fjellTab3[5, 4:14] = rypeSorTot$imputed
fjellTab4[1, 4:14] = rypeNordTot$se_imp
fjellTab4[2, 4:14] = rypeNordTot$se_imp
fjellTab4[3, 4:14] = rypeNordTot$se_imp
fjellTab4[4, 4:14] = rypeNordTot$se_imp
fjellTab4[5, 4:14] = rypeNordTot$se_imp# For loop that gets the additive and multiplicative slopes, as well as the trend class for each model
model_list = list(fjellNordMod, fjellMidtMod, fjellOstMod, fjellVestMod, fjellSorMod)
tmp_add = c()
tmp_add_se = c()
tmp_mul = c()
tmp_mul_se = c()
tmp_p = c()
tmp_meaning = c()
for(i in 1:length(model_list)){
tmp_add = append(tmp_add, overall(model_list[[i]], which = "imputed")$slope$add)
tmp_add_se = append(tmp_add_se, overall(model_list[[i]], which = "imputed")$slope$se_add)
tmp_mul = append(tmp_mul, overall(model_list[[i]], which = "imputed")$slope$mul)
tmp_mul_se = append(tmp_mul_se, overall(model_list[[i]], which = "imputed")$slope$se_mul)
tmp_p = append(tmp_p, overall(model_list[[i]], which = "imputed")$slope$p)
tmp_meaning = append(tmp_meaning, overall(model_list[[i]], which = "imputed")$slope$meaning)
}
fjellTab1$Slope_add = tmp_add
fjellTab1$Slope_add_se = tmp_add_se
fjellTab1$Slope_mul = tmp_mul
fjellTab1$Slope_mul_se = tmp_mul_se
fjellTab1$p_trend = tmp_p
fjellTab1$Trend_class = tmp_meaning
fjellTab_df = rbind(fjellTab1, fjellTab2, fjellTab3, fjellTab4)
fjellTab_df = fjellTab_df %>% arrange(factor(Region, levels = c("nord", "midt", "øst", "vest", "sør")))
fjellTab_df[, 4:14] = round(fjellTab_df[, 4:14], digits = 3)
fjellTab_df[, 16:20] = round(fjellTab_df[, 16:20], digits = 3)
fjellTab_df = fjellTab_df %>% mutate(Slope_add = ifelse(Rectype != 1, NA, Slope_add), Slope_add_se = ifelse(Rectype != 1, NA, Slope_add_se), Slope_mul = ifelse(Rectype != 1, NA, Slope_mul), Slope_mul_se = ifelse(Rectype != 1, NA, Slope_mul_se), p_trend = ifelse(Rectype != 1, NA, p_trend), Trend_class = ifelse(Rectype != 1, NA, Trend_class))
fjellTab_df[1:8,]## Species SpeciesID Rectype 2010 2011 2012 2013 2014 2015
## 1 fjellrype 1540 1 1 1.951 1.357 1.709 2.376 2.212
## 11 fjellrype 1540 2 0 0.994 0.664 0.842 1.133 1.010
## 12 fjellrype 1540 3 58 113.000 79.000 99.000 138.000 129.000
## 13 fjellrype 1540 4 26 29.000 18.000 23.000 27.000 16.000
## 1.1 fjellrype 1540 1 1 1.621 2.317 0.750 1.798 1.782
## 1.11 fjellrype 1540 2 0 1.078 1.462 0.528 1.139 1.144
## 1.12 fjellrype 1540 3 17 27.000 38.000 12.000 30.000 29.000
## 1.13 fjellrype 1540 4 26 29.000 18.000 23.000 27.000 16.000
## 2016 2017 2018 2019 2020 Region Slope_add Slope_add_se
## 1 2.317 2.478 1.565 1.716 1.652 nord 0.029 0.026
## 11 1.049 1.126 0.719 0.790 0.784 nord NA NA
## 12 135.000 144.000 91.000 100.000 96.000 nord NA NA
## 13 16.000 18.000 12.000 14.000 17.000 nord NA NA
## 1.1 3.263 2.543 3.034 1.993 2.446 midt 0.083 0.034
## 1.11 2.006 1.580 1.866 1.286 1.519 midt NA NA
## 1.12 54.000 42.000 50.000 33.000 40.000 midt NA NA
## 1.13 16.000 18.000 12.000 14.000 17.000 midt NA NA
## Slope_mul Slope_mul_se p_trend Trend_class
## 1 1.029 0.027 0.298 Uncertain
## 11 NA NA NA <NA>
## 12 NA NA NA <NA>
## 13 NA NA NA <NA>
## 1.1 1.087 0.037 0.036 Moderate increase (p<0.05)
## 1.11 NA NA NA <NA>
## 1.12 NA NA NA <NA>
## 1.13 NA NA NA <NA># Export the index and totals table
# write_xlsx(fjellTab_df, path = "P:/41201042_okologisk_tilstand_fastlandsnorge_2020_dataanaly/fjell2021/data/TRIM Fjellrype/fjellrypeIndex.xlsx")Lag en bootstrapping for-løkke som sampler lognormale indeksverdier og standardavvik for å finne gjennomsnittlig indeks for hvert år.
### Take the mean of 10000 samples for fjellrype in each region
lnorm_fjellNord = normal2Lognormal(fjellTab1[1, 4:14], fjellTab2[1, 4:14])
lnorm_fjellMidt = normal2Lognormal(fjellTab1[2, 4:14], fjellTab2[2, 4:14])
lnorm_fjellOst = normal2Lognormal(fjellTab1[3, 4:14], fjellTab2[3, 4:14])
lnorm_fjellVest = normal2Lognormal(fjellTab1[4, 4:14], fjellTab2[4, 4:14])
lnorm_fjellSor = normal2Lognormal(fjellTab1[5, 4:14], fjellTab2[5, 4:14])
fjellNordSmp = c()
fjellMidtSmp = c()
fjellOstSmp = c()
fjellVestSmp = c()
fjellSorSmp = c()
tmp1 = c()
tmp2 = c()
tmp3 = c()
tmp4 = c()
tmp5 = c()
k = 10000
for(i in 1:length(unique(fjellNord$Year))){
for(j in 1:k){
tmp1 = append(tmp1, rlnorm(n = 1, meanlog = lnorm_fjellNord$mean[[i]],
sdlog = lnorm_fjellNord$sd[[i]]))
tmp2 = append(tmp2, rlnorm(n = 1, meanlog = lnorm_fjellMidt$mean[[i]],
sdlog = lnorm_fjellMidt$sd[[i]]))
tmp3 = append(tmp3, rlnorm(n = 1, meanlog = lnorm_fjellOst$mean[[i]],
sdlog = lnorm_fjellOst$sd[[i]]))
tmp4 = append(tmp4, rlnorm(n = 1, meanlog = lnorm_fjellVest$mean[[i]],
sdlog = lnorm_fjellVest$sd[[i]]))
tmp5 = append(tmp5, rlnorm(n = 1, meanlog = lnorm_fjellSor$mean[[i]],
sdlog = lnorm_fjellSor$sd[[i]]))
fjellNordSmp = append(fjellNordSmp, mean(tmp1))
fjellMidtSmp = append(fjellMidtSmp, mean(tmp2))
fjellOstSmp = append(fjellOstSmp, mean(tmp3))
fjellVestSmp = append(fjellVestSmp, mean(tmp4))
fjellSorSmp = append(fjellSorSmp, mean(tmp5))
tmp1 = c()
tmp2 = c()
tmp3 = c()
tmp4 = c()
tmp5 = c()
}
}
fjellNordSmp = matrix(data = fjellNordSmp, nrow = k, ncol = length(unique(fjellNord$Year)), dimnames = list(c(), c(as.character(2010:2020))))
fjellMidtSmp = matrix(data = fjellMidtSmp, nrow = k, ncol = length(unique(fjellMidt$Year)), dimnames = list(c(), c(as.character(2010:2020))))
fjellOstSmp = matrix(data = fjellOstSmp, nrow = k, ncol = length(unique(fjellOst$Year)), dimnames = list(c(), c(as.character(2010:2020))))
fjellVestSmp = matrix(data = fjellVestSmp, nrow = k, ncol = length(unique(fjellVest$Year)), dimnames = list(c(), c(as.character(2010:2020))))
fjellSorSmp = matrix(data = fjellSorSmp, nrow = k, ncol = length(unique(fjellSor$Year)), dimnames = list(c(), c(as.character(2010:2020))))
# Mean of the 10000 samples
fjellSmpMean = rbind(colMeans(fjellNordSmp), colMeans(fjellMidtSmp), colMeans(fjellOstSmp), colMeans(fjellVestSmp), colMeans(fjellSorSmp))
fjellSmpMean = as_tibble(fjellSmpMean) %>% mutate(Region = c("nord", "midt", "øst", "vest", "sør"))
head(fjellSmpMean)## # A tibble: 5 x 12
## `2010` `2011` `2012` `2013` `2014` `2015` `2016` `2017` `2018` `2019` `2020`
## <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 1 1.94 1.35 1.72 2.37 2.22 2.33 2.47 1.56 1.71 1.64
## 2 1 1.61 2.31 0.749 1.80 1.79 3.27 2.54 3.04 2.00 2.44
## 3 1 3.30 1.87 1.15 6.78 3.13 8.38 3.59 2.91 4.02 4.36
## 4 1 0.283 1.44 0.789 1.03 1.38 0.752 1.30 1.44 1.08 1.29
## 5 1 1.17 1.30 0.515 0.504 1.85 0.279 0.534 0.720 0.455 2.41
## # ... with 1 more variable: Region <chr>fjellSmp = rbind(fjellNordSmp, fjellMidtSmp, fjellOstSmp, fjellVestSmp, fjellSorSmp)
fjellSmp = as_tibble(fjellSmp) %>% mutate(Region = c(rep("nord", 10000), rep("midt", 10000), rep("øst", 10000), rep("vest", 10000), rep("sør", 10000)))
head(fjellSmp)## # A tibble: 6 x 12
## `2010` `2011` `2012` `2013` `2014` `2015` `2016` `2017` `2018` `2019` `2020`
## <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 1 0.712 1.72 1.26 4.51 1.20 2.74 0.930 0.863 1.05 1.47
## 2 1 2.08 1.18 2.15 2.69 1.74 3.59 2.46 1.68 0.929 2.45
## 3 1 2.42 1.29 1.90 1.48 2.02 1.40 1.19 2.06 1.01 1.32
## 4 1 2.94 0.868 1.60 1.63 2.95 1.80 1.48 1.41 3.01 0.990
## 5 1 0.613 1.62 2.28 2.99 0.982 1.34 2.22 2.82 1.45 2.00
## 6 1 1.05 1.05 2.52 2.20 2.40 1.11 2.65 0.696 1.46 2.67
## # ... with 1 more variable: Region <chr># write_xlsx(fjellSmpMean, path = "P:/41201042_okologisk_tilstand_fastlandsnorge_2020_dataanaly/fjell2021/data/TRIM Fjellrype/fjellrypeSmpMean.xlsx")
# write_xlsx(fjellSmp, path = "P:/41201042_okologisk_tilstand_fastlandsnorge_2020_dataanaly/fjell2021/data/TRIM Fjellrype/fjellrypeSmp.xlsx")