Trouble interpreting null model approach for artefactual periodicity test

[Holly Gamblin]

Hi there, 

I'm interested in assessing space use and resource selection of Arctic and red foxes during the breeding season, and I have a data set containing GPS data collected at a frequency of 1 fix/15 minutes between April - August 2023. I think it'd be really interesting to look at periodic patterns of space use, but I'm a little concerned that my sampling schedule produces artefactual periodicity. The sampling schedule rotates on a 4-day cycle, collecting 1 fix/15 min for 6 hours each day so that in 4 days all 24 hours are covered (i.e., day 1 hours = 000-0600, day 2 hours = 0600-1200, day 3 hours = 1200-1800, day 4 hours = 1800-2400).

When I did a visual assessment of variograms/periodograms for two foxes, Randy and Woody, I had trouble interpreting what I was seeing. It appears Woody exhibits range residency from the variogram, although Randy is a different story. I'm not sure why the best fitting model for Randy varies so much from the variogram, but I don't suspect Randy is resident. To give some background on Randy's movements, it appears Randy was captured and collared during either an exploratory movement or a dispersal, as the next two days after capture it moves and settles ~25 km from the capture location. It stays in this area for ~ 1 month and then disperses northwest of the study area, where it appears to have remained and established a home range. I'm not sure why the variogram doesn't reflect that pattern. 

The sampling schedule periodograms (red) seem very similar to Woody and Randy's periodograms (black), so I suspected this might be indicative of artefactual periodicity. I attempted to use the code in Additional File 3 from Péron et al. (2016), but I got the same output for both foxes when I ran the period.test function: 

> period.test(Woody, Woody.OUF, n = 100) # period (days) = 14745600/(86400) = 170, p = 0
> period.test(Randy, Randy.OUF, n = 100) # period (days) = 14745600/(86400) = 170, p = 0

Here is a link to the R script and data I used to generate all models, figures, and tests: https://github.com/gamblinh/FoxData_Periodicity

Is it a problem with my code? Am I overlooking some glaring issue in the model diagnostics? 

Also, as more of a basic question, how exactly does one interpret a periodogram? I tried following the Péron et al. (2016) and the Calabrese et al. (2016) paper (and associated appendices, such as Appendix C), and I'm still struggling to interpret the graph. I understand vertical lines and peaks indicate the period of repeated space use, but how exactly do I interpret the graph? What are those little black ticks on the x-axis? How do I tease out the exact hour/day/month that the periodic pattern occurs? My apologies for this basic question.

Thank you so much for your help!

Woody variogram with OUF model overlayed

Woody periodogram 

Randy variogram with OUF model overlayed

Randy periodogram

[Christen Fleming]

Hi Holly,

Since that paper was published, ctmm.select() was updated with an appropriate model selection criterion for that model (based on AIC/BIC + MSPE). Presently, you should just be able to call ctmm.select() and go with the output of that.

The x-axis of the variogram is "time lag" between two times and not absolute time. It will not reflect the schedule over which the data were collected, but here Randy's variogram seems to indicate resident phase(s) of the data about ~1 week long, within a much larger dispersal range on a seasonal scale, while Woody's variogram is mostly resident day-to-day, but slightly shifted over the seasonal scale.

If the data are evenly sampled, then the peaks of the periodogram correspond to the periods and harmonics (integer fractions of the period) of repeating space use. The tick marks correspond to meaningful periods (e.g., a day) and their potential harmonics (e.g., 1/2 a day, 1/3 a day, 1/4 a day, ...). I will update the help file with this information. Here you have peaks corresponding to the 4-day, 1-day and 6-hr scheduling, so you can't really see anything in the periodogram of the data (black) that isn't in the periodogram of the schedule (red).

Best,

Chris

[Holly Gamblin]

Hi, 

Thank you so much for your explanation! I was worried the peaks were more related to the sampling schedule than true periodicity. If there is artefactual periodicity from the schedule, is there an alternative way to get at the question of periodic patterns of space use? 

I did have one follow-up question in regard to null hypothesis testing to confirm artefactual periodicity. In my original code I had used the ctmm.select() function, which indicated OUF was the best model for both foxes. The models in the variograms are from the ctmm.select() function. When I ran the period.test function with the fitted OUF models, the values for Woody and Randy were identical and it returned P = 0, which seems suspicious considering the periodograms seem to indicate artefactual periodicity. I've pasted a shortened version of the code below for an example (full code and data available at https://github.com/gamblinh/FoxData_Periodicity): 

# RF Woody # 

#Calculate the variogram

Woody.vg<-variogram(Woody)

plot(Woody.vg)

 

#Calculate the periodogram

Woody.pg <- periodogram(Woody)

plot(Woody.pg, diagnostic = TRUE)

 

#Use the sliders provided by variogram.fit to specify starting values.

variogram.fit(Woody.vg)

#save the parameters in the manipulator window to GUESS

 

#Automatically fit the range-resident models via maximum likelihood

Woody.fit.mods<-ctmm.select(Woody,CTMM=GUESS, verbose=TRUE)

summary(Woody.fit.mods)

 

# OUF Anisotropic is top model for Woody

Woody.OUF<-Woody.fit.mods[[1]]

summary(Woody.OUF)

plot(Woody.vg, CTMM = Woody.OUF)

 

# Note that AF Randy OUF calculated using similar code to Woody 

############################################################

# Function from Appendix

############################################################

period.test=function(data,model,period=NULL,n=100,deltat=60,...) {

  # initiate output vector

  RES=NULL

  # compute periodogram of the true data

  per=periodogram(data,...)

  # set period of interest at the periodogram peak if not provided by user

  if(is.null(period)) period=1/per$f[which.max(per$LSP)]

  else period=period[1]

  # MC loop

  for (k in 1:n) {

    # simulate dataset using model provided by user and time grid from the true data

    s=simulate(object=model,t=data$t-data$t[1])

    # compute periodogram of the simulated dataset

    s.per=periodogram(s,...)

    # compute whether simulated value is higher than true value

    RES=c(RES, max(s.per$LSP[abs(1/s.per$f - period)<=deltat])<

            max(per$LSP[abs(1/per$f - period)<=deltat]) )

  }

  return(list(period=period,P.value=1-sum(RES)/sum(!is.na(RES))))

}

####################################################################

#returns a list with two items: tested period in seconds and P-value of that period

period.test(Woody, Woody.OUF, n = 1000) # period (days) = 14745600/(3600*24) = 170,  p = 0

period.test(Randy, Randy.OUF, n = 1000) # period(days) = 14745600/(3600*24) = 170, p = 0

[Holly Gamblin]

As one more follow-up to the previous code regarding the reported output from the period.test, it appears the period (in days) for both foxes is 170 days, which is greater than the number of days we have the data for (~ April 18 - August 31 = 135 days). Why would the detected period be longer than the total period of the track? The identical output that extends beyond the survey period for both foxes is puzzling. 

[Christen Fleming]

Hi Holly,

It looks like you aren't fitting any periodic models in your code. I would take a look a the periodogram vignette: https://ctmm-initiative.github.io/ctmm/articles/periodogram.html

This will also run your through the newer model selection, rather than the older periodogram-based NHT.

Best,

Chris