Rein v3 - Usikkerhet i bestandstall
Start
Dette er en felles indikator for villrein og tamreinsbestander Utarbeidelsen av metodikken er gjort av Erik Fremstad. I denne varsjonen av indikatoren legger vi til en usikkerhet i bestandsestimatet - litt større for villrein enn for tamrein. Det er usikkert hvor presist dette er, men man kan gjerne tolke usikkerhet til også å omfavne usikkerheten i referanseverdien
Import
Datasett med bestandstall, areal, og TRIverdier. Topographic roughenss index brukes til å modellere referansetilstand. For hvert reinområde andgis areal per region.
dat <- read_excel("../data/rein.xlsx")
head(dat)## # A tibble: 6 × 13
## type navn navn2 Totalreal fjellareal bestand `TRI-fjell` E S W
## <chr> <chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 Villr… Sete… Sete… 7067. 5119. 3500 364. NA 3212. 1907.
## 2 Villr… Sete… Sete… 2433. 660. 2000 437. NA 660. NA
## 3 Villr… Skau… Skau… 654. 604. 60 573. NA NA 604.
## 4 Villr… Våmu… Våmu… 449. 51.2 240 382. NA 51.2 NA
## 5 Villr… Brat… Brat… 417. 213. 500 326. NA 213. NA
## 6 Villr… Blef… Blef… 264. 64.4 140 324. 29.2 35.2 NA
## # … with 3 more variables: C <dbl>, N <dbl>, totalt <dbl>Kolonnen med total areal (helt til høyre) stemmer nesten med total fjellareal, så jeg kan bruke de vekselsvis.
Rydder
dat <- dplyr::select(dat, -navn2)Så setter jeg alle NA til 0. Dette inkluderer totalbestand for Trollheimen. Trollheimen har tom celle for totalbestand. Siden det er skrevet 0 i cellen for fjellbestand så antar jeg det skal være null i den andre også.
dat[is.na(dat)] <- 0Fordeling av areal mellom regioner.
dat$propN <- dat$N/dat$totalt
dat$propC <- dat$C/dat$totalt
dat$propE <- dat$E/dat$totalt
dat$propW <- dat$W/dat$totalt
dat$propS <- dat$S/dat$totaltMelt
Nå deler jeg opp villreinområdene i så mange regioner som de inngår i
setDT(dat)
dat2 <- melt(dat,
measure.vars = c("N", "C", "E", "W", "S"),
value.name = "area",
variable.name= "reg",
na.rm=T)Dette gir fem rader per område. Så tar jeg bort de radene der det ikke finnes noe (fjell)areal
temp <- dat2[dat2$type=="ingen rein",]
dat2 <- dat2[dat2$area>0,]
dat2 <- rbind(dat2, temp)
rm(temp)Det er noen rader med type = ingen rein som vi må kombinere på et vis.
temp <- dat2[dat2$type=="ingen rein",]
dat2 <- dat2[dat2$type!="ingen rein",]temp2 <- temp[temp$navn == "Fjell uten rein - nord" &
temp$reg %in% c("N"),]
temp3 <- temp[temp$navn == "Fjell uten rein - sør" &
temp$reg %in% c("E", "S", "W", "C"),]
temp4 <- temp[temp$navn =="Fjellområder uten rein",]
temp4 <- temp4[!duplicated(temp4$reg)]
temp <- rbind(temp2, temp3)
temp$area <- temp4$area[match(temp$reg, temp4$reg)]dat2 <- rbind(dat2, temp)Så kan jeg regne ut referansetettheter.
dat2$ref <- 0.75*1.0759*exp(-0.001*dat2$`TRI-fjell`)Så legger jeg på en usikkerhet rundt betsnadsestimatene. 10%CV for villrein og 5%CV for tamrein.
Her er de observerte dyretetthetene. Siden tettheten er antatt lik i fjell og ikke-fjell trenger vi ikke gå veien innom å regne ut bestandstall
dat2$tetthet <- dat2$bestand/dat2$Totalreal
#Endrer NaN til 0
dat2$tetthet[dat2$type=="ingen rein"] <- 0
dat2$ID <- paste0(dat2$navn, dat2$reg)Jeg legger på en CV på 10% for villrein og 5% for tamrein. Dette gir en spredning i verdier som ser slik ut:
temp <- seq(0,150, 1)
CV10 <- rnorm(150, 1, .1)
ggplot()+
geom_point(mapping = aes(y = temp*CV10,
x = 1:length(temp)))
nsim <- 1000
refs <- c()
IDs <- c()
reg <- c()
area <- c()
ind <- c()
indtemp <- c()
for(i in dat2$ID){
temp <- dat2[dat2$ID == i,]
ifelse(temp$type == "Tamrein",
temp2 <- rnorm(nsim, 1, 0.05),
temp2 <- rnorm(nsim, 1, 0.1))
indtemp <- (temp$bestand/temp$Totalreal)*temp2
IDs <- c(IDs, rep(i, nsim))
refs <- c(refs, rep(temp$ref, nsim))
reg <- c(reg, rep(as.character(temp$reg), nsim))
area <- c(area, rep(temp$area, nsim))
ind <- c(ind, indtemp)
}sims <- data.table(
reg = reg,
ID = IDs,
area = area,
refs = refs,
ind = ind)Det blir NaN der vi har null rein.
unique(sims$ID[sims$ind== "NaN"])## [1] "Fjell uten rein - nordN" "Fjell uten rein - sørC"
## [3] "Fjell uten rein - sørE" "Fjell uten rein - sørW"
## [5] "Fjell uten rein - sørS"Disse kan vi sette til null.
sims$ind[sims$ind == "NaN"] <- 0ggplot(data = sims, aes(x = reg, y = ind))+
geom_violin()
ggplot(data = sims, aes(x = ID, y = ind))+
geom_violin(trim=F)+
coord_flip()
summary(sims$ind)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.0000 0.2897 0.6078 0.8336 1.1702 3.7272Det ser bra ut. Variasjonen øker når gjennomsnittet øker.
sims$diff <- sims$ref-sims$ind
dotchart(sample(sims$diff, 1000, replace = F))
Her må vi ha en tosidig skalering. Erik har laget et forslag til skalering som ser slik ut
temp <- data.frame(x = c(0, 0.6, 1, 2,10),
y = c(0, 0.6, 1, 0.6, 0))
(mod <- ggplot(data = temp, aes(x = x, y = y))+
geom_line(size=1.3, linetype="dashed", colour="grey30")+
geom_point(size=5)+
theme_bw()+
ylab("Skalert verdi")+
xlab("Obs./Ref.")
)
sims$rat <- sims$ind/sims$ref
dotchart(sample(sims$rat, 1000, replace = F))
rm(dat2, dat)sims$val <- NA
sims$val[sims$rat<1] <- sims$rat[sims$rat<1]
sims$val[sims$rat %between% c(1, 2)] <- 1-(0.4*((sims$rat[sims$rat %between% c(1, 2)])-1))
sims$val[sims$rat >2] <- (1-((sims$rat[sims$rat >2]-2)/(10-2)))/1.6
(mod + geom_point(data=sims, aes(x = rat, y = val),
size=3, alpha=.005, stroke=3))
summary(sims$val)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.0000 0.4578 0.6672 0.6425 0.9218 1.0000ggplot(data = sims, aes(x = ID, y = val))+
geom_violin(trim=F)+
coord_flip()
ggplot(data = sims, aes(x = reg, y = val))+
geom_point(position = position_dodge2(0.4),
aes(size = area),
alpha = 0.01)+
scale_size_continuous(range = c(1, 10))
Fint. Så trenger jeg et arealvektet gjennomsnitt per region. Bruker 100*500 resamplings.
nsim <- 500
vals <- c()
reg <- c()
temp2 <- c()
for(i in unique(sims$reg)){
temp <- sims[sims$reg==i,]
for(i in 1:nsim){
temp3 <- c()
wgt <- c()
for(m in 1:100){
temp2 <- temp[runif(1,0,nrow(temp)),]
temp3 <- c(temp3, temp2$val)
wgt <- c(wgt, temp2$area)
}
vals <- c(vals, weighted.mean(temp3, wgt))
reg <- c(reg, temp$reg[1])
}
}
df <- data.frame(reg = reg,
val = vals,
year = 2020,
X = "rein")ggplot(data = df, aes(x = reg, y = val))+
geom_violin(trim=F, fill="grey")+
theme_bw()
Så må vi ha en verdi for hele Norge Fjellareal per region:
wgt <- readRDS("../data/fjellareal.rds")
wgt$Fjellareal2 <- wgt$Fjellareal/max(wgt$Fjellareal)
wgt$reg <- c("N", "C", "E", "W", "S")regions <- c("N", "C", "W", "S", "E")
temp <- c()
temp2 <- c()
for(i in 1:1000){
temp <- c(
sample(df$val[df$reg == "N"], 1),
sample(df$val[df$reg == "C"], 1),
sample(df$val[df$reg == "E"], 1),
sample(df$val[df$reg == "W"], 1),
sample(df$val[df$reg == "S"], 1)
)
temp2 <- c(temp2, weighted.mean(temp, wgt$Fjellareal2))
}
Norge <- data.frame(
reg = rep("Norge", 1000),
val = temp2,
year = 2020,
X = "rein"
)
df <- rbind(df, Norge)df$reg <- as.factor(df$reg)Plotting
temp2 <- c("N", "C", "E", "S", "W", "Norge")
temp3 <- c("Nord-Norge", "Midt-Norge", "Østlandet", "Sørlandet", "Vestlandet", "Norge")df2 <- df
df <- aggregate(data = df,
val~reg,
FUN = function(x) c(quantile(x, c(0.025, .5, 0.975))))
df <- do.call(data.frame, df)
names(df) <- c("reg", "low", "med", "upp")(gg <- ggplot(data = df, aes(x = factor(reg, levels = temp2), y = med))+
geom_bar(stat="identity", colour="black", size=1.2,
fill="grey80")+
geom_errorbar(aes(ymax = upp, ymin=low), width=.5, size=2)+
theme_bw(base_size = 20)+
ylab("Tetthet av rein\nskalert mot referansetilstand")+
xlab("")+
theme(
axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)
)+
scale_x_discrete(breaks = temp2,
labels = temp3)+
scale_y_continuous(breaks = seq(0,1,.2),
limits = c(0,1))
) 
png("../output/indicatorPlots/skalert/rein.png",
units="in", width=6, height=8, res=300)
gg
dev.off()write.csv(df2, "../output/indicator_values/rein.csv", row.names = F)