/
warning_example2.R
68 lines (56 loc) · 2.1 KB
/
warning_example2.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
## Example using the saddle-node simulation
reps <- 100
cpu <- 2
# Simulate a dataset from the full individual, nonlinear model
T<- 2000
n_pts <- 501
require(stochPop)
pars = c(Xo = 570, e = 0.5, a = 100, K = 1000, h = 200,
i = 0, Da = .03, Dt = 0)
sn <- saddle_node_ibm(pars, times=seq(0,T, length=n_pts))
X <- ts(sn$x1,start=c(sn$time[1], sn$time[2]-sn$time[1]))
png("saddlenode_data.png")
plot(X, cex.lab=2, cex.axis=2, lwd=3, xlab="time", ylab="pop")
dev.off()
# Initialize some parameter estimates as initial guesses for fitting
theta <- mean(X)
sigma <- sd(X)
alpha <- sigma
# Intialize the three models we'll test
ou <- list(alpha=alpha, theta=theta, sigma=sigma)
wa <- list(alpha_0=alpha, theta=theta, sigma=sigma, beta=0)
cpt <- list(alpha_1=alpha, alpha_2=alpha, theta=theta, sigma=sigma, t_shift=time(X)[length(X)]/2)
class(ou) <- c("OU", "list")
class(wa) <- c("warning", "list")
class(cpt) <- c("changePt", "list")
# Fit the models, consider MLE call for OU & warning
ou <- update(ou, X = X)
ou <- update(ou, X = X)
wa <- update(wa, X = X)
wa <- update(wa, X = X)
cpt <- update(cpt, X = X)
cpt <- update(cpt, X = X)
model_list <- list(ou=ou, wa=wa, cpt=cpt)
model_boots <- bootstrapLR(model_list, rep=reps, cpu=cpu)
save(list=ls(), file= "warning_example2.Rdat")
## Plot the results of the model choice
png("model_choice.png", 1600, 400)
par(mfrow=c(2,3))
#par(mfrow=c(2,1))
LRplot(model_boots, test=2,null=1, main= "1 vs 2")
LRplot(model_boots, test=3, null=1, main="1 vs 3")
LRplot(model_boots, 3,2, main = "2 vs 3" )
LRplot(model_boots, 1,2)
LRplot(model_boots, 1,3)
LRplot(model_boots, 2,3)
dev.off()
## And the bootstraps of the parameter estimates
png("ou_boot.png",1600,400);
plot_bootstrap(model_boots, model=1, cex.lab=2, cex.axis=2, lwd=3);
dev.off()
png("wa_boot.png",1600,400); plot_bootstrap(model_boots, model=2, cex.lab=2, cex.axis=2, lwd=3); dev.off()
png("cpt_boot.png",1600,400); plot_bootstrap(model_boots, model=3, cex.lab=2, cex.axis=2, lwd=3); dev.off()
pt <- power_test(ou, wa, nboot=nboot, cpu=cpu, values=c(.001, .o1, .5, seq(1,15, by=3)) )
png("power.png")
plot(pt$values, pt$power)
dev.off()