Skip to content

Commit

Permalink
include tests without OU
Browse files Browse the repository at this point in the history
  • Loading branch information
@NiklasHohmann
NiklasHohmann committed May 22, 2024
1 parent 9705368 commit aca9f9d
Showing 1 changed file with 102 additions and 12 deletions.
114 changes: 102 additions & 12 deletions code/test.paleots.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -3,6 +3,7 @@
#Original script by Melanie J Hopkins 2/26/24.
# Modified by Niklas Hohmann


compile<-array(dim = c(4,100,9))

#### Parameters for simulation ####
Expand All @@ -15,17 +16,20 @@ intra_pop_var = 0.1 # variance of trait values of pupulation around mean trait v
stasis_var = 1 #variance of stasis model
no_of_specimens = 100 # specimens found at each sampling location
grw_step_var = 0.1 # variance for GRW model
set.seed(1) # set seed for reproducibility

#### stasis with test for OU ####

cat("Running stasis simulations \n")
for (j in 1:9){
print(j)
for (i in 1:100){

test<-paleoTS::sim.Stasis(ns=ns[j],
theta = 0,
omega = 0.5,
vp = intra_pop_var,
nn = rep(no_of_specimens, ns[j]))
theta = 0,
omega = 0.5,
vp = intra_pop_var,
nn = rep(no_of_specimens, ns[j]))
fts<-paleoTS::fitSimple(test,model = 'Stasis')
fto<-paleoTS::fitSimple(test,model = 'OU')
ftb<-paleoTS::fitSimple(test,model = 'URW')
Expand All @@ -39,11 +43,11 @@ compile.group<-data.frame(
x=matrix(compile,ncol=1),
y=c(rep(c('Stasis','OU','URW','GRW'),900)),
z=c(rep(5,400),rep(10,400),rep(15,400),rep(20,400),rep(25,400),
rep(35,400),rep(50,400),rep(100,400),rep(200,400)),
rep(35,400),rep(50,400),rep(100,400),rep(200,400)),
stringsAsFactors = FALSE
)
compile.group$y<-ordered(compile.group$y,levels=c('Stasis','OU','URW','GRW'))
jpeg("figs/test_stasis.jpeg")
jpeg("figs/test_stasis_with_ou.jpeg")
boxplot(x~y+z,data = compile.group,
col=c('lightgreen','darkgoldenrod2','dodgerblue3','firebrick3'),
xaxt = 'n',
Expand All @@ -56,17 +60,17 @@ axis(1,at=c(2.5,6.5,10.5,14.5,18.5,22.5,26.5,30.5,34.5),
mtext('green = stasis, yellow = OU, blue = URW, red = GRW',side=3)
dev.off()


#### URW with test for OU ####
cat("Running URW simulations \n")
for (j in 1:9){
print(j)
for (i in 1:100){

test = paleoTS::sim.GRW(ns = ns[j], # time series length
ms = 0, # mean step, set to 0 for URW model
vs = grw_step_var, # step variance
vp = intra_pop_var, # intrapopulation variance
nn = rep(no_of_specimens, ns[j]))
ms = 0, # mean step, set to 0 for URW model
vs = grw_step_var, # step variance
vp = intra_pop_var, # intrapopulation variance
nn = rep(no_of_specimens, ns[j]))
fts<-paleoTS::fitSimple(test,model = 'Stasis')
fto<-paleoTS::fitSimple(test,model = 'OU')
ftb<-paleoTS::fitSimple(test,model = 'URW')
Expand All @@ -84,7 +88,7 @@ compile.group<-data.frame(
)
compile.group$y<-ordered(compile.group$y,levels=c('Stasis','OU','URW','GRW'))

jpeg("figs/test_urw.jpeg")
jpeg("figs/test_urw_with_ou.jpeg")
boxplot(x~y+z,data = compile.group,
col=c('lightgreen','darkgoldenrod2','dodgerblue3','firebrick3'),
xaxt = 'n',
Expand All @@ -96,4 +100,90 @@ axis(1,at=c(2.5,6.5,10.5,14.5,18.5,22.5,26.5,30.5,34.5),
labels = c(5,10,15,20,25,35,50,100,200))
mtext('green = stasis, yellow = OU, blue = URW, red = GRW',side=3)
dev.off()

compile<-array(dim = c(3,100,9))

#### Stasis without test for OU ####

cat("Running stasis simulations \n")
for (j in 1:9){
print(j)
for (i in 1:100){

test<-paleoTS::sim.Stasis(ns=ns[j],
theta = 0,
omega = 0.5,
vp = intra_pop_var,
nn = rep(no_of_specimens, ns[j]))
fts<-paleoTS::fitSimple(test,model = 'Stasis')
ftb<-paleoTS::fitSimple(test,model = 'URW')
ftt<-paleoTS::fitSimple(test,model = 'GRW')
compile[,i,j]<-paleoTS::compareModels(fts,ftb,ftt,silent = TRUE)$modelFits$Akaike.wt
}
}


compile.group<-data.frame(
x=matrix(compile,ncol=1),
y=c(rep(c('Stasis','URW','GRW'),900)),
z=c(rep(5,300),rep(10,300),rep(15,300),rep(20,300),rep(25,300),
rep(35,300),rep(50,300),rep(100,300),rep(200,300)),
stringsAsFactors = FALSE
)
compile.group$y<-ordered(compile.group$y,levels=c('Stasis','URW','GRW'))
jpeg("figs/test_stasis.jpeg")
boxplot(x~y+z,data = compile.group,
col=c('lightgreen','dodgerblue3','firebrick3'),
xaxt = 'n',
xlab = 'Time Series Length',
main = "AICc weight under simulated stasis",
'AICc weight'
)
axis(1,at=c(2.5,6.5,10.5,14.5,18.5,22.5,26.5,30.5,34.5),
labels = c(5,10,15,20,25,35,50,100,200))
mtext('green = stasis, blue = URW, red = GRW',side=3)
dev.off()



#### URW without test for OU ####

cat("Running URW simulations \n")
for (j in 1:9){
print(j)
for (i in 1:100){

test = paleoTS::sim.GRW(ns = ns[j], # time series length
ms = 0, # mean step, set to 0 for URW model
vs = grw_step_var, # step variance
vp = intra_pop_var, # intrapopulation variance
nn = rep(no_of_specimens, ns[j]))
fts<-paleoTS::fitSimple(test,model = 'Stasis')
ftb<-paleoTS::fitSimple(test,model = 'URW')
ftt<-paleoTS::fitSimple(test,model = 'GRW')
compile[,i,j]<-paleoTS::compareModels(fts,ftb,ftt,silent = TRUE)$modelFits$Akaike.wt
}
}

compile.group<-data.frame(
x=matrix(compile,ncol=1),
y=c(rep(c('Stasis','URW','GRW'),900)),
z=c(rep(5,300),rep(10,300),rep(15,300),rep(20,300),rep(25,300),
rep(35,300),rep(50,300),rep(100,300),rep(200,300)),
stringsAsFactors = FALSE
)
compile.group$y<-ordered(compile.group$y,levels=c('Stasis','URW','GRW'))

jpeg("figs/test_urw.jpeg")
boxplot(x~y+z,data = compile.group,
col=c('lightgreen','dodgerblue3','firebrick3'),
xaxt = 'n',
xlab = 'Time Series Length',
main = "AICc weight under simulated GRW",
'AICc weight'
)
axis(1,at=c(2.5,6.5,10.5,14.5,18.5,22.5,26.5,30.5,34.5),
labels = c(5,10,15,20,25,35,50,100,200))
mtext('green = stasis, blue = URW, red = GRW',side=3)
dev.off()
cat("Done! Results are in figs/")

0 comments on commit aca9f9d

Please sign in to comment.