moveHMM - použití HMM na samici geparda štíhlého

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).

Instalce balíčku, načtení dat a jejich úprava

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

Model pro geparda po vypuštění

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

Vizualizace “prebData”

“plot(”prebData)" zobrazí:

mapu trajektorie jedinců
časové řady kroků a úhlů
histogramy kroků a úhlů

plot(data, compact = T)

Model fitHMM

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

Konfidenční intervaly

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

Vizualizace modelu

plot(m2, plotCI = T)

## Decoding states sequence... DONE

Pravděpodobnosti stavu

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

Pravděpodobnost stacionárncího stavu

plotStationary(m2, plotCI=TRUE)

Stacionární stav bude nejpravděpodobnější mezi 15-20°C.

Kontrola modelu

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)))

Model pro březího geparda

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

Konfidenční intervaly

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

Vizualizace modelu

plot(m4_D, plotCI = T)

## Decoding states sequence... DONE

Pravděpodobnosti stavu

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

Kontrola modelu

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)))

Model pro geparda s mláďaty v poporodní fázi

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

Konfidenční intervaly

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

Vizualizace modelu

plot(m2_PD, plotCI = T)

## Decoding states sequence... DONE

Pravděpodobnosti stavu

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

Pravděpodobnost stacionárního stavu

plotStationary(m2_PD, plotCI=TRUE)

Kontrola modelu

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)))

Použití balíčku trajr na vyhlazení trajektorie a její nasledné zobrazení pomocí leaflet

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

moveHMM

Kadlec Ivo

18 2 2021

moveHMM - použití HMM na samici geparda štíhlého

Instalce balíčku, načtení dat a jejich úprava

Model pro geparda po vypuštění

Vizualizace “prebData”

Model fitHMM

Konfidenční intervaly

Vizualizace modelu

Pravděpodobnosti stavu

Pravděpodobnost stacionárncího stavu

Kontrola modelu

Model pro březího geparda

Konfidenční intervaly

Vizualizace modelu

Pravděpodobnosti stavu

Kontrola modelu

Model pro geparda s mláďaty v poporodní fázi

Konfidenční intervaly

Vizualizace modelu

Pravděpodobnosti stavu

Pravděpodobnost stacionárního stavu

Kontrola modelu

Použití balíčku trajr na vyhlazení trajektorie a její nasledné zobrazení pomocí leaflet