Net squared displacement (NSD) looks at the movement vectors of animals to determine their use of the landscape (Bunnefeld et al. 2011, Papworth et al. 2012). Bunnefeld et al. (2011) determined a novel method to ideintfy movement patters using NSD which is the straight line distance between an animals' srating location and subsequent locations. Movements were then categorized into one of 5 categories based on the top model that describes movement for each individual. In the code for this section, we have updated the code to include the 5 movement equations along with R code provided in Papworth et al. (2012) to enable determination of migratory, mixed migratory, disperser, home range, or nomadic movement behavior. Figure 1
in Bunnefeld et al. (2011) is a helpful to interpret the output of this code that identifies the pattern of NSD over time that is determined by which of the 5 movement behaviors the animal follows. Those interested in this section should read the 2 papers cited for more details and specifics of the methods.

  1. Exercise 3.6 - Download and extract zip folder into your preferred location
  2. Set working directory to the extracted folder in R under File - Change dir...
  3. First we need to load the packages needed for the exercise
    library(trip)
    library(stringr)
    library(adehabitatHR)
    library(lattice)
    library(gmodels)
    library(spatstat)
    library(maptools)

  4. Now open the script "NSDscript.R" and run code directly from the script
    muleys<-read.csv("DCmuleysedited.csv", header=T, sep=",")
    str(muleys)
    #Remove outlier locations, if needed, by first plotting coords below then editing
    #this line accordingly
    newmuleys <-subset(muleys, muleys$Long > -110.50 & muleys$Lat > 37.3
    & muleys$Long < -107)
    muleys <- newmuleys
    #Make a spatial data frame of locations after removing outliers
    coords<-data.frame(x = muleys$Long, y = muleys$Lat)
    crs<-"+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"
    head(coords)
    plot(coords)
    deer.spdf <- SpatialPointsDataFrame(coords= coords, data = muleys,
    proj4string = CRS(crs))
    head(deer.spdf)
    class(deer.spdf)
    proj4string(deer.spdf)
    plot(deer.spdf,col=deer.spdf$id)
    muleys$NewDate<-as.POSIXct(muleys$GPSFixTime, format="%Y.%m.%d %H:%M:%S", origin="1970-01-01")
    #TIME DIFF NECESSARY IN BBMM CODE
    timediff <- diff(muleys$NewDate)*60*60
    # remove first entry without any difference
    muleys <- muleys[-1,]
    muleys$timelag <-as.numeric(abs(timediff))
    summary(muleys$timelag)
    #Remove locations greater than 24 hours apart in time
    muleys <- subset(muleys, muleys$timelag < 18000)

  5. The key to NSD is proper delineation of movement periods so we can explore a few alternatives here. First we will define deer and year based simply on the Year the location was recorded. Not very biologically meaningful but for simplicity we will start with calender year. Be sure to skip step 6 below and continue on with step 7 to the end of the exercise.
    muleys$Year <- format(muleys$NewDate, "%Y")
    muleys <- subset(muleys, muleys$Year != "NA")
    muleys$YearBurst <- c(paste(muleys$id,muleys$Year,sep="_"))
    muleys$YearBurst <- as.factor(muleys$YearBurst)
    str(muleys)
    summary(muleys$YearBurst)
    muleys <- subset(muleys, table(muleys$YearBurst)[muleys$YearBurst] > 100)
    muleys$YearBurst <- factor(muleys$YearBurst)

  6. Alternatively, we could assign a Year as when we would expect mule deer to disperse or migrate from summer to winter range. In our Colorado example, the mule deer were captured on 30 September 2011 (winter range) so we will start there and separate out a summer season from May to August using the code below. Be sure to use skip step 5 above and run step 6 instead but do not run both.
    range(muleys$NewDate)
    muleys$Year <- NULL
    muleys$Year[muleys$NewDate > "2011-09-30 00:30:00" & muleys$NewDate
    < "2012-04-01 23:00:00"] <- 2011
    muleys$Year[muleys$NewDate > "2012-03-31 00:30:00" & muleys$NewDate
    < "2012-10-01 23:00:00"] <- 2012
    muleys$Year <- as.factor(muleys$Year)
    muleys$YearBurst <- c(paste(muleys$id,muleys$Year,sep="_"))
    muleys$YearBurst <- as.factor(muleys$YearBurst)
    table(muleys$YearBurst)]

  7. Next we are going to rename our dataset to simply follow along with previous code
    unless you want to change d1 to muleys throughout
    d1 <- muleys
    #Code separate each animal into a shapefile or text file to use as a "List"
    # get input file
    indata <- d1
    innames <- unique(d1$YearBurst)
    innames <- innames[1:12]#needs to be number of unique IDs
    outnames <- innames
    # begin loop to calculate home ranges
    for (i in 1:length(innames)){
    data <- indata[which(indata$YearBurst==innames[i]),]
    if(dim(data)[1] != 0){
    #data <-data[c(-21)]
    # export the point data into a shp file
    data.xy = data[c("X", "Y")]
    coordinates(data.xy) <- ~X+Y
    sppt <- SpatialPointsDataFrame(coordinates(data.xy),data)
    proj4string(sppt) <- CRS("+proj=utm +zone=12 +datum=WGS84")
    #writePointsShape(sppt,fn=paste(outnames[i],sep="/"),factor2char=TRUE)
    #sppt <-data[c(-22,-23)]
    write.table(sppt, paste(outnames[i],"txt",sep="."), sep="\t", quote=FALSE,
    row.names=FALSE)
    write.table(paste(outnames[i],"txt",sep="."), sep="\t", quote=FALSE,
    row.names=FALSE, col.names=FALSE, "In_list.txt", append=TRUE)
    #The write.table line above should only be run once to create the In_list.txt
    #file otherwise it rights all animals each time
    #Note: There are 3 lines of code that are not active that can be activated
    #to export point shapefiles for all resulting animals but need to remove
    #as.POSIX column first
    }}
    ##############################################################################
    ##############################################################################
    ##
    #Code below to get NSD and best movement model for each bird
    ##
    ##############################################################################
    ##############################################################################
    date()
    # Reads the List file of GPS datasets
    List<-read.table("In_list.txt",sep="\t",header=F)
    head(List) #List contains the filenames of the deer datasets
    # Generation of results vectors
    ID <- rep(0,nrow(List))
    LOCS <- rep(0,nrow(List))
    MIGR <- rep(0,nrow(List))
    MIXM <- rep(0,nrow(List))
    DISP <- rep(0,nrow(List))
    HORA <- rep(0,nrow(List))
    NOMA <- rep(0,nrow(List))
    ID <- rep(0,nrow(List))
    LOCS <- rep(0,nrow(List))

    MIGR <- rep(0,nrow(List))
    MIXM <- rep(0,nrow(List))
    DISP <- rep(0,nrow(List))
    HORA <- rep(0,nrow(List))
    NOMA <- rep(0,nrow(List))
    AICC_1 <- rep(0,nrow(List))
    AICC_2 <- rep(0,nrow(List))
    AICC_3 <- rep(0,nrow(List))
    AICC_4 <- rep(0,nrow(List))
    AICC_5 <- rep(0,nrow(List))
    minAIC <- rep(0,nrow(List))
    d_AICC_1 <- rep(0,nrow(List))
    d_AICC_2 <- rep(0,nrow(List))
    d_AICC_3 <- rep(0,nrow(List))
    d_AICC_4 <- rep(0,nrow(List))
    d_AICC_5 <- rep(0,nrow(List))
    LL_AICC_1 <- rep(0,nrow(List))
    LL_AICC_2 <- rep(0,nrow(List))
    LL_AICC_3 <- rep(0,nrow(List))
    LL_AICC_4 <- rep(0,nrow(List))
    LL_AICC_5 <- rep(0,nrow(List))
    sumLL_AICC <- rep(0,nrow(List))
    wi_AICC_1 <- rep(0,nrow(List))
    wi_AICC_2 <- rep(0,nrow(List))
    wi_AICC_3 <- rep(0,nrow(List))
    wi_AICC_4 <- rep(0,nrow(List))
    wi_AICC_5 <- rep(0,nrow(List))
    for(i in 1:nrow(List)) {
    coords<-read.table(as.character(List[i,]),sep="\t",header=T)
    coords$DT<-as.POSIXct(coords$NewDate, format="%Y-%m-%d %H:%M:%S")
    ##make a data.frame of coordinates. Here the raw values are divided
    #by 1000 so that trajectories are calculated using km as the unit of
    #measurement not meters
    coord<-data.frame((coords$Y/1000),(coords$X/1000))
    # make ltraj: a trajectory of all the relocations
    d2<-as.ltraj(coord,coords$DT, coords$YearBurst, #separate your data by individual.
    burst=coords$YearBurst, #burst is used to creat subdivisions within an individual.
    typeII=TRUE) #typeII can be TRUE: radio-track data, or FALSE: not time
    #recorded, such as tracks in the snow
    #[1] "2007-06-05 16:00:00 EDT" "2015-03-12 14:00:00 EDT"
    #you can now make your trajectory regular
    #first create a reference start time
    #refda <- strptime("00:00:00", "%H:%M:%S")#all relocations should be altered
    #to occur at 30 seconds past each minute

    #you can now make your trajectory regular, as radio tracks tend to lose
    #a few seconds / minutes with each relocation
    #firstly add "NA" for each missing location in your trajectory
    #d3<-setNA(d2,refda,
    #as.POSIXct("2007-06-01 06:00:00 EDT"), #any time before earliest timedate
    #7200, #stating there should be a location every 2 hours
    #tol=7200, #how many time units to search each side of expected location
    #units="sec") #specifying the time units
    #you can now make your trajectory regular
    #firstly create a reference start time
    refda <- strptime("00:00:30", "%H:%M:%S")
    ##NOTE: The refda and d3 code above was not run because it results in too many
    #relocations as "NA" that get removed below. Not quite sure the reason behind
    # it being included? We can now make your trajectory regular
    d4<-sett0(d2, refda,
    10800, #stating the interval at which relocations should be
    correction.xy =c("none"), #if "cs" performs location correction based on the
    #assumption the individual moves at a constant speed
    tol=10800, #how many time units to search either side of an expected location
    units = "sec") #specifying the time units
    #to view your regular trajectory of points with NA's
    summary(d4)
    #now calculating NSD for each point
    datansd<-NULL
    for(n in 1:length(summary(d4)[,1])) #stating that NSD should be
    #calculated separately for each burst
    {
    nsdall<-d4[[n]][,8] #extracting the NSD for each location
    nsdtimeall<-d4[[n]][,3] #extracting the time for each location
    nsdtimestartzero<-d4[[n]][,3]-d4[[n]][1,3]
    #extracting the time since trip start for each location
    nsdid<-rep(as.vector(summary(d4)[n,1]),
    length.out=summary(d4)[n,3])
    #extracting the individual associated with each location
    nsdtrip<-rep(as.vector(summary(d4)[n,2]),length.out=summary(d4)[n,3])
    #extracting the trip associated with each location
    datansd1<-data.frame(nsdall,nsdtimeall,nsdtimestartzero,nsdid,nsdtrip)
    #joining all these variables together in a data frame
    datansd<-rbind(datansd,datansd1)
    #joining all the data frames together
    }
    datansd$zero1<-as.numeric(unclass(datansd$nsdtimestartzero))
    # making seconds since trip start numeric
    datansd$zerostart<-datansd$zero1/60
    #changing the time since trip start from seconds to minutes
    datansd$minslitr2<-as.numeric(strftime(as.POSIXlt(datansd$nsdtimeall),
    format="%M"))
    #making a vector of the hour of the day a location occured

    datansd$hdaylitr2<-as.numeric(strftime(as.POSIXlt(datansd$nsdtimeall),
    format="%H"))
    #making a vector of the minute in an hour a location occured
    datansd$minsday<-((datansd$hdaylitr2*60)+datansd$minslitr2)
    #calculating the minute in the day a location occured
    summary(datansd)
    datansd1<-na.omit(datansd) #remove NA's
    datansd1$coordinates<-coord #add the coordinates for each point
    #you now have the dataframe you need (datansd) to start analysis
    #NSD
    #table(datansd1$nsdid)
    #Now we can start modelling NSD using nlme.
    #Equations are from Bunnefeld at al (2011) A model-driven approach to quantify
    #migration patterns: #individual, regional and yearly differences.
    #Journal of Animal Ecology 80: 466 - 476
    #First we are going to model the data using nls, a least squares method,
    #the simplest method and first method in Bunnefeld et al. 2011 (i.e., MIGRATION)
    #that uses a double sigmoid or s-shaped function.
    ###########################
    ##
    ## MIGRATION
    ##
    ###########################
    m1<-tryCatch(nls(nsdall ~ asym /(1+exp((xmidA-zerostart)/scale1)) +
    (-asym / (1 + exp((xmidB-zerostart)/scale2))), #Equation 1 in Bunnefeld et al. 2011
    start = c(asym=15000000,xmidA=200000,xmidB=450000,scale1=1000,scale2=1000)
    #these are the starting values for each parameter of the equation
    ,data=na.omit(datansd1)),error=function(e)99999) #this is the data
    summary(m1) #this will print a summary of the converged model
    #NOTE: The error function is simply to prevent the loop from crashing if model
    #does not converge
    ###########################
    ##
    ## MIXED MIGRATORY
    ##
    ###########################
    m2 <-tryCatch(nls(nsdall ~ asymA /(1+exp((xmidA-zerostart)/scale1)) +
    (-asymB / (1 + exp((xmidB-zerostart)/scale2))), #Equation 2 in Bunnefeld et al. 2011
    start = c(asymA=15000000,asymB=10000000, xmidA=200000,xmidB=450000,scale1=1000,
    scale2=1000)
    #these are the starting values for each parameter of the equation

    ,data=na.omit(datansd1)),error=function(e)99999) #this is the data
    summary(m2) #this will print a summary of the converged model
    ###########################
    ##
    ## DISPERSAL
    ##
    ###########################
    m3 <-tryCatch(nls(nsdall ~ asym /(1+exp((xmid-zerostart)/scale)), #Equation 3 in
    #Bunnefeld et al. 2011
    start = c(asym=15000000,xmid=200000,scale=1000)
    #these are the starting values for each parameter of the equation
    ,data=na.omit(datansd1)),error=function(e)99999) #this is the data
    summary(m3) #this will print a summary of the converged model
    ###########################
    ##
    ## HOME RANGE
    ##
    ###########################
    m4 <- tryCatch(nls(nsdall ~ intercept, data=na.omit(datansd1),
    start = list(intercept = 0)),error=function(e)99999)
    #Equation 4 in Bunnefeld et al. 2011 where c is a constant
    summary(m4) #this will print a summary of the converged model
    ###########################
    ##
    ## NOMADIC
    ##
    ###########################
    m5 <- tryCatch(nls(nsdall ~ beta*zerostart,start=c(beta=1),
    data=na.omit(datansd1)),error=function(e)99999) #Equation 5 in Bunnefeld et al. 2011
    #where beta is a constant and t the number of days since initial start date
    # (i.e., 1 June of each year)
    summary(m5) #this will print a summary of the converged model
    #Below we are going to set up the AIC table
    ID[i] <- paste(unique(as.factor(datansd$nsdid)))
    LOCS[i] <- nrow(coords)
    MIGR[i] <- print(tryCatch(AIC(m1),error=function(e)0))
    MIXM[i] <- print(tryCatch(AIC(m2),error=function(e)0))
    DISP[i] <- print(tryCatch(AIC(m3),error=function(e)0))
    HORA[i] <- print(tryCatch(AIC(m4),error=function(e)0))
    NOMA[i] <- print(tryCatch(AIC(m5),error=function(e)0))
    AICC_1[i] <- print(tryCatch(AIC(m1),error=function(e)99999))

    AICC_3[i] <- print(tryCatch(AIC(m3),error=function(e)99999))
    AICC_4[i] <- print(tryCatch(AIC(m4),error=function(e)99999))
    AICC_5[i] <- print(tryCatch(AIC(m5),error=function(e)99999))
    minAIC[i] <- min(AICC_1[i],AICC_2[i],AICC_3[i],AICC_4[i],AICC_5[i])
    d_AICC_1[i] <- (AICC_1[i] - minAIC[i])
    d_AICC_2[i] <- (AICC_2[i] - minAIC[i])
    d_AICC_3[i] <- (AICC_3[i] - minAIC[i])
    d_AICC_4[i] <- (AICC_4[i] - minAIC[i])
    d_AICC_5[i] <- (AICC_5[i] - minAIC[i])
    LL_AICC_1[i] <- exp(-0.5*d_AICC_1[i])
    LL_AICC_2[i] <- exp(-0.5*d_AICC_2[i])
    LL_AICC_3[i] <- exp(-0.5*d_AICC_3[i])
    LL_AICC_4[i] <- exp(-0.5*d_AICC_4[i])
    LL_AICC_5[i] <- exp(-0.5*d_AICC_5[i])
    sumLL_AICC[i] <- sum(LL_AICC_1[i],LL_AICC_2[i],LL_AICC_3[i],LL_AICC_4[i],LL_AICC_5[i])
    wi_AICC_1[i] <- LL_AICC_1[i]/sumLL_AICC[i]
    wi_AICC_2[i] <- LL_AICC_2[i]/sumLL_AICC[i]
    wi_AICC_3[i] <- LL_AICC_3[i]/sumLL_AICC[i]
    wi_AICC_4[i] <- LL_AICC_4[i]/sumLL_AICC[i]
    wi_AICC_5[i] <- LL_AICC_5[i]/sumLL_AICC[i]
    filename<-paste(substr(List[i,],1,8),"png",sep=".")
    #NOTE:Numbers after "List[i,] need to encompass possible lengths of output name
    #(i.e., D19.txt is 6 characters)
    png(filename,height=20,width=30,units="cm",res=600)
    #graphical exploration of the data will help you find sensible starting values
    #for each of the parameters asym, xmidA, xmidB, scale1 and scale2.
    #to graph nsd against time, use:
    #xyplot(nsdall~zerostart|nsdtrip,data=datansd)
    #str(nsdtest)
    #now plot the data with the predicted curve
    nsdplot <- xyplot(nsdall ~ zerostart/3600, data=datansd1,
    col="grey", #color for the observed locations
    type='b', # 'b' shows the locations as dots, with a line connecting
    #successive locations. Can also be 'p' for just the locations, or 'l' for just
    #the line between locations
    ylab=expression(paste('Net squared displacement ',' ', (km^2))), #y axis label
    xlab="Hours after trip start")
    plot(nsdplot)
    dev.off()
    }
    #Create table of AIC values with lower AIC identifying best model

    RESULT<-cbind(ID,LOCS,MIGR,MIXM,DISP,HORA,NOMA)
    colnames(RESULT)<- c("ID","LOCS","MIGR","MIXM","DISP","HORA","NOMA")#
    RESULT
    #write.table(RESULT,"OUT_AIC_NSD.txt",sep="\t")#Output to excel if needed
    #Create table of raw values to calculate AICweights
    Migratory <- rbind(ID,AICC_1,d_AICC_1,LL_AICC_1,wi_AICC_1)
    MixedMig <- rbind(ID,AICC_2,d_AICC_2,LL_AICC_2,wi_AICC_2)
    Disperser <- rbind(ID,AICC_3,d_AICC_3,LL_AICC_3,wi_AICC_3)
    HomeRange <- rbind(ID,AICC_4,d_AICC_4,LL_AICC_4,wi_AICC_5)
    Nomadic <- rbind(ID,AICC_5,d_AICC_5,LL_AICC_5,wi_AICC_5)
    RESULT2 <- rbind(Migratory,MixedMig,Disperser,HomeRange,Nomadic)
    RESULT2