jmsigner / amt

38 stars 13 forks source link

NULLs for step length (sl_) and turn angle (ta_) distributions in fit_issf() output #79

Closed kpiecora closed 2 years ago

kpiecora commented 2 years ago

Hi,

My fit_issf() output has NULLs instead of returning step length and turn angle distributions. I modeled my code from "Appendix B: SSF Examples" by Fieberg et al and am able to run the fisher example all the way through without any issues. When I attempt "sl_distr()" it returns the following error: "Error in UseMethod("sl_distr") : no applicable method for 'sl_distr' applied to an object of class "data.frame", yet no error when doing the same with the example fisher data frame. Below is a small subset of data that reproduces the NULLs I'm seeing with my full df.

Any suggestions?

### Load libraries

library(amt)
#> 
#> Attaching package: 'amt'
#> The following object is masked from 'package:stats':
#> 
#>     filter
library(reprex)

### Load data

dat <- data.frame(
  stringsAsFactors = FALSE,
  id = c(2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L),
  burst_ = c(1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L),
  x1_ = c(737577.6672,737577.6672,
          737577.6672,737577.6672,737577.6672,737577.6672,
          737579.5066,737579.5066,737579.5066,737579.5066),
  x2_ = c(737579.5066,736838.2145,
          737706.4363,737306.0804,737528.273,737692.0972,737757.4911,
          737507.2334,739510.3222,737107.5565),
  y1_ = c(4585896.991,4585896.991,
          4585896.991,4585896.991,4585896.991,4585896.991,4585908.79,
          4585908.79,4585908.79,4585908.79),
  y2_ = c(4585908.79,4588583.422,
          4586008.427,4585353.744,4586070.49,4585681.683,4585328.705,
          4584701.401,4585904.357,4586007.453),
  sl_ = c(11.94075569,2786.341859,
          170.2920436,607.3524,180.39332,243.8272835,606.7755773,
          1209.549363,1930.820668,482.1529877),
  ta_ = c(2.902262515,0.423265156,
          -0.702772055,2.832661421,0.432018573,-2.498424097,
          -2.689224085,-0.357494216,1.270794122,-2.074590289),
  t1_ = c("2/19/2021","2/19/2021",
          "2/19/2021","2/19/2021","2/19/2021","2/19/2021",
          "2/19/2021","2/19/2021","2/19/2021","2/19/2021"),
  t2_ = c("2/19/2021","2/19/2021",
          "2/19/2021","2/19/2021","2/19/2021","2/19/2021",
          "2/19/2021","2/19/2021","2/19/2021","2/19/2021"),
  dt_ = c(2.5, 2.5, 2.5, 2.5, 2.5, 2.5, 2.5, 2.5, 2.5, 2.5),
  step_id_ = c(2496L,2496L,2496L,2496L,
               2496L,2496L,2497L,2497L,2497L,2497L),
  case_ = c(1L, 0L, 0L, 0L, 0L, 0L, 1L, 0L, 0L, 0L),
  PercCrop = c(51.4739229,0,68.93424036,
               20.1814059,41.26984127,60.77097506,40.36281179,
               0.453514739,5.442176871,0),
  PercWetl = c(0, 1.133786848, 0, 0, 0, 0, 0, 0, 0.680272109, 0),
  cos_ta_ = c(-0.971496984730054,
              0.911752673860636,0.763053444078309,-0.952659062644458,
              0.908122414247981,-0.80019948708203,-0.899414331808547,
              0.936776607666766,0.295522312981637,-0.482751595231058),
  log_sl_ = c(2.47995739658145,
              7.93248485284213,5.13751486667119,6.40910918273146,5.19513957811811,
              5.49646012013771,6.40815899831238,7.0980031419498,
              7.56570040810029,6.1782614655565)
)

# Model

m0 <- dat %>% 
  fit_issf(case_ ~ PercCrop + PercWetl + sl_ + log_sl_ + cos_ta_ + strata(step_id_), model = TRUE)
#> Warning in coxexact.fit(X, Y, istrat, offset, init, control, weights =
#> weights, : Ran out of iterations and did not converge

summary(m0)
#> Call:
#> coxph(formula = Surv(rep(1, 10L), case_) ~ PercCrop + PercWetl + 
#>     sl_ + log_sl_ + cos_ta_ + strata(step_id_), data = data, 
#>     model = TRUE, method = "exact")
#> 
#>   n= 10, number of events= 2 
#> 
#>                coef  exp(coef)   se(coef)      z Pr(>|z|)
#> PercCrop  4.497e-01  1.568e+00  4.603e+02  0.001    0.999
#> PercWetl -2.512e+01  1.234e-11  3.765e+04 -0.001    0.999
#> sl_       2.307e-02  1.023e+00  1.813e+01  0.001    0.999
#> log_sl_  -1.010e+01  4.119e-05  1.330e+04 -0.001    0.999
#> cos_ta_  -5.327e+00  4.859e-03  1.342e+04  0.000    1.000
#> 
#>          exp(coef) exp(-coef) lower .95 upper .95
#> PercCrop 1.568e+00  6.378e-01 0.000e+00       Inf
#> PercWetl 1.234e-11  8.102e+10 0.000e+00       Inf
#> sl_      1.023e+00  9.772e-01 3.797e-16 2.758e+15
#> log_sl_  4.119e-05  2.428e+04 0.000e+00       Inf
#> cos_ta_  4.859e-03  2.058e+02 0.000e+00       Inf
#> 
#> Concordance= 1  (se = 0 )
#> Likelihood ratio test= 6.36  on 5 df,   p=0.3
#> Wald test            = 0  on 5 df,   p=1
#> Score (logrank) test = 6.66  on 5 df,   p=0.2

Created on 2022-11-02 with reprex v2.0.2

jmsigner commented 2 years ago

How did you create dat? The function random_steps() fits a distribution to the step lengths. If you did not use this, then the step-length distribution is not saved. Maybe you could give a few more details on your work flow prior to fit_issf().

kpiecora commented 2 years ago

I created it with random_steps() first, but then saved as a csv to extract habitat variables in ArcGIS (my computer couldn't handle extracting from rasters within R). Which means that when I brought the file back into R it was only a data frame. Would the tibble originally created with random_steps() have contained the distributions which were then lost when I exported it as a csv?

jmsigner commented 2 years ago

Yes, when exporting the to .csv file all attributes are lost. You can extract the step-length (and also turn-angle) distribution prior to the export using the same function and save them to a variable for later.

kpiecora commented 2 years ago

Makes sense. Thanks for helping me get to the bottom of it!