-
Notifications
You must be signed in to change notification settings - Fork 3
/
Copy pathctmm_akde.R
147 lines (105 loc) · 3.34 KB
/
ctmm_akde.R
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
##########
# AKDE
##########
library(ctmm)
help("akde")
help("bandwidth")
# load buffalo data
data(buffalo)
projection(buffalo) <- median(buffalo)
names(buffalo)
# here we will work with Pepper
DATA <- buffalo$Pepper
COL <- color(DATA,by="time")
plot(DATA,col=COL,main="Pepper")
# this dataset has problems
dt.plot(DATA)
# selected autocorrelation model
GUESS <- ctmm.guess(DATA,interactive=FALSE)
FIT <- ctmm.select(DATA,GUESS,trace=3)
# save(FIT,file="pepper.rda")
# I've already run this
load("data/pepper.rda")
summary(FIT)
# analogous IID model
IID <- ctmm.fit(DATA)
summary(IID)
# regular KDE
KDE <- akde(DATA,IID)
# default AKDE
AKDE <- akde(DATA,FIT)
# optimally weighted AKDE
wAKDE <- akde(DATA,FIT,weights=TRUE)
# you only need this with irregular sampling - can be slow
# unweighted AKDE places too much density on oversampled times
# Pepper's optimal weights
plot(DATA$timestamp,wAKDE$weights,xlab="time",ylab="weight",main="Optimal Weights")
plot(DATA$timestamp,wAKDE$weights,xlab="time",ylab="weight",main="Optimal Weights",ylim=c(0,0.005))
# matching extent for plotting
EXT <- extent(list(KDE,AKDE,wAKDE))
plot(DATA,KDE,ext=EXT,main="KDE")
# note CIs, grid, etc...
summary(KDE)
plot(DATA,AKDE,ext=EXT,main="AKDE")
summary(AKDE)
plot(DATA,wAKDE,ext=EXT,main="optimally weighted AKDE")
summary(wAKDE)
# Over-smoothing bias
osAKDE <- akde(DATA,FIT,weights=TRUE,debias=FALSE)
plot(DATA,osAKDE,ext=EXT,main="uncorrected wAKDE")
###########################
# Home-range meta-analysis
###########################
help("meta")
FITS <- list()
for(i in 1:length(buffalo))
{
GUESS <- ctmm.guess(buffalo[[i]],interactive=FALSE)
FITS[[i]] <- ctmm.select(buffalo[[i]],GUESS,trace=3)
}
names(FITS) <- names(buffalo)
# save(FITS,file="data/buffalo.rda")
load("data/buffalo.rda")
# calculate AKDES on a consistent grid
AKDES <- akde(buffalo,FITS,weights=TRUE)
# save(AKDES,file="data/buffalo_akdes.rda")
load("data/buffalo_akdes.rda")
# color to be spatially distinct
COL <- color(AKDES,by='individual')
# plot AKDEs
plot(AKDES,col.UD=COL,col.level=COL,col.grid=NA,level=NA,main="African buffalo AKDEs")
# Mean buffalo HR "the old way"
AREA <- vector("numeric", length = length(AKDES))
for(i in 1:length(AKDES))
{ AREA[i] <- summary(AKDES[[i]],units=FALSE)$CI[2] }
AREA
mean(AREA) # mean
sqrt(var(AREA)/length(AREA)) # SE
help('meta',package="ctmm")
# meta-analysis of buffalo home-range areas
meta(AKDES,col=c(COL,'black'),sort=TRUE)
# model selection: Dirac-delta > inverse-Gaussian
# force inverse-Gaussian population distribution
meta(AKDES,plot=FALSE,IC=NA)
# since CoV isn't a selected feature, its underestimated here
# comparing sub-groups
BUFFALO <- list(South=AKDES[1:3],North=AKDES[4:6])
META <- meta(BUFFALO)
META
META['South/','/North',]
# more general meta-analytic regressions
help("Log")
# then you can use the 'metafor' R package
#########################
# Population density
#########################
# this is a straight mean of the individual densities that doesn't model population variance
help("mean.UD")
# note the 'sample' argument for correct CIs
# straight mean - for a population of 6 buffalo
MEAN <- mean(AKDES,sample=FALSE)
plot(buffalo,MEAN,col=COL,main="Mean African buffalo AKDE")
# this is a population kernel density estimate (paper coming)
help("pkde")
PKDE <- pkde(buffalo,AKDES)
plot(buffalo,PKDE,col=COL,main="African buffalo PKDE")