Load all the data

require(smartwolf)
data(smartwolf2)

allwolves <- list(
    Viki = subset(smartwolf, ID == "Viki"),
    Niki = subset(smartwolf, ID == "Niki"),
    Miki = subset(smartwolf, ID == "Miki"),
    LU = subset(smartwolf, ID == "LentuanUros"),
    YN = subset(smartwolf, ID == "YlaneNaaras"),
    H11 = subset(smartwolf, ID2 == "Hessu 2011"),
    H13 = subset(smartwolf, ID2 == "Hessu 2013"),
    Julla = subset(smartwolf, ID == "Julla"),
    Vellu = subset(smartwolf, ID == "Vellu")
)

ids <- names(allwolves)

#try(setwd("c:/users/guru/box sync/WolfResearch/smartwolves/analysis/code"))
#try(setwd("c:/users/elie/box sync/WolfResearch/smartwolves/analysis/code"))

fits <- lapply(ids, function(n){ load(paste0("../results/",n,".fit.robj")); wolf.fit})
names(fits) <- ids
n.zones <- llply(fits, function(fit){ 
  a <- (fit$lambdafit$wolf.fits[[1]]$fit$model %>% row.names)[1:10]
  sum(grepl("1.", a[1:10]))
})

Get the best boundary / predation model for Viki

fits$Viki$estimates

Estimate for Viki, and visualize

wolf.data <- subset(smartwolf, ID == "Viki")
wolf.zones <- createZones(wolf.data, n.zones = 6, plotme=FALSE)
wolf.data$Zone <- wolf.zones$zone.vector
lag.max <- fits$Viki$estimates$lambda[[1]]
tau.max <- fits$Viki$estimates$tau[[1]]

choice.table.best <- buildChoiceTable(wolf.data, lag = lag.max, tau.p = tau.max)
best.fit <-  mlogit(choice ~ 1|(boundary + predation.normalized), data = choice.table.best) 

palette(rich.colors(6))

n <- max(choice.table.best$trip)
z <- with(choice.table.best, zone[choice == 1])


m <- t(best.fit$probabilities)
col.matrix <- m * NA
col.matrix[cbind(as.numeric(z), 1:ncol(m))] <- z

m.cumsum <- rbind(rep(0, n), apply(m, 2, cumsum))
n.zones <- max(as.numeric(z))

barplot(m, space = 0, width = 1, bor = NA, col = alpha(1:6,.3), ylim = c(0,1.1))
for(i in 1:n.zones)
  rect(0:(n-1), m.cumsum[i,], 1:n, m.cumsum[i+1,], 
       col = col.matrix[i,], bor = "white")

axis(1)

matplot(t(m), type = "l", pch= NA, col = 1:6, lty = 1, lwd = 2)
points(1:n, rep(0.5, n), pch = 19, cex = 1, col = z)

Assess score of model

z.obs <- with(choice.table.best, zone[choice == 1]) %>% as.numeric
null.probs <- prop.table(table(z.obs))
(p.null <- mean(null.probs[z.obs]))

## [1] 0.1826221

model.probs <- m[cbind(as.numeric(z), 1:ncol(m))]
mean(model.probs)

## [1] 0.2620451

Model improvement index:

mean(model.probs) / mean(null.probs[z.obs])

## [1] 1.434904

Monte carlo confidence interval of proportion correct:

samples <- 1:1e4
p.correct <- sapply(samples, function(x) mean(rbinom(length(z.obs), size = 1, p = model.probs)))
hist(p.correct/p.null, breaks = 30, col="grey", bor = "darkgrey", prob = TRUE)
abline(v = 1, lwd = 2)

quantile(p.correct/p.null, c(0.025,.5,.975))

##      2.5%       50%     97.5% 
## 0.8589474 1.3957895 2.0400000

Monte carlo confidence interval of improvement-over-null-index

hist(p.correct, breaks = 30, col="grey", bor = "darkgrey", prob = TRUE)
abline(v = p.null, lwd = 2)

quantile(p.correct, c(0.025,.5,.975))

##      2.5%       50%     97.5% 
## 0.1568627 0.2549020 0.3725490

Predicting from model

Load all the data

Get the best boundary / predation model for Viki

Assess score of model