Data jsou stažena z https://www.movebank.org/. Data pochází ze GPS sledování samice reintrodukovaného geparda žijícího v Jižní Africe v národním parku Pilanesberg. Data se začali sbírat v roce 2014, kdy samice neměla žádné mláďata a ani nebyla březí (v datech jako cheetah_1). V listopadu toho stejného roku zabřezla (Cheetah_2) a v roce 2015 porodila 4 mláďata samce a jedno za chvíli umřelo (Cheetah_3).
library(moveHMM)
Vstupní data musí mít správný formát, aby bylo možné data zpracovat a poté je analyzovat. Data-frame musí obsahovat 2 sloupce:
sloupec s názvem “x”: Easting nebo zeměpisná délka (longitude)
sloupec s názvem “y”: Northing nebo zeměpisná šířka (latitude)
Důležitým sloupec je “ID”. Pokud není ve vstupních datech žádný sloupec s názvem “ID”, tak balíček bude předpokládat, že daná data jsou pouze pro jednoho jedince.
gepard <- read.table("Gepard_1.csv", header = T, sep = ";", dec=",")
gepard$time = as.POSIXct(gepard$time)
gepard_D <- read.table("Gepard_2.csv", header = T, sep = ";", dec=",")
gepard_D$time = as.POSIXct(gepard_D$time)
gepard_PD <- read.table("Gepard_3.csv", header = T, sep = ";", dec=",")
gepard_PD$time = as.POSIXct(gepard_PD$time)
Převod do souřadnicového systému UTM:
library(sf)
floor((180+gepard$long[1])/6) + 1
gepard_sf <- st_as_sf(gepard, coords = c("long", "lat"), crs = 4326)
gepard_utm <- st_transform(gepard_sf, crs = 32635)
gepard <- cbind(gepard_utm, st_coordinates(gepard_utm))
gepard_D_sf <- st_as_sf(gepard_D, coords = c("long", "lat"), crs = 4326)
gepard_D_utm <- st_transform(gepard_D_sf, crs = 32635)
gepard_D <- cbind(gepard_D_utm, st_coordinates(gepard_D_utm))
gepard_PD_sf <- st_as_sf(gepard_PD, coords = c("long", "lat"), crs = 4326)
gepard_PD_utm <- st_transform(gepard_PD_sf, crs = 32635)
gepard_PD <- cbind(gepard_PD_utm, st_coordinates(gepard_PD_utm))
Náhled na vstupní data:
## Simple feature collection with 6 features and 6 fields
## geometry type: POINT
## dimension: XY
## bbox: xmin: 507733.9 ymin: -2793073 xmax: 509791.7 ymax: -2792354
## projected CRS: WGS 84 / UTM zone 35N
## ID TEMP Height time X Y
## 1 Cheetah_1 13.5 1161 2014-04-12 04:00:33 509791.7 -2792354
## 2 Cheetah_1 25.5 1164 2014-04-12 22:00:33 507754.0 -2793037
## 3 Cheetah_1 22.5 1167 2014-04-13 04:00:33 507754.0 -2793037
## 4 Cheetah_1 30.5 1167 2014-04-13 10:00:33 507937.4 -2792934
## 5 Cheetah_1 28.0 1161 2014-04-13 16:00:33 507752.0 -2793039
## 6 Cheetah_1 22.0 1164 2014-04-13 22:00:33 507733.9 -2793073
## geometry
## 1 POINT (509791.7 -2792354)
## 2 POINT (507754 -2793037)
## 3 POINT (507754 -2793037)
## 4 POINT (507937.4 -2792934)
## 5 POINT (507752 -2793039)
## 6 POINT (507733.9 -2793073)
Data obsahují jak sloupce se souřadnicemi, tak i sloupce s možnými kovariátami:
TEMP = venkovní teplota (°C)
Height = nadmořská výška (m n. m.)
Navíc je tady i sloupec ID (označení “jedince”) a čas sběru dat (time)
V datech lze vidět, že některá měření chybí a je potřeba je doplnit, tak aby bylo každé měření po 6 hodinách.
gepard_1a=data.frame(time=seq.POSIXt(from = min(gepard$time), by= "6 hour",length.out = 788))
gepard_1b=data.frame(time=seq.POSIXt(from = gepard$time[788],to=max(gepard$time), by= "6 hour"))
geparad_data1=rbind(gepard_1a,gepard_1b)
gepard = merge(geparad_data1,gepard,by="time",all.x = TRUE)
#gepard_D (březí)
gepard_2a=data.frame(time=seq.POSIXt(from = min(gepard_D$time), to =max(gepard_D$time), by= "6 hour"))
gepard_D = merge(gepard_2a,gepard_D,by="time",all.x = TRUE)
#gepard_PD (s mláďaty)
gepard_3a=data.frame(time=seq.POSIXt(from = min(gepard_PD$time), to =max(gepard_PD$time), by= "6 hour"))
gepard_PD = merge(gepard_3a,gepard_PD,by="time",all.x = TRUE)
Hodnoty “x” a “y” jsou v datech vyjádřeny v metrech, proto je potřeba data přetransformovat, aby byla délka v kilometrech.
gepard$X <- gepard$X/1000
gepard$Y <- gepard$Y/1000
gepard_D$X <- gepard_D$X/1000
gepard_D$Y <- gepard_D$Y/1000
gepard_PD$X <- gepard_PD$X/1000
gepard_PD$Y <- gepard_PD$Y/1000
Před použití následující funkce „PrebData“ je potřeba pro správné fungovaní baličku spustit následující funkci:
get_closest_id <- function(i, ids=gepard$ID){
id <- ids[i]
ii <- i
while (is.na(id) & ii > 1){
ii <- ii - 1
id <- ids[ii]
}
ii <- i
while (is.na(id) & ii < length(ids)){
ii <- ii + 1
id <- ids[ii]
}
return(id)
}
for (i in which(is.na(gepard$ID))) {
gepard$ID[i] <- get_closest_id(i)
}
Funkce “PrebData” umožňuje vymodelovat délky kroků a jejich úhel otočení a je vstupním objektem pro tvorbu modelu.
Duležité je zvolit správné argumenty funkce:
type: určuje, zda-li jsou souřadnice v easting/northing (type=“UTM”), a nebo v zeměpisné šířce/délce (type=“LL”), který je výchozí
coordNames: pouřívá se v případě, když jsou sloupce se souřadnicemi nejsou popsané “x” a “x”.
Funkce “PrebData”:
data <- prepData(gepard,type="UTM",coordNames=c("X","Y"))
## Warning in prepData(gepard, type = "UTM", coordNames = c("X", "Y")): There are
## 36 missing covariate values. Each will be replaced by the closest available
## value.
Varování upozorňuje na data bez kovariáty a říká, že nahradí nejbližší hodnotou.
## ID step angle x y time TEMP
## 1 Cheetah_1 NA NA 509.7917 -2792.354 2014-04-12 04:00:33 13.5
## 2 Cheetah_1 NA NA NA NA 2014-04-12 10:00:33 13.5
## 3 Cheetah_1 NA NA NA NA 2014-04-12 16:00:33 13.5
## 4 Cheetah_1 0.0000000 NA 507.7540 -2793.037 2014-04-12 22:00:33 25.5
## 5 Cheetah_1 0.2102280 NA 507.7540 -2793.037 2014-04-13 04:00:33 22.5
## 6 Cheetah_1 0.2130708 -3.137161 507.9374 -2792.934 2014-04-13 10:00:33 30.5
## Height geometry
## 1 1161 POINT (509791.7 -2792354)
## 2 1161 POINT EMPTY
## 3 1161 POINT EMPTY
## 4 1164 POINT (507754 -2793037)
## 5 1167 POINT (507754 -2793037)
## 6 1167 POINT (507937.4 -2792934)
Po použití funkce summary lze vidět v hlavičce, kolik je sledování pro daného jedince a ve spodní části je zobrazeno minimum, medián, průměr kovariát
library(dplyr)
## Warning: package 'dplyr' was built under R version 4.0.3
data %>% select(-c(geometry)) %>%
summary()
## Movement data for 1 animal:
## Cheetah_1 -- 802 observations
##
## Covariate(s):
## time
## Min. 25%
## "2014-04-12 04:00:33 CEST" "2014-06-01 05:30:33 CEST"
## Median Mean
## "2014-07-21 07:00:33 CEST" "2014-07-21 07:51:52 CEST"
## 75% Max.
## "2014-09-09 08:30:33 CEST" "2014-10-31 10:00:33 CET"
##
## TEMP
## Min. 25% Median Mean 75% Max.
## 0.00000 17.50000 23.50000 23.23192 29.37500 39.00000
##
## Height
## Min. 25% Median Mean 75% Max.
## 1088.000 1132.000 1161.000 1177.014 1201.000 1455.000
“plot(”prebData)" zobrazí:
mapu trajektorie jedinců
časové řady kroků a úhlů
histogramy kroků a úhlů
plot(data, compact = T)
Fit HMM je funkce, která vytvoří model pro skryté Markovovy modely, používající numericko-logarimtickou funkci.
Základní argumenty pro model:
data: moveData
nbStates: počet stavů HMM
stepPar0: počáteční vektor délky
anglePar0: počáteční vektor úhlů
formula: regresní formula pro kovariáty (výchozí: 1)
stepDist: název rozdělení pro délku kroků (výchozí: gama)
angleDist: název rozdělní pro úhel (výchozí: vm)
angleMean: vektro úhlů (výchozí: NULL)
stationary: argument FALSE použije se pokud existují kovariáty (výchozí: FALSE)
Dvoustavový model: Parametry pro délku kroků:
m0 = průměrná délka kroku
sigma0 = odchylka kroku
zeromass0 - použiju pokud se v datech nachází 0 délka kroku
StepPa0 = c(mu,sigma0,zeromass0)
mu0 <- c(0.5,2)
sigma0 <- c(0.5,1)
#kontrola jestli je potřeba zero-mass
whichzero <- which(data$step == 0)
length(whichzero)/nrow(data)
## [1] 0.02244389
#vyšlo číslo 0.02244389, které říká jaký je poměr 0 hodnot v datech
zeromass0 <- c(0.1,0.05)
stepPar0 <- c(mu0,sigma0,zeromass0)
Parametry pro velikost úhlů:
angleMean0 = průměrná velikost úhlů
kappa0 = odchylka pro velikost úhlu
anglePar0 <- c(angleMean0, kappa0)
angleMean0 <- c(pi,0)
kappa0 <- c(0.7,1.5)
anglePar0 <- c(angleMean0, kappa0)
Zvolené rozdělení:
stepDist - rozdělení pro délku kroku (“gamma” = gamma rozdělení)
angleDist - rozdělení pro úhly otáčení (“vm” = von Misesovo rozdělení)
stepDist <- "gamma"
angleDist <- "vm"
m1 <- fitHMM(data=data, nbStates = 2, stepPar0 = stepPar0, anglePar0 = anglePar0,
formula = ~1, stepDist=stepDist,
angleDist =angleDist)#žádná kovariata
m2 <- fitHMM(data=data, nbStates = 2, stepPar0 = stepPar0, anglePar0 = anglePar0,
formula = ~TEMP, stepDist=stepDist,
angleDist =angleDist)#kovariata -> teplota
m3 <- fitHMM(data=data, nbStates = 2, stepPar0 = stepPar0, anglePar0 = anglePar0,
formula = ~Height, stepDist=stepDist,
angleDist =angleDist)#"kovariata -> nadmořská výška
m4 <- fitHMM(data=data, nbStates = 2, stepPar0 = stepPar0, anglePar0 = anglePar0,
formula = ~TEMP+Height, stepDist=stepDist,
angleDist =angleDist)# kovarita TEMP + Height
Tří stavový model:
mu0 <- c(0.5,2,4)
sigma0 <- c(0.5,1,1.5)
whichzero <- which(data$step == 0)
length(whichzero)/nrow(data)
## [1] 0.02244389
zeromass0 <- c(0.1,0.05,0.01)
stepPar0 <- c(mu0,sigma0,zeromass0)
angleMean0 <- c(pi,pi/2,0)
kappa0 <- c(0.7,1.5,2)
anglePar0 <- c(angleMean0, kappa0)
m1_3 <- fitHMM(data=data, nbStates = 3, stepPar0 = stepPar0, anglePar0 = anglePar0,
formula = ~1, stepDist=stepDist,
angleDist =angleDist)# žádná kovariata
m3_3 <- fitHMM(data=data, nbStates = 3, stepPar0 = stepPar0, anglePar0 = anglePar0,
formula = ~Height, stepDist=stepDist,
angleDist =angleDist)# kovariata -> nadmořská výška
m4_3 <- fitHMM(data=data, nbStates = 3, stepPar0 = stepPar0, anglePar0 = anglePar0,
formula = ~TEMP+Height, stepDist=stepDist,
angleDist =angleDist)# kovarita TEMP + Height
Model s jedním stavem:
mu0 <- c(2)
sigma0 <- c(1)
whichzero <- which(data$step == 0)
length(whichzero)/nrow(data)
## [1] 0.02244389
zeromass0 <- c(0.05)
stepPar0 <- c(mu0,sigma0,zeromass0)
angleMean0 <- c(pi/2)
kappa0 <- c(1.5)
anglePar0 <- c(angleMean0, kappa0)
m0 <- fitHMM(data=data, nbStates = 1, stepPar0 = stepPar0, anglePar0 = anglePar0,
formula = ~1, stepDist=stepDist,
angleDist =angleDist)
Test AIC
print(AIC(m1,m2,m3,m4,m1_3,m3_3,m4_3,m0))
## Model AIC
## 1 m2 3348.546
## 2 m1_3 3354.094
## 3 m4 3357.687
## 4 m1 3360.814
## 5 m3 3366.062
## 6 m3_3 3368.870
## 7 m4_3 3416.735
## 8 m0 3453.976
Jako nejlepší model byl vybrán nakonec m2, protože má ze všech modelů nejnižší hodnotu AIC
m2
## Value of the maximum log-likelihood: -1659.273
##
## Step length parameters:
## ----------------------
## state 1 state 2
## mean 0.31415220 1.76045768
## sd 0.50432602 1.29899567
## zero-mass 0.01419297 0.04005698
##
## Turning angle parameters:
## ------------------------
## state 1 state 2
## mean -2.92956218 0.004325341
## concentration 0.09265775 0.776972965
##
## Regression coeffs for the transition probabilities:
## --------------------------------------------------
## 1 -> 2 2 -> 1
## intercept 0.7821853 -1.13322631
## TEMP -0.1107295 0.01595603
##
## Initial distribution:
## --------------------
## [1] 1.327189e-06 9.999987e-01
CI(m2)
## $stepPar
## $stepPar$lower
## state 1 state 2
## mean 0.223120888 1.51845178
## sd 0.358965707 1.12890130
## zero-mass 0.003069412 0.01454566
##
## $stepPar$upper
## state 1 state 2
## mean 0.44232346 2.041034
## sd 0.70854885 1.494719
## zero-mass 0.06307772 0.105521
##
##
## $anglePar
## $anglePar$lower
## state 1 state 2
## mean -5.36278696 -0.2706063
## concentration 0.02378886 0.5098092
##
## $anglePar$upper
## state 1 state 2
## mean -0.8072294 0.2976098
## concentration 0.2612791 1.0599034
##
##
## $beta
## $beta$lower
## 1 -> 2 2 -> 1
## intercept -0.4872198 -2.59643054
## TEMP -0.1760870 -0.04219851
##
## $beta$upper
## 1 -> 2 2 -> 1
## intercept 2.05159036 0.32997791
## TEMP -0.04537194 0.07411056
V tabulce CI lze vidět:
interval po střední délku kroku -> pro stav 1 = (0.223207117-0.359104588) m; pro stav 2 = (1.51905094-2.041851) m
interval pro variabliní délku kroku -> pro stav 1 = 0.359104588 m; pro stav 2 = 1.12934926 m
střední velikost úhlu -> pro stav 1 = (-5.35580872 do -0.8082987) -> pro stav 2 = (-0.2698473 -0.2968336)
koncetrace úhlů -> pro stav 1 = 0.1948008
-> pro stav 2 = 0.61492293
pravděpodobnost přechodu z jednoho stavu do druhého a na zpět
plot(m2, plotCI = T)
## Decoding states sequence... DONE
sp <- stateProbs(m2)
head(sp)
## [,1] [,2]
## [1,] 1.301543e-06 0.999998698
## [2,] 2.751819e-01 0.724818136
## [3,] 3.643776e-01 0.635622361
## [4,] 4.894820e-01 0.510517997
## [5,] 8.992283e-01 0.100771676
## [6,] 9.905691e-01 0.009430944
plotStates(m2)
## Decoding states sequence... DONE
## Computing states probabilities... DONE
plotStationary(m2, plotCI=TRUE)
Stacionární stav bude nejpravděpodobnější mezi 15-20°C.
res <- pseudoRes(m2)
plotPR(m2)
## Computing pseudo-residuals... DONE
Test normality
library(tseries)
## Warning: package 'tseries' was built under R version 4.0.3
#test normality kroků
## Jarqueův a Berryho test normality
print(jarque.bera.test(res$step[which(!is.na(res$step))]))
##
## Jarque Bera Test
##
## data: res$step[which(!is.na(res$step))]
## X-squared = 18.171, df = 2, p-value = 0.0001133
## Shapirro-Wilkův test normality
shapiro.test(res$stepRes)
##
## Shapiro-Wilk normality test
##
## data: res$stepRes
## W = 0.97946, p-value = 6.616e-09
#test normality pro úhly
print(jarque.bera.test(res$angle[which(!is.na(res$angle))]))
##
## Jarque Bera Test
##
## data: res$angle[which(!is.na(res$angle))]
## X-squared = 357.67, df = 2, p-value < 2.2e-16
shapiro.test(res$angleRes)
##
## Shapiro-Wilk normality test
##
## data: res$angleRes
## W = 0.9639, p-value = 2.777e-12
Test korelace
par(mfrow=c(1,2))
acf(na.omit((data$step)))
acf(na.omit((data$angle)))
get_closest_id <- function(i, ids=gepard_D$ID){
id <- ids[i]
ii <- i
while (is.na(id) & ii > 1){
ii <- ii - 1
id <- ids[ii]
}
ii <- i
while (is.na(id) & ii < length(ids)){
ii <- ii + 1
id <- ids[ii]
}
return(id)
}
for (i in which(is.na(gepard_D$ID))) {
gepard_D$ID[i] <- get_closest_id(i)
}
data_D <- prepData(gepard_D,type="UTM",coordNames=c("X","Y"))
## Warning in prepData(gepard_D, type = "UTM", coordNames = c("X", "Y")): There
## are 4 missing covariate values. Each will be replaced by the closest available
## value.
plot(data_D, compact = T)
Model:
# Parametry 2 stavy
mu0 <- c(0.2,1)
sigma0 <- c(0.5,1)
whichzero <- which(data_D$step == 0)
length(whichzero)/nrow(data_D)
## [1] 0.0375
zeromass0 <- c(0.1,0.05)
stepPar0 <- c(mu0,sigma0,zeromass0)
angleMean0 <- c(pi,0)
kappa0 <- c(0.7,1.5)
anglePar0 <- c(angleMean0, kappa0)
stepDist <- "gamma"
angleDist <- "vm"
# Model 2 stavy
m2_D <- fitHMM(data=data_D, nbStates = 2, stepPar0 = stepPar0, anglePar0 = anglePar0,
formula = ~TEMP, stepDist=stepDist,
angleDist =angleDist)
m1_D <- fitHMM(data=data_D, nbStates = 2, stepPar0 = stepPar0, anglePar0 = anglePar0,
formula = ~Height, stepDist=stepDist,
angleDist =angleDist)
m3_D <- fitHMM(data=data_D, nbStates = 2, stepPar0 = stepPar0, anglePar0 = anglePar0,
formula = ~TEMP+Height, stepDist=stepDist,
angleDist =angleDist)# kovarita TEMP + Height
m4_D <- fitHMM(data=data_D, nbStates = 2, stepPar0 = stepPar0, anglePar0 = anglePar0,
formula = ~1, stepDist=stepDist,
angleDist =angleDist)
# Parametry 3 stavy
mu0 <- c(0.5,2,4)
sigma0 <- c(0.5,1,1.5)
whichzero <- which(data$step == 0)
length(whichzero)/nrow(data)
## [1] 0.02244389
zeromass0 <- c(0.1,0.05,0.01)
stepPar0 <- c(mu0,sigma0,zeromass0)
angleMean0 <- c(pi,pi/2,0)
kappa0 <- c(0.7,1.5,2)
anglePar0 <- c(angleMean0, kappa0)
m1_3_D <- fitHMM(data=data_D, nbStates = 3, stepPar0 = stepPar0, anglePar0 = anglePar0,
formula = ~1, stepDist=stepDist,
angleDist =angleDist)#žádná kovariata
m3_3_D <- fitHMM(data=data_D, nbStates = 3, stepPar0 = stepPar0, anglePar0 = anglePar0,
formula = ~Height, stepDist=stepDist,
angleDist =angleDist)#"kovariata -> nadmořská výška
m2_3_D <- fitHMM(data=data_D, nbStates = 3, stepPar0 = stepPar0, anglePar0 = anglePar0,
formula = ~TEMP, stepDist=stepDist,
angleDist =angleDist)#"kovariata -> nadmořská výška
m4_3_D <- fitHMM(data=data_D, nbStates = 3, stepPar0 = stepPar0, anglePar0 = anglePar0,
formula = ~TEMP+Height, stepDist=stepDist,
angleDist =angleDist)# kovarita TEMP + Height
# Parametry 1 stav
mu0 <- c(1)
sigma0 <- c(0.5)
whichzero <- which(data_D$step == 0)
length(whichzero)/nrow(data_D)
## [1] 0.0375
zeromass0 <- c(0.1)
stepPar0 <- c(mu0,sigma0,zeromass0)
angleMean0 <- c(pi)
kappa0 <- c(0.7)
anglePar0 <- c(angleMean0, kappa0)
stepDist <- "gamma"
angleDist <- "vm"
# Model 1 stav
m0_D <- fitHMM(data=data_D, nbStates = 1, stepPar0 = stepPar0, anglePar0 = anglePar0,
formula = ~TEMP, stepDist=stepDist,
angleDist =angleDist)
AIC
print(AIC(m2_D,m0_D,m1_D,m4_D, m1_3_D,m2_3_D,m3_3_D, m4_3_D))
## Model AIC
## 1 m4_D -97.39544
## 2 m0_D 253.74059
## 3 m2_D 265.41650
## 4 m1_D 266.98672
## 5 m1_3_D 274.46164
## 6 m2_3_D 301.74059
## 7 m3_3_D 301.74059
## 8 m4_3_D 313.74059
m4_D
## Value of the maximum log-likelihood: 61.69772
##
## Step length parameters:
## ----------------------
## state 1 state 2
## mean 0.003888028 5.684876e-01
## sd 0.002003555 7.439935e-01
## zero-mass 0.067995624 5.768532e-23
##
## Turning angle parameters:
## ------------------------
## state 1 state 2
## mean 2.7691230 -2.962851
## concentration 0.4894726 0.349154
##
## Regression coeffs for the transition probabilities:
## --------------------------------------------------
## 1 -> 2 2 -> 1
## intercept -1.670741 -1.311827
##
## Transition probability matrix:
## -----------------------------
## [,1] [,2]
## [1,] 0.8416745 0.1583255
## [2,] 0.2121813 0.7878187
##
## Initial distribution:
## --------------------
## [1] 8.815086e-09 1.000000e+00
CI(m4_D)
## Warning in CI(m4_D): Some of the parameter estimates seem to lie close to the boundaries of their parameter space.
## The associated CIs are probably unreliable (or might not be computable).
## $stepPar
## $stepPar$lower
## state 1 state 2
## mean 0.003446767 4.330280e-01
## sd 0.001596398 5.632207e-01
## zero-mass 0.035577869 5.768532e-23
##
## $stepPar$upper
## state 1 state 2
## mean 0.004385781 7.463216e-01
## sd 0.002514556 9.827875e-01
## zero-mass 0.126089830 5.768532e-23
##
##
## $anglePar
## $anglePar$lower
## state 1 state 2
## mean 2.1574773 -3.9943596
## concentration 0.2367741 0.1122786
##
## $anglePar$upper
## state 1 state 2
## mean 3.3617254 -1.9924064
## concentration 0.7872902 0.6735214
##
##
## $beta
## $beta$lower
## 1 -> 2 2 -> 1
## intercept -2.192615 -1.949317
##
## $beta$upper
## 1 -> 2 2 -> 1
## intercept -1.148867 -0.6743373
plot(m4_D, plotCI = T)
## Decoding states sequence... DONE
sp <- stateProbs(m4_D)
head(sp)
## [,1] [,2]
## [1,] 3.636281e-18 1.0000000
## [2,] 3.721885e-220 1.0000000
## [3,] 0.000000e+00 1.0000000
## [4,] 5.948822e-08 0.9999999
## [5,] 1.281626e-09 1.0000000
## [6,] 0.000000e+00 1.0000000
plotStates(m4_D)
## Decoding states sequence... DONE
## Computing states probabilities... DONE
res <- pseudoRes(m4_D)
## Note: Some angles are equal to pi, and the corresponding pseudo-residuals are not included
plotPR(m4_D)
## Computing pseudo-residuals...
## Note: Some angles are equal to pi, and the corresponding pseudo-residuals are not included
## DONE
Test normality
library(tseries)
#test normality kroků
## Jarqueův a Berryho test normality
print(jarque.bera.test(res$step[which(!is.na(res$step))]))
##
## Jarque Bera Test
##
## data: res$step[which(!is.na(res$step))]
## X-squared = 6.6987, df = 2, p-value = 0.03511
## Shapirro-Wilkův test normality
shapiro.test(res$stepRes)
##
## Shapiro-Wilk normality test
##
## data: res$stepRes
## W = 0.98045, p-value = 0.002479
#test normality pro úhly
print(jarque.bera.test(res$angle[which(!is.na(res$angle))]))
##
## Jarque Bera Test
##
## data: res$angle[which(!is.na(res$angle))]
## X-squared = 115.46, df = 2, p-value < 2.2e-16
shapiro.test(res$angleRes)
##
## Shapiro-Wilk normality test
##
## data: res$angleRes
## W = 0.93971, p-value = 1.354e-07
Test korelace
par(mfrow=c(1,2))
acf(na.omit((data_D$step)))
acf(na.omit((data_D$angle)))
get_closest_id <- function(i, ids=gepard_PD$ID){
id <- ids[i]
ii <- i
while (is.na(id) & ii > 1){
ii <- ii - 1
id <- ids[ii]
}
ii <- i
while (is.na(id) & ii < length(ids)){
ii <- ii + 1
id <- ids[ii]
}
return(id)
}
for (i in which(is.na(gepard_PD$ID))) {
gepard_PD$ID[i] <- get_closest_id(i)
}
data_PD <- prepData(gepard_PD,type="UTM",coordNames=c("X","Y"))
## Warning in prepData(gepard_PD, type = "UTM", coordNames = c("X", "Y")): There
## are 34 missing covariate values. Each will be replaced by the closest available
## value.
plot(data_PD, compact = T)
Model:
# Parametry 2 stavy
mu0 <- c(1,2)
sigma0 <- c(0.5,1)
whichzero <- which(data_D$step == 0)
length(whichzero)/nrow(data_D)
## [1] 0.0375
zeromass0 <- c(0.1,0.05)
stepPar0 <- c(mu0,sigma0,zeromass0)
angleMean0 <- c(pi,0)
kappa0 <- c(0.7,1.5)
anglePar0 <- c(angleMean0, kappa0)
stepDist <- "gamma"
angleDist <- "vm"
# Model 2 stavy
m2_PD <- fitHMM(data=data_PD, nbStates = 2, stepPar0 = stepPar0, anglePar0 = anglePar0,
formula = ~TEMP, stepDist=stepDist,
angleDist =angleDist)
m1_PD <- fitHMM(data=data_PD, nbStates = 2, stepPar0 = stepPar0, anglePar0 = anglePar0,
formula = ~Height, stepDist=stepDist,
angleDist =angleDist)
m3_PD <- fitHMM(data=data_PD, nbStates = 2, stepPar0 = stepPar0, anglePar0 = anglePar0,
formula = ~Height+TEMP, stepDist=stepDist,
angleDist =angleDist)
m4_PD <- fitHMM(data=data_PD, nbStates = 2, stepPar0 = stepPar0, anglePar0 = anglePar0,
formula = ~1, stepDist=stepDist,
angleDist =angleDist)
# Parametry 3 stavy
mu0 <- c(0.5,2,4)
sigma0 <- c(0.5,1,1.5)
whichzero <- which(data$step == 0)
length(whichzero)/nrow(data)
## [1] 0.02244389
zeromass0 <- c(0.1,0.05,0.01)
stepPar0 <- c(mu0,sigma0,zeromass0)
angleMean0 <- c(pi,pi/2,0)
kappa0 <- c(0.7,1.5,2)
anglePar0 <- c(angleMean0, kappa0)
# Model 3 stavy
m1_3_PD <- fitHMM(data=data_PD, nbStates = 3, stepPar0 = stepPar0, anglePar0 = anglePar0,
formula = ~1, stepDist=stepDist,
angleDist =angleDist)#žádná kovariata
m3_3_PD <- fitHMM(data=data_PD, nbStates = 3, stepPar0 = stepPar0, anglePar0 = anglePar0,
formula = ~Height, stepDist=stepDist,
angleDist =angleDist)#"kovariata -> nadmořská výška
m2_3_PD <- fitHMM(data=data_PD, nbStates = 3, stepPar0 = stepPar0, anglePar0 = anglePar0,
formula = ~TEMP, stepDist=stepDist,
angleDist =angleDist)#"kovariata -> teplota
m4_3_PD <- fitHMM(data=data_PD, nbStates = 3, stepPar0 = stepPar0, anglePar0 = anglePar0,
formula = ~TEMP+Height, stepDist=stepDist,
angleDist =angleDist)# kovarita TEMP + Height
# Parametry 1 stav
mu0 <- c(0.5)
sigma0 <- c(0.5)
whichzero <- which(data_D$step == 0)
length(whichzero)/nrow(data_D)
## [1] 0.0375
zeromass0 <- c(0.1)
stepPar0 <- c(mu0,sigma0,zeromass0)
angleMean0 <- c(pi)
kappa0 <- c(0.7)
anglePar0 <- c(angleMean0, kappa0)
stepDist <- "gamma"
angleDist <- "vm"
# Model 1 stav
m0_PD <- fitHMM(data=data_PD, nbStates = 1, stepPar0 = stepPar0, anglePar0 = anglePar0,
formula = ~TEMP, stepDist=stepDist,
angleDist =angleDist)
AIC
print(AIC(m2_PD,m1_PD,m4_PD,m0_PD,m1_3_PD,m2_3_PD,m3_3_PD,m4_3_PD))
## Model AIC
## 1 m2_PD 467.6710
## 2 m0_PD 487.3513
## 3 m4_PD 503.3513
## 4 m1_PD 507.3513
## 5 m1_3_PD 514.6504
## 6 m2_3_PD 535.3513
## 7 m3_3_PD 535.3513
## 8 m4_3_PD 547.3513
m2_PD
## Value of the maximum log-likelihood: -218.8355
##
## Step length parameters:
## ----------------------
## state 1 state 2
## mean 0.13209784 1.284182e+00
## sd 0.18587667 8.402803e-01
## zero-mass 0.01099041 6.373658e-07
##
## Turning angle parameters:
## ------------------------
## state 1 state 2
## mean -2.9651941 -1.4476360
## concentration 0.1801042 0.3700251
##
## Regression coeffs for the transition probabilities:
## --------------------------------------------------
## 1 -> 2 2 -> 1
## intercept 14.2203231 11.5369488
## TEMP -0.5674046 -0.3356778
##
## Initial distribution:
## --------------------
## [1] 1 0
CI(m2_PD)
## Warning in CI(m2_PD): Some of the parameter estimates seem to lie close to the boundaries of their parameter space.
## The associated CIs are probably unreliable (or might not be computable).
## $stepPar
## $stepPar$lower
## state 1 state 2
## mean 0.088705201 1.009076e+00
## sd 0.120913691 6.230168e-01
## zero-mass 0.001541984 6.370849e-07
##
## $stepPar$upper
## state 1 state 2
## mean 0.19671719 1.634290e+00
## sd 0.28574212 1.133310e+00
## zero-mass 0.07404036 6.376468e-07
##
##
## $anglePar
## $anglePar$lower
## state 1 state 2
## mean -5.37328150 -2.74611578
## concentration 0.04917076 0.09589523
##
## $anglePar$upper
## state 1 state 2
## mean -0.7989543 1.996592
## concentration 0.5288348 0.889244
##
##
## $beta
## $beta$lower
## 1 -> 2 2 -> 1
## intercept 4.1309013 0.9470310
## TEMP -0.9498143 -0.6413018
##
## $beta$upper
## 1 -> 2 2 -> 1
## intercept 24.309745 22.12686659
## TEMP -0.184995 -0.03005385
plot(m2_PD, plotCI = T)
## Decoding states sequence... DONE
sp <- stateProbs(m2_PD)
head(sp)
## [,1] [,2]
## [1,] 1.000000e+00 0.000000e+00
## [2,] 8.133254e-06 9.999919e-01
## [3,] 9.836618e-04 9.990163e-01
## [4,] 9.987902e-01 1.209764e-03
## [5,] 9.999920e-01 8.032164e-06
## [6,] 1.095280e-02 9.890472e-01
plotStates(m2_PD)
## Decoding states sequence... DONE
## Computing states probabilities... DONE
plotStationary(m2_PD, plotCI=TRUE)
res <- pseudoRes(m2_PD)
plotPR(m2_PD)
## Computing pseudo-residuals... DONE
Test normality
library(tseries)
#test normality kroků
print(jarque.bera.test(res$step[which(!is.na(res$step))]))
##
## Jarque Bera Test
##
## data: res$step[which(!is.na(res$step))]
## X-squared = 2.2862, df = 2, p-value = 0.3188
## Shapirro-Wilkův test normality
shapiro.test(res$stepRes)
##
## Shapiro-Wilk normality test
##
## data: res$stepRes
## W = 0.97851, p-value = 0.03119
#test normality pro úhly
print(jarque.bera.test(res$angle[which(!is.na(res$angle))]))
##
## Jarque Bera Test
##
## data: res$angle[which(!is.na(res$angle))]
## X-squared = 402.61, df = 2, p-value < 2.2e-16
shapiro.test(res$angleRes)
##
## Shapiro-Wilk normality test
##
## data: res$angleRes
## W = 0.89378, p-value = 4.03e-08
Test korelace
par(mfrow=c(1,2))
acf(na.omit((data_PD$step)))
acf(na.omit((data_PD$angle)))
library(trajr)
trj1<-read.table("Cheetah_1.csv", dec=",",sep = ";", head=T)
trj1$time=as.Date(trj1$time)
trj2<-read.table("Cheetah_2.csv", dec=",",sep = ";", head=T)
trj2$time=as.Date(trj2$time)
trj3<-read.table("Cheetah_3.csv", dec=",",sep = ";", head=T)
trj3$time=as.Date(trj3$time)
Ch_1 <- TrajFromCoords(trj1, xCol = "x", yCol = "y", timeCol = "time")
Ch_2 <- TrajFromCoords(trj2, xCol = "x", yCol = "y", timeCol = "time")
Ch_3 <- TrajFromCoords(trj3, xCol = "x", yCol = "y", timeCol = "time")
p=9
n=p + 3 - p%%2
Cheetah_1 <- TrajSmoothSG(Ch_1, p, n)
Cheetah_2 <- TrajSmoothSG(Ch_2, p, n)
Cheetah_3 <- TrajSmoothSG(Ch_3, p, n)
library(leaflet)
map = leaflet() %>%
addTiles() %>%
addPolylines(data = Cheetah_1, lat = ~y, lng = ~x, color = "orange",weight =2, group = )%>%
addPolylines(data = Cheetah_2, lat = ~y, lng = ~x, color = "blue", weight= 2, group = )%>%
addPolylines(data = Cheetah_3, lat = ~y, lng = ~x, color = "green", weight = 2, group = )%>%
addLegend(position = "topright", labels = c("Samice po vypuštění", "Těhotná samice", "Samice s mláďaty"),
colors = c("orange","blue","green"))
map