look at 3rd axis in PCA

R/2-muir-analysis.Rmd

@@ -183,7 +183,7 @@ patchwork <-
   plot_layout(guides = "collect") +
   plot_annotation(tag_levels = "A") &
   theme(panel.grid = element_blank()) &
-  labs(fill = "Precip. (mm)")
+  labs(fill = "Precip.\n(mm/yr)")
 
 patchwork
 ggsave(here("out", "muir-biplots.png"), width = 9, height = 4)

R/5-model-stats.R

@@ -33,6 +33,16 @@ none.summary <-
   full_join(none.pcr.summary) %>% full_join(none.pcr.auto.summary)
 # write_excel_csv(none.summary, here("out", "none summary.csv"))
 
+### how often is each axis significant?
+none.anova <-
+  none.pcr.auto %>%
+  future_map("anova") %>%
+  bind_rows(.id = "dataset") %>% 
+  group_by(term) %>% 
+  summarize(n = n(),
+            n.sig = sum(p.value < 0.05),
+            `%p < 0.05` = n.sig/n*100)
+
 rm(none.pls, none.pcr, none.pls.auto, none.pcr.auto)
 
 
@@ -53,6 +63,16 @@ apparent.summary <-
   full_join(apparent.pcr.summary) %>% full_join(apparent.pcr.auto.summary)
 # write_excel_csv(apparent.summary, here("out", "apparent summary.csv"))
 
+### how often is each axis significant?
+apparent.anova <-
+  apparent.pcr.auto %>%
+  future_map("anova") %>%
+  bind_rows(.id = "dataset") %>% 
+  group_by(term) %>% 
+  summarize(n = n(),
+            n.sig = sum(p.value < 0.05),
+            `%p < 0.05` = n.sig/n*100)
+
 rm(apparent.pls, apparent.pcr, apparent.pls.auto, apparent.pcr.auto)
 
 
@@ -74,6 +94,16 @@ hidden.summary <-
   full_join(hidden.pcr.summary) %>% full_join(hidden.pcr.auto.summary)
 # write_excel_csv(hidden.summary, here("out", "hidden summary.csv"))
 
+### how often is each axis significant?
+hidden.anova <-
+  hidden.pcr.auto %>%
+  future_map("anova") %>%
+  bind_rows(.id = "dataset") %>% 
+  group_by(term) %>% 
+  summarize(n = n(),
+            n.sig = sum(p.value < 0.05),
+            `%p < 0.05` = n.sig/n*100)
+
 rm(hidden.pls, hidden.pcr, hidden.pls.auto, hidden.pcr.auto)
 
 
@@ -117,5 +147,6 @@ outtable <-
 
 write_excel_csv(outtable, here("out", "model statistics table.csv"), na = "-")
 
-
-
+none.anova
+apparent.anova
+hidden.anova

R/6-simulation-figures.Rmd

@@ -405,6 +405,7 @@ hidden.pls.bp <-
 hidden.pls.bp
 ```
 
+
 ## Combine Biplots
 
 I think I want 4 panels, PCA on the top, PLS on the bottom, apparent on the left, hidden on the right (to match figure 3)

R/pcr.R

@@ -7,6 +7,9 @@ library(rsample)
 library(purrr)
 library(here)
 library(lmtest)
+library(broom)
+library(car)
+
 source(here("R", "ropls_helpers.R"))
 
 #' Do PCA logistic regression (PCA-LR) using ropls::opls() and logistic glm.
@@ -91,8 +94,8 @@ pca_lr <- function(data, X_vars, Y_var, reg_pcs = "max", CV = 7, ...){
                  summarize(R2_tjur = diff(means), .groups = "drop_last") %>%
                  as.numeric(),
                p.value = lik.test$`Pr(>Chisq)`[2])
-
-  return(list(pca = pca, scores = scores, glm = m, mod.stats = mod.stats, data = data))
+  marginal <- broom::tidy(car::Anova(m))
+  return(list(pca = pca, scores = scores, glm = m, mod.stats = mod.stats, data = data, anova = marginal))
 }