.Rhistory

#IC calculator for 0D model. IC=1 gives AICc, while IC=2 gives BIC
############################################################################
icMod0D <- function(N, S, C, loD, lpD, IC){
n0 <- 2*S
nP <- 2*(N-1)*S
ic <- switch(IC,
#AICc
n0*(log(loD) + n0/(n0-2)) + nP*(log(lpD) + nP/(nP-2)) - 2*C,
#BIC
log(N)*2 - 2*(C - (n0/2)*(log(loD) + 1) - (nP/2)*(log(lpD) + 1))
)
return(ic)
}
############################################################################
#IC calculator for UD model. IC=1 gives AICc, while IC=2 gives BIC
############################################################################
icModUD <- function(N, S, C, loD, lpD, loDU, lpDU, IC){
n0 <- 2*S
nP <- 2*(N-1)*S
ic <- switch(IC,
#AICc
n0*(log(loDU)+((n0+2)*(n0-2))/(n0*(n0-4)))+nP*(log(lpDU)+nP/(nP-2))-2*C,
#BIC
log(N)*4-2*(C-(n0/2)*(log(loDU) + loD/loDU)-(nP/2)*(log(lpDU)+lpD/lpDU))
)
return(ic)
}
##################################################
#INDEX ESTIMATION
##################################################
############################################################################
#Calculate the bias of an MCI estimator
############################################################################
biasEta <- function(eta, P, biasP, M, biasM, covPM, varM){
eta*(biasP/P - biasM/M - covPM/(P*M) + varM/M^2)
}
############################################################################
#MCI bias correction and interval estimation based on the Fisher transform
############################################################################
zTrans <- function(eta, etaBias, varEta){
z <- atanh(eta)
dz <- -(1-eta^2)*etaBias
varzdz <- varEta * (1/(1-eta^2) + 2 * eta * etaBias)^2
zCiL <- z + dz - 1.96*sqrt(varzdz)
zCiU <- z + dz + 1.96*sqrt(varzdz)
etaDb <- tanh(z + dz)
etaDbCiL <- tanh(zCiL)
etaDbCiU <- tanh(zCiU)
return(c(etaDb, etaDbCiL, etaDbCiU))
}
############################################################################
#Estimate the bias-corrected set of MCIs with CIs from the 00 model
############################################################################
eta00 <- function(msRes, cntTm=NA){
#Everything is zero in this model
etaDif <- c(0, 0, 0)
etaDft <- c(0, 0, 0)
etaTot <- c(0, 0, 0)
#Package results for return
#return(data.frame(cbind(etaDif, etaDft, etaTot)))
#return(list(data.frame(cbind(etaDif, etaDft, etaTot)), msRes))
if(is.na(cntTm)){
return(list(data.frame(cbind(etaDif, etaDft, etaTot)), msRes))
} else {
return(list(data.frame(cbind(etaDif, etaDft, etaTot)), msRes, cntTm))
}
}
############################################################################
#Estimate the bias-corrected set of MCIs with CIs from the U0 model
############################################################################
etaU0 <- function(N, S, mt, Tt, vU2, sgU, msRes, cntTm=NA){
#Index num and denom components
Pdif <- 0
Pdft <- mt*vU2
M <- 2*sgU + mt*vU2
#Index component covariance Jacobian
cmp.jac <- matrix(c(0, mt, 2, mt), nrow=2, ncol=2, byrow=T)
#Model degrees of freedom
d <- 2*(N*S-1)
#Model component variances
varSgU <- 2 * sgU^2 / d
varVU2 <- 4 * vU2 * sgU / (N * Tt) + 4 * (sgU / (N * Tt))^2
#Model sufficient statistic variances
sfs.var <- matrix(c(varSgU, 0, 0, varVU2), nrow=2, ncol=2, byrow=T)
#Index estimator component covariances by applying the Jacobian to the suff   stat vars
cmp.cov <- cmp.jac %*% sfs.var %*% t(cmp.jac)
#Biases of index estimator components
biasPdif <- 0
biasPdft <- 2*(mt/(N*Tt))*sgU
biasM <- 2*(mt/(N*Tt))*sgU
#Biased index estimators
etaDifBiased <- 0
etaDftBiased <- Pdft/M
etaTotBiased <- etaDftBiased
#Jacobian for transforming index estimator components into index covariances
ind.jac <- (1/M)*matrix(c(1, -etaDftBiased, 1, -etaTotBiased), nrow=2, ncol
=2, byrow=T)
#Index covariances
ind.cov <- ind.jac %*% cmp.cov %*% t(ind.jac)
#etaDif estimator bias
biasEtaDif <- 0
#etaDft estimator bias
biasEtaDft <- biasEta(etaDftBiased, Pdft, biasPdft, M, biasM, cmp.cov[1, 2],
cmp.cov[2,2])
#etaTot estimator bias
biasEtaTot <- biasEtaDft
#Debiased index estimates plus 95% CI end points
etaDif <- c(0, 0, 0)
etaDft <- zTrans(etaDftBiased, biasEtaDft, ind.cov[1,1])
etaTot <- zTrans(etaTotBiased, biasEtaTot, ind.cov[2,2])
#Package results for return
if(is.na(cntTm)){
return(list(data.frame(cbind(etaDif, etaDft, etaTot)), msRes))
} else {
return(list(data.frame(cbind(etaDif, etaDft, etaTot)), msRes, cntTm))
}
}
############################################################################
#Estimate the bias-corrected set of MCIs with CIs from the 0D model
############################################################################
eta0D <- function(N, S, loD, lpD, msRes, cntTm=NA){
#Index num and denom components
Pdif <- (2/N)*(loD - lpD)
Pdft <- 0
M <- (2/N)*(loD + (N - 1)*lpD)
#Index component covariance Jacobian
cmp.jac <- matrix(c(2/N, -2/N, 2/N, 2*(N-1)/N), nrow=2, ncol=2, byrow=T)
#Model degrees of freedom
no <- 2*S
np <- 2*(N-1)*S
#Model component variances
varLoD <- 2 * loD^2 / no
varLpD <- 2 * lpD^2 / np
#Model sufficient statistic variances
sfs.var <- matrix(c(varLoD, 0, 0, varLpD), nrow=2, ncol=2, byrow=T)
#Index estimator component covariances by applying the Jacobian to the        sufficient statistic vars
cmp.cov <- cmp.jac %*% sfs.var %*% t(cmp.jac)
#Biases of index estimator components
biasPdif <- 0
biasPdft <- 0
biasM <- 0
#Biased index estimators
etaDifBiased <- Pdif/M
etaDftBiased <- 0
etaTotBiased <- etaDifBiased
#Jacobian for transforming index estimator components into index covariances
ind.jac <- (1/M)*matrix(c(1, -etaDifBiased, 1, -etaTotBiased), nrow=2, ncol
=2, byrow=T)
#Index covariances
ind.cov <- ind.jac %*% cmp.cov %*% t(ind.jac)
#etaDif estimator bias
biasEtaDif <- biasEta(etaDifBiased, Pdif, biasPdif, M, biasM, cmp.cov[1, 2],
cmp.cov[2,2])
#etaDft estimator bias
biasEtaDft <- 0
#etaTot estimator bias for DU model
biasEtaTot <- biasEtaDif
#Debiased index ests plus 95% CI end points
etaDif <- zTrans(etaDifBiased, biasEtaDif, ind.cov[1,1])
etaDft <- c(0, 0, 0)
etaTot <- zTrans(etaTotBiased, biasEtaTot, ind.cov[2,2])
#Package results for return
#return(data.frame(cbind(etaDif, etaDft, etaTot)))
#return(list(data.frame(cbind(etaDif, etaDft, etaTot)), msRes))
if(is.na(cntTm)){
return(list(data.frame(cbind(etaDif, etaDft, etaTot)), msRes))
} else {
return(list(data.frame(cbind(etaDif, etaDft, etaTot)), msRes, cntTm))
}
}
############################################################################
#Estimate the bias-corrected set of MCIs with CIs from the UD model
############################################################################
etaUD <- function(N, S, mt, Tt, vU2, loDU, lpDU, msRes, cntTm=NA){
#Index num and denom components
Pdif <- (2/N)*(loDU - lpDU)
Pdft <- mt*vU2
M <- (2/N)*(loDU + (N - 1)*lpDU) + mt*vU2
#Index component covariance Jacobian
cmp.jac <- matrix(c(2/N, -2/N, 0, 0, 0, mt, 2/N, 2*(N-1)/N, mt), nrow=3,
ncol=3, byrow=T)
#Model degrees of freedom
no <- 2*S
np <- 2*(N-1)*S
#Model component variances
varLoDU <- 2 * loDU^2 / (no - 2)
varLpDU <- 2 * lpDU^2 / np
varVU2 <- 4 * vU2 * loDU / (N * Tt) + 4 * (loDU / (N * Tt))^2
#Sufficient statistic variances
sfs.var <- matrix(c(varLoDU, 0, 0, 0, varLpDU, 0, 0, 0, varVU2), nrow=3,
ncol=3, byrow=T)
#Index estimator component covariances by applying the Jacobian to the suff   stat vars
cmp.cov <- cmp.jac %*% sfs.var %*% t(cmp.jac)
#Biases of index estimator components
biasPdif <- 0
biasPdft <- 2*(mt/(N*Tt))*loDU
biasM <- 2*(mt/(N*Tt))*loDU
#Biased index estimators for DU model
etaDifBiased <- Pdif/M
etaDftBiased <- Pdft/M
etaTotBiased <- etaDifBiased + etaDftBiased
#Jacobian for transforming index estimator components into index covariances
ind.jac <- (1/M)*matrix(c(1, 0, -etaDifBiased, 0, 1, -etaDftBiased, 1, 1,
-etaTotBiased), nrow=3, ncol=3, byrow=T)
#Index covariances
ind.cov <- ind.jac %*% cmp.cov %*% t(ind.jac)
#etaDif estimator bias
biasEtaDif <- biasEta(etaDifBiased, Pdif, biasPdif, M, biasM, cmp.cov[1, 3],
cmp.cov[3,3])
#etaDft estimator bias
biasEtaDft <- biasEta(etaDftBiased, Pdft, biasPdft, M, biasM, cmp.cov[2, 3],
cmp.cov[3,3])
#etaTot estimator bias for DU model
biasEtaTot <- biasEta(etaDftBiased, Pdif + Pdft, biasPdif + biasPdft, M,
biasM, cmp.cov[1, 3] + cmp.cov[2, 3], cmp.cov[3,3])
#Debiased index estimates and 95% CI end points
etaDif <- zTrans(etaDifBiased, biasEtaDif, ind.cov[1,1])
etaDft <- zTrans(etaDftBiased, biasEtaDft, ind.cov[2,2])
etaTot <- zTrans(etaTotBiased, biasEtaTot, ind.cov[3,3])
#Package results for return
if(is.na(cntTm)){
return(list(data.frame(cbind(etaDif, etaDft, etaTot)), msRes))
} else {
return(list(data.frame(cbind(etaDif, etaDft, etaTot)), msRes, cntTm))
}
}
############################################################################
#A function that returns an information criterion for associate with
#a given piece of data. Performs model selection on the given piece
#of data, and gives the result from the lowest IC model
#Can choose between AICc (IC=1) or BIC (IC=2)
############################################################################
getIC <- function(data, IC=1){
#Eliminate any row that has an NA anywhere in it
data <- data[complete.cases(data), ]
#Check that W is still respected after removing NAs
#NOTE: THIS FUNCTION SHOULD BE REVISED TO EXPLICITLY ACCEPT W
#AS AN ARGUMENT
if(dim(data)[1] >= 2){
#Get the dimensions of the no-NA dataset
nCol <- ncol(data)
N <- (nCol-1)/2
S <- nrow(data) - 1
#Times are expected to be in POSIXct format
#if(class(data[,1])[1]!="POSIXct"){
# data[,1] <- as.POSIXct(strptime(data[, 1], dateFormat))
#}
tms <- data[ ,1]
ti <- diff(tms)
units(ti) <- "secs"
ti <- as.numeric(ti)
#Get the regular timestep from the data
mt <- median(ti)
#Calculate total T
Tt <- sum(ti)
#Calculate the position of the center of the current window
#winCnt <- i - (W - 1)/2
#td1 <- diff(c(t0, tmsAll[floor(winCnt)]))
#td2 <- diff(c(t0, tmsAll[ceiling(winCnt)]))
#td <- (td1 + td2)/2
#cntTm <- t0 + td
#cnt[cntWins] <- cntTm
#Get the X and Y positions for all individuals
x <- data[, 2:(N+1)]
y <- data[, (N + 2):(2*N + 1)]
#Calculate the x and y position shifts for all individuals
x1 <- x[1:S,]
x2 <- x[2:(S + 1),]
dx <- x2 - x1
y1 <- y[1:S,]
y2 <- y[2:(S + 1),]
dy <- y2 - y1
#Mean vector is just 0's for the no drift models
v0 <- rep(0, N)
#Estimate the X and Y uniform drift values
vxU <- (1/(N*Tt)) * c(sum( dx ))
vyU <- (1/(N*Tt)) * c(sum( dy ))
#Calculate the squared drift
vU2 <- vxU^2 + vyU^2
#Create vectors of the X and Y uniform drift values of length N
vyU <- rep(vyU, N)
vxU <- rep(vxU, N)
#Estimate the sufficient statistics s and r for the no drift models
s0 <- (1/2)*(s(N, S, dx, v0, ti) + s(N, S, dy, v0, ti))
r0 <- (1/2)*(r(N, S, dx, v0, ti) + r(N, S, dy, v0, ti))
#Estimate the sufficient statistics s and r for the uniform drift              models
sU <- (1/2)*(s(N, S, dx, vxU, ti) + s(N, S, dy, vyU, ti))
rU <- (1/2)*(r(N, S, dx, vxU, ti) + r(N, S, dy, vyU, ti))
##################################################
#Model-specific parameter estimates
#Lambda ests for correlated diffusion, uniform drift
loDU <- N*S*rU / (S - 1)
lpDU <- N*(sU - rU) / (N - 1)
#Variance est for uncorrelated diffusion, uniform drift
sgU <- (N*S / (N*S-1)) * sU
#Lambda ests for correlated diffusion, no drift
loD <- N*r0
lpD <- N*(s0 - r0) / (N - 1)
#Variance est for uncorrelated diffusion, no drift
sg <- s0
##################################################
##################################################
#Model selection
#Model independent constant
C <- -N*sum(log(2*pi*ti))
#No drift, uncorrelated diffusion
ic00 <- icMod00(N, S, C, sg, IC)
#Uniform drift, uncorrelated diffusion
icU0 <- icModU0(N, S, C, sg, sgU, IC)
#No drift, correlated diffusion
ic0D <- icMod0D(N, S, C, loD, lpD, IC)
#Uniform drift, correlated diffusion
icUD <- icModUD(N, S, C, loD, lpD, loDU, lpDU, IC)
#Identify the selected model
icVals <- c(ic00, icU0, ic0D, icUD)
#aicVals <- c(aic00, aicU0)
msRes <- which.min(icVals)
#msRes <- 4
return(icVals[msRes])
} else {
#If W is violated after removing NAs, return NA
return(NA)
}
}
############################################################################
#A function to find the single best partition within a chunk of data
#based on IC differences. Returns NA if no such partition exists.
############################################################################
icPart <- function(data, W=2, IC=1){
#Get the length of the dataset passed to the function
nPart <- nrow(data)
#Check if it is possible to introduce a partition without violating
#the window constraint. If not, return NA
if(nPart >= (2*W-1)){
pCnt <- 0
partIC <- NULL
#Loop over all possible partition points
for(j in W:(nPart-W+1)){
pCnt <- pCnt + 1
#Check if there is any missing data at the focal partition point.
#If so, skip it and return NA, if not, evaluate it.
if(sum(!is.na(data[j, ]))==length(data[j, ])){
#Split the data into two pieces at the focal partition point j
part1 <- data[1:j, ]
part2 <- data[j:nPart, ]
#Get the IC vals for each piece and add them together
ic1 <- getIC(part1, IC)
ic2 <- getIC(part2, IC)
partIC[pCnt] <- ic1 + ic2
} else {
partIC[pCnt] <- NA
}
} #end j for
#Find the lowest IC value associated with the candidate partition points.
sPartLoc <- which.min(partIC)
minIC <- partIC[sPartLoc]
#Check if the lowest IC partition point is better than leaving the        #data chunk intact.
if((getIC(data, IC) - minIC) > 0){
#If the best partition point beats the intact data chunk, return its
#abosulte location and associated IC val.
loc <- (W-1) + sPartLoc
return(data.frame(loc=loc, ic=minIC))
} else {
#If it is not possible to beat the IC of the intact data chunk,
#return NAs.
return(data.frame(loc=NA, ic=NA))
} #end minAIC if
} else {
#If it is not possible to introduce a partition point without
#violating W, return NAs.
return(data.frame(loc=NA, ic=NA))
} #end nPart if
}
############################################################################
#A function to find partitions in the dataset via IC differences
#and a window-based stopping rule
############################################################################
findPrts <- function(data, W, IC=1){
#Set the absolute minimum window width to the level that, below which
#models can no longer be fit due to insufficient data
MW <- 2
#Get the length of the dataset
dNd <- dim(data)[1]
#Start with full dataset IC, then modify to reflect changes via partitioning
dIcNew <- getIC(data, IC)
dIcOld <- dIcNew + 1
#Give the start and end of the dataset as the initial "partition points"
pPS <- c(1, dNd)
while(dIcNew < dIcOld){
#Update the looping criterion
dIcOld <- dIcNew
cIC <- NULL
cPS <- NULL
cMW <- NULL
#Loop over partitions
for(i in 1:(length(pPS)-1)){
#Get the start and end points of the current partitions
pSt <- pPS[i]
pNd <- pPS[i+1]
#Try to partition the current chunk of data
prt <- icPart(data[pSt:pNd, ], MW, IC)
#If the attempted partitioning worked (ie, did not produce NAs)
#Update the IC diffs, partition location lists, and min width of
#the resulting new partition.
#If the partitioning failed, return all NAs.
if(!is.na(prt$ic)){
icFC <- getIC(data[pSt:pNd, ], IC)
cIC <- c(cIC, icFC - prt$ic)
cps <- pSt + prt$loc - 1
cPS <- c(cPS, cps)
cMW <- c(cMW, min(c(pNd-cps+1, prt$loc)))
} else {
cIC <- c(cIC, NA)
cPS <- c(cPS, NA)
cMW <- c(cMW, NA)
}
}
#Find the biggest AIC difference
#PROBABLY NEED TO DEAL WITH CASE OF MULTIPLE IDENTICAL MAX DELTA IC VALS     #HERE
mIC <- which.max(cIC)
#If this position if not an NA, and the resulting partitions do not
#violate the min window size, then accept the new partition and
#update the partition locations list, and also the IC value for
#the whole dataset.
if(length(mIC) == 1){
if(!is.na(cMW[mIC]) & (cMW[mIC] >= W)){
pPS <- sort(c(pPS, cPS[mIC]))
dIcNew <- dIcOld - cIC[mIC]
} else {
dIcNew <- dIcOld
}
} else {
dIcNew <- dIcOld
}
}
#Return the list of partition points
return(pPS)
}
ptm <- proc.time()
nRep <- 8
Si <- 50
Ns <- c(5, 10)
#Ns <- c(5, 10, 15, 20, 25)
sig1 <- 2
sig2 <- 6
#rho0 <- 0
rho1 <- 0.5
rho2 <- 0.25
mt <- 86400
for(N in Ns){
#print(N)
#Create mean vectors
#muX0 <- rep(0, N)
#muY0 <- rep(0, N)
#PS1
muX1 <- rep(2, N)
muY1 <- rep(2, N)
muX2 <- rep(5, N)
muY2 <- rep(-3, N)
#cov0 <- covMat(rho0, sig0, N)
cov1 <- covMat(rho1, sig1, N)
cov2 <- covMat(rho2, sig2, N)
x <- foreach(j=1:nRep) %dopar% {  #.combine='rbind'
xdata1 <- sim(N, Si, muX1, muX2, cov1, cov2, "X", miss=0)
ydata1 <- sim(N, Si, muY1, muY2, cov1, cov2, "Y")
xdata2 <- sim(N, Si, muX1, muX2, cov1, cov2, "X", miss=0)
ydata2 <- sim(N, Si, muY1, muY2, cov1, cov2, "Y")
xdata3 <- sim(N, Si, muX1, muX2, cov1, cov2, "X", miss=0)
ydata3 <- sim(N, Si, muY1, muY2, cov1, cov2, "Y")
xdata4 <- sim(N, Si, muX1, muX2, cov1, cov2, "X", miss=0)
ydata4 <- sim(N, Si, muY1, muY2, cov1, cov2, "Y")
xdata5 <- sim(N, Si, muX1, muX2, cov1, cov2, "X", miss=0.25)
ydata5 <- sim(N, Si, muY1, muY2, cov1, cov2, "Y")
xdata6 <- sim(N, Si, muX1, muX2, cov1, cov2, "X", miss=0.25)
ydata6 <- sim(N, Si, muY1, muY2, cov1, cov2, "Y")
xdata7 <- sim(N, Si, muX1, muX2, cov1, cov2, "X", miss=0.25)
ydata7 <- sim(N, Si, muY1, muY2, cov1, cov2, "Y")
xdata8 <- sim(N, Si, muX1, muX2, cov1, cov2, "X", miss=0.25)
ydata8 <- sim(N, Si, muY1, muY2, cov1, cov2, "Y")
ts <- tsp(mt, 7*Si)
#Create the combined dataset in the same format as real data
dataAll1 <- data.frame(timestamp=ts, xdata1, ydata1)
dataAll2 <- data.frame(timestamp=ts, xdata2, ydata2)
dataAll3 <- data.frame(timestamp=ts, xdata3, ydata3)
dataAll4 <- data.frame(timestamp=ts, xdata4, ydata4)
dataAll5 <- data.frame(timestamp=ts, xdata5, ydata5)
dataAll6 <- data.frame(timestamp=ts, xdata6, ydata6)
dataAll7 <- data.frame(timestamp=ts, xdata7, ydata7)
dataAll8 <- data.frame(timestamp=ts, xdata8, ydata8)
dataSim1 <- as.corrData(dataAll1, timeformat = "%m/%d/%y")
dataSim2 <- as.corrData(dataAll2, timeformat = "%m/%d/%y")
dataSim3 <- as.corrData(dataAll3, timeformat = "%m/%d/%y")
dataSim4 <- as.corrData(dataAll4, timeformat = "%m/%d/%y")
dataSim5 <- as.corrData(dataAll5, timeformat = "%m/%d/%y")
dataSim6 <- as.corrData(dataAll6, timeformat = "%m/%d/%y")
dataSim7 <- as.corrData(dataAll7, timeformat = "%m/%d/%y")
dataSim8 <- as.corrData(dataAll8, timeformat = "%m/%d/%y")
res1 <- findPrts(dataSim1, W=15)
res2 <- findPrts(dataSim2, W=25)
res3 <- findPrts(dataSim3, W=35)
res4 <- findPrts(dataSim4, W=45)
res5 <- findPrts(dataSim5, W=15)
res6 <- findPrts(dataSim6, W=25)
res7 <- findPrts(dataSim7, W=35)
res8 <- findPrts(dataSim8, W=45)
resOut <- list(res1, res2, res3, res4, res5, res6, res7, res8)
resOut
}
filename <- paste("simPS2_", "N", toString(N), ".Rdata", sep="")
save(x, file = filename)
}
proc.time() - ptm