Merge pull request #1 from MindTheGap-ERC/remove_OU

Remove ou

README.md

@@ -6,7 +6,9 @@ Examine paleoTS model selection performance with time series length
 
 ## Introduction
 
-Hohmann et al 2023 ([DOI: 10.1101/2023.12.18.572098](https://doi.org/10.1101/2023.12.18.572098)), supplementary code ([DOI: 10.5281/zenodo.10390267](https://doi.org/10.5281/zenodo.10390267)) found unexpected behavior in the model selection performance of the paleoTS package, version 0.5.3 ([CRAN.R-project.org/package=paleoTS](https://CRAN.R-project.org/package=paleoTS)). Here, this behavior is examined using the simulation tools implemented in the paleoTS package.
+Hohmann et al 2023 ([DOI: 10.1101/2023.12.18.572098](https://doi.org/10.1101/2023.12.18.572098)), supplementary code ([DOI: 10.5281/zenodo.10390267](https://doi.org/10.5281/zenodo.10390267)) found unexpected behavior in the model selection performance of the paleoTS package, version 0.5.3 ([CRAN.R-project.org/package=paleoTS](https://CRAN.R-project.org/package=paleoTS)) when including Ornstein-Uhlenbeck processes in the set of modes of evolution tested for. Here, this behavior is examined using the simulation tools implemented in the paleoTS package.
+
+The code simulates stasis and undirected random walks, and examines how AICc weights for different modes of evolution (stasis, (un)directed random walk) change with time series length.
 
 ## Authors
 
@@ -52,8 +54,10 @@ Results of the analysis are already stored in the repository under _figs/_. You
 * code : folder for code
   * test.paleots.R : R script examining paleoTS model selection with time series length
 * figs : folder for figures
-  * test_stasis.jpeg : AICc weights under stasis model with increasing time series length
-  * test_urw.jpeg : AICc weights under undirected random walk (URW) model with increasing time series length
+  * test_stasis_with_ou.jpeg : AICc weights under stasis model with increasing time series length, including OU in the tested modes
+  * test_stasis_without_ou.jpeg : AICc weights under stasis model with increasing time series length, not including OU in the tested modes
+  * test_urw_with_ou.jpeg : AICc weights under undirected random walk (URW) model with increasing time series length, including OU in the tested modes
+  * test_urw_without_ou.jpeg : AICc weights under undirected random walk (URW) model with increasing time series length, not including OU in the tested modes
 * renv : folder for `renv` package
 * .Rprofile : R session info
 * .gitignore : untracked files

code/test.paleots.R

@@ -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 ####
@@ -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")
+cat("Running stasis simulations 1/2 \n")
 for (j in 1:9){
-  print(j)
+  cat(paste0("time series length ", ns[j],"\n"))
   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')
@@ -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',
@@ -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()
 
-
-cat("Running URW simulations \n")
+#### URW with test for OU #### 
+cat("Running URW simulations 1/2 \n")
 for (j in 1:9){
-  print(j)
+  cat(paste0("time series length ", ns[j],"\n"))
   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')
@@ -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',
@@ -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  2/2 \n")
+for (j in 1:9){
+  cat(paste0("time series length ", ns[j],"\n"))
+  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_without_ou.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 2/2 \n")
+for (j in 1:9){
+  cat(paste0("time series length ", ns[j],"\n"))
+  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_without_ou.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/")