[R-sig-eco] Brownian Bridge for migratory paths

[Tsewang Namgail]

Hi everyone,

I am using the Brownian Bridge Model to demarcate migratory paths of
birds. I used the following script used for mammals. Using this I
didn't get bridges but simply blotches of high use areas without any
link between them. Can anyone guide me as to how can I get continuous
bridges between all the bird locations. Any help would be highly
appreciated.

Thanks
Tsewang

BrownianBridge=function(X,Y,Time,LocationError,cell.size=50,max.lag)
        {options(digits=20)

#Alternative
#BBMM=brownian.bridge(x=locations$X,y=locations$Y,time.lag=locations$time.lag,location.error=20,cell.size=50)

#X = locsTN$X
#Y = locsTN$Y
#Time = locsTN$Clock
#LocationError = locsTN$TeleError
#cell.size = 500

#max.lag = abs(min(diff(Time))+1)
LocationError = LocationError/sqrt(pi/2)

# Create grid
xmin=(min(X) - 0.1*(max(X)-min(X)))
xmax=(max(X) + 0.1*(max(X)-min(X)))
ymin=(min(Y) - 0.1*(max(Y)-min(Y)))
ymax=(max(Y) + 0.1*(max(Y)-min(Y)))
x = seq(xmin, xmax, by=cell.size)
y = seq(ymin, ymax, by=cell.size)
x = rep(x, length(y))
y = sort(rep(y, length(x)/length(y)))
grid = data.frame(x, y)
cat("Size of grid =", length(unique(grid$y)), "X",
length(unique(grid$x)), fill=TRUE)

	n.locs = length(X)
    cat("Number of animal locations =", n.locs, fill=TRUE)

    #--------------------------------------------------------------
	#Calculate Brownian motion variance "BMvar"
    #LocationError = rep(LocationError, length(X))
    T.jump = alpha = ztz = ob = loc.error.1 = loc.error.2 = NULL
    i = 2
    while(i < n.locs){
    	if((Time[i+1]-Time[i-1])/ 300> max.lag) {i = i + 1}
    	else {
    		ob = c(ob, i)
    		t = Time[i+1]-Time[i-1]
    		T.jump = c(T.jump, t)
    		a = (Time[i]-Time[i-1]) / t
    		alpha = c(alpha, a)
    		u = c(X[i-1], Y[i-1]) + a*(c(X[i+1], Y[i+1]) -
                c(X[i-1], Y[i-1]))
    		ztz = c(ztz, (c(X[i], Y[i]) - u)%*%(c(X[i], Y[i]) - u))
    		loc.error.1 = c(loc.error.1, LocationError[i-1])
    		loc.error.2 = c(loc.error.2, LocationError[i+1])
    	}
    	i = i + 2
    }

    #Likelihood
    likelihood <- function(var){
    	v = T.jump*alpha*(1-alpha)*var +
            ((1-alpha)^2)*(loc.error.1^2) +
    		(alpha^2)*(loc.error.2^2)
    	l = (1/(2*pi*v))*exp(-ztz/(2*v))
    	-sum(log(l), na.rm=T)
    }
    BMvar = nlminb(start=10000, likelihood, lower=10)$par
    cat("Brownian Motion Variance (meters^2) =", round(BMvar),
        fill=TRUE)
    BMvar = rep(BMvar, length(X))

	#Estimate Brownian bridge
	#Note: 5 minute time step (dt interval in Eq.5 Horne et al. 2007) is used.
	Time.Diff = c(diff(Time), NA)
	Max.Time.Diff = max(diff(Time))
	T.Total = Time[n.locs] - Time[1]
	probability = NULL
	int = 0
	for(i in 1:(n.locs-1)){
		#if(Time.Diff[i] <= max.lag){
		if(Time.Diff[i] <= Max.Time.Diff){
			theta = NULL
			tm = 0
			while(tm <= Time.Diff[i]){
				alpha = tm/Time.Diff[i]
				mu.x = X[i] + alpha*(X[i+1] - X[i])
				mu.y = Y[i] + alpha*(Y[i+1] - Y[i])
				sigma.2 = Time.Diff[i]*alpha*(1-alpha)*BMvar[i] +
					((1-alpha)^2)*(LocationError[i]^2) +
					(alpha^2)*(LocationError[i+1]^2)
				ZTZ = (grid$x - mu.x)^2 + (grid$y - mu.y)^2
				theta = (1/sqrt(2*pi*sigma.2))*exp(-ZTZ/(2*sigma.2))
				int = int + theta
				tm = tm + 1800
			}
		}
	}