• R语言绘制PCA双标图、碎石图、变量载荷图和变量贡献图

    1、原论文数据双标图

    代码：

    1. setwd("D:/Desktop/0000/R") #更改路径
    2.  
    3. #导入数据
    4. df <- read.table("Input data.csv", header = T, sep = ",")
    5.  
    6. # -----------------------------------
    7. #所需的包:
    8. packages <- c("ggplot2", "tidyr", "dplyr", "readr", "ggrepel", "cowplot", "factoextra")
    9. #安装你尚未安装的R包
    10. installed_packages <- packages %in% rownames(installed.packages())
    11. if (any(installed_packages == FALSE)) {
    12.   install.packages(packages[!installed_packages])
    13. }
    14. invisible(lapply(packages, library, character.only = TRUE))
    15.  
    16. # -----------------------------------
    17. # 设置一些颜色、文字的基础设置
    18. # Colors:
    19. CatCol <- c(
    20.   CSH = "#586158", DBF = "#C46B39", EBF = "#4DD8C0", ENF = "#3885AB", GRA = "#9C4DC4",
    21.   MF = "#C4AA4D", OSH = "#443396", SAV = "#CC99CC", WET = "#88C44D", WSA = "#AB3232"
    22. )
    23. Three_colorblind <- c("#A8AD6F", "#AD6FA8", "#6FA8AD") #c("#809844", "#4f85b0", "#b07495")
    24. graph_elements_dark <- "black"
    25. plot_elements_light <- "gray75"
    26. plot_elements_dark <- "gray25"
    27.  
    28. # Transparency:
    29. boot_alpha_main <- 0.9
    30. boot_alpha_small <- 0.05
    31.  
    32. # Text:
    33. # if (n_pcs > 3) {x_angle <- 270; x_adjust <- 0.25} else {x_angle <- 0; x_adjust <- 0} # option to change orientation of x axis text
    34. x_angle <- 0; x_adjust <- 0
    35. title_text <- 9 # Nature Communications: max 7 pt; cowplot multiplier: 1/1.618; 7 pt : 1/1.618 = x pt : 1; x = 7 / 1/1.618; x = 11.326 (round up to integer)
    36. subtitle_text <- 9
    37. normal_text <- 9 # Nature Communications: min 5 pt; cowplot multiplier: 1/1.618; 5 pt : 1/1.618 = x pt : 1; x = 5 / 1/1.618; x = 8.09 (round up to integer)
    38.  
    39. # Element dimensions:
    40. plot_linewidth <- 0.33
    41. point_shape <- 18
    42. point_size <- 1.5
    43.  
    44. # Initialize figure lists:
    45. p_biplot <- list(); p_r2 <- list(); p_load <- list(); p_contr <- list(); col_ii <- list()
    46.  
    47. # Labels:
    48. veg_sub_labels <- c("All Sites", "All Forests", "Evergreen Needle-Forests") 
    49.  
    50. # -----------------------------------
    51. #选择PCA所需的数据
    52. codes_4_PCA <- c("SITE_ID", "IGBP", "GPPsat", "wLL", "wNmass", "wLMA", "RECOmax") # 选择需要的列数据
    53. #执行筛选
    54. df_subset <- df %>%
    55.   dplyr::select(all_of(codes_4_PCA))
    56. #运行PCA。dplyr::select(-species)：将不需要的列数据去除
    57. pca_result <- FactoMineR::PCA(df_subset %>% dplyr::select(-SITE_ID, -IGBP), scale.unit = T, ncp = 10, graph = F)
    58.  
    59. # -----------------------------------
    60. #绘图
    61. p1<- fviz_pca_biplot(pca_result,
    62.                      axes = c(1, 2),
    63.                      col.ind = df_subset$IGBP, #"grey50",
    64.                      # col.ind = NA, #plot_elements_light, #"white",
    65.                      geom.ind = "point",
    66.                      palette = CatCol,#'futurama',
    67.                      label = "var",
    68.                      col.var = plot_elements_dark,
    69.                      labelsize = 3,
    70.                      repel = TRUE,
    71.                      pointshape = 16,
    72.                      pointsize = 2,
    73.                      alpha.ind = 0.67,
    74.                      arrowsize = 0.5)
    75.  
    76. # -----------------------------------
    77. # 它是ggplot2对象，我们在此基础上进一步修改一下标注。
    78. p1<-p1+
    79.   labs(title = "",
    80.        x = "PC1",
    81.        y = "PC2",
    82.        fill = "IGBP") +
    83.   guides(fill = guide_legend(title = "")) +
    84.   theme(title = element_blank(),
    85.         text = element_text(size = normal_text),
    86.         axis.line = element_blank(),
    87.         axis.ticks = element_blank(),
    88.         axis.title = element_text(size = title_text, face = "bold"),
    89.         axis.text = element_text(size = normal_text),
    90.         #plot.margin = unit(c(0, 0, 0, 0), "cm"),
    91.         # legend.position = "none"
    92.         legend.text = element_text(size = subtitle_text),
    93.         legend.key.height = unit(5, "mm"),
    94.         legend.key.width = unit(2, "mm")
    95.   )
    96. p1

    参考：Leaf-level coordination principles propagate to the ecosystem scale (https://doi.org/10.1038/s41467-023-39572-5)、主成分分析图。

    2、我选用iris数据进行重新绘制测试双标图

    代码：

    1. setwd("D:/Desktop/0000/R") #更改路径
    2.  
    3. #导入数据
    4. df <- read.table("iris1.csv", header = T, sep = ",")
    5.  
    6. # -----------------------------------
    7. #所需的包:
    8. packages <- c("ggplot2", "tidyr", "dplyr", "readr", "ggrepel", "cowplot", "factoextra")
    9. #安装你尚未安装的R包
    10. installed_packages <- packages %in% rownames(installed.packages())
    11. if (any(installed_packages == FALSE)) {
    12.   install.packages(packages[!installed_packages])
    13. }
    14. invisible(lapply(packages, library, character.only = TRUE))
    15.  
    16. # -----------------------------------
    17. # 设置一些颜色、文字的基础设置
    18. # Colors:
    19. CatCol <- c(
    20.   setosa = "#586158", versicolor = "#C46B39", virginica = "#4DD8C0") # 设置类别颜色
    21. Three_colorblind <- c("#A8AD6F", "#AD6FA8", "#6FA8AD") #c("#809844", "#4f85b0", "#b07495")
    22. graph_elements_dark <- "black"
    23. plot_elements_light <- "gray75"
    24. plot_elements_dark <- "gray25"
    25.  
    26. # Transparency:
    27. boot_alpha_main <- 0.9
    28. boot_alpha_small <- 0.05
    29.  
    30. # Text:
    31. # if (n_pcs > 3) {x_angle <- 270; x_adjust <- 0.25} else {x_angle <- 0; x_adjust <- 0} # option to change orientation of x axis text
    32. x_angle <- 0; x_adjust <- 0
    33. title_text <- 9 # Nature Communications: max 7 pt; cowplot multiplier: 1/1.618; 7 pt : 1/1.618 = x pt : 1; x = 7 / 1/1.618; x = 11.326 (round up to integer)
    34. subtitle_text <- 9
    35. normal_text <- 9 # Nature Communications: min 5 pt; cowplot multiplier: 1/1.618; 5 pt : 1/1.618 = x pt : 1; x = 5 / 1/1.618; x = 8.09 (round up to integer)
    36.  
    37. # Element dimensions:
    38. plot_linewidth <- 0.33
    39. point_shape <- 18
    40. point_size <- 1.5
    41.  
    42. # Initialize figure lists:
    43. p_biplot <- list(); p_r2 <- list(); p_load <- list(); p_contr <- list(); col_ii <- list()
    44.  
    45. # Labels:
    46. veg_sub_labels <- c("All Sites", "All Forests", "Evergreen Needle-Forests") 
    47.  
    48. # -----------------------------------
    49. #选择PCA所需的数据
    50. codes_4_PCA <- c("sepal_length", "sepal_width", "petal_length", "petal_width", "species") # 选择需要的列数据
    51. #执行筛选
    52. df_subset <- df %>%
    53.   dplyr::select(all_of(codes_4_PCA))
    54. #运行PCA。dplyr::select(-species)：将不需要的列数据去除
    55. pca_result <- FactoMineR::PCA(df_subset %>% dplyr::select(-species), scale.unit = T, ncp = 10, graph = F)
    56.  
    57. # -----------------------------------
    58. #绘图
    59. p1<- fviz_pca_biplot(pca_result,
    60.                      axes = c(1, 2),
    61.                      col.ind = df_subset$species, #"grey50",
    62.                      # col.ind = NA, #plot_elements_light, #"white",
    63.                      geom.ind = "point",
    64.                      palette = CatCol,#'futurama',
    65.                      label = "var",
    66.                      col.var = plot_elements_dark,
    67.                      labelsize = 3,
    68.                      repel = TRUE,
    69.                      pointshape = 16,
    70.                      pointsize = 2,
    71.                      alpha.ind = 0.67,
    72.                      arrowsize = 0.5)
    73.  
    74. # -----------------------------------
    75. # 它是ggplot2对象，我们在此基础上修改一下标注。
    76. p1<-p1+
    77.   labs(title = "",
    78.        x = "PC1",
    79.        y = "PC2",
    80.        fill = "IGBP") +
    81.   guides(fill = guide_legend(title = "")) +
    82.   theme(title = element_blank(),
    83.         text = element_text(size = normal_text),
    84.         axis.line = element_blank(),
    85.         axis.ticks = element_blank(),
    86.         axis.title = element_text(size = title_text, face = "bold"),
    87.         axis.text = element_text(size = normal_text),
    88.         #plot.margin = unit(c(0, 0, 0, 0), "cm"),
    89.         # legend.position = "none"
    90.         legend.text = element_text(size = subtitle_text),
    91.         legend.key.height = unit(5, "mm"),
    92.         legend.key.width = unit(2, "mm")
    93.   )
    94. p1

    3、iris数据进行绘制碎石图、变量载荷图、变量贡献图

    代码：

    1. #加载包
    2. library(dplyr) #用于数据预处理
    3. library(tidyr) #用于数据预处理
    4. library(stringr) #用于字符串处理
    5. library(modelr) #用于自助法重抽样
    6. library(FactoMineR) #用于PCA
    7. library(ade4) #用于PCA
    8. library(factoextra) #用于PCA结果提取及绘图
    9. #所需的包:
    10. packages <- c("ggplot2", "tidyr", "dplyr", "readr", "ggrepel", "cowplot", "factoextra")
    11. #安装你尚未安装的R包
    12. installed_packages <- packages %in% rownames(installed.packages())
    13. if (any(installed_packages == FALSE)) {
    14.   install.packages(packages[!installed_packages])
    15. }
    16. invisible(lapply(packages, library, character.only = TRUE))
    17.  
    18. setwd("D:/Desktop/0000/R") #更改路径
    19. # 加载数据
    20. df <- read.csv("iris.csv",header = T, row.names = 1) # row.names = 1: 第一列为标签，这时赋值给df时就没有这列了
    21.  
    22. #重抽样
    23. set.seed(123) #设置随机种子
    24. tt=99 #设置重抽样的次数。iris[,-5]：表示去除第5列，因为这是类别
    25. df_boot <- iris[,-5] %>% modelr::bootstrap(n = tt) #重抽样，结果是一个列表，包含499个数据框
    26.  
    27. #使用循环对每一个数据集进行PCA
    28. #初始化3个空变量
    29. N_PCS <- tibble() #使用维数检验保留的PC数量
    30. pca_stats <- tibble() #变量的贡献和载荷
    31. R2 <- c() #解释方差占比
    32.  
    33. #使用循环对每一个数据集进行PCA
    34. #初始化3个空变量
    35. N_PCS <- tibble() #使用维数检验保留的PC数量
    36. pca_stats <- tibble() #变量的贡献和载荷
    37. R2 <- c() #解释方差占比
    38.  
    39. #循环
    40. for (j in 1:tt) {
    41.   ##提取第j次bootstrap的数据
    42.   dat <- df_boot %>% 
    43.     slice(j) %>% # 选择第j行
    44.     pull(strap) %>% # 提取列表
    45.     as.data.frame() # 提取数据集
    46.   
    47.   #使用FactoMineR包执行PCA
    48.   pca_result <- FactoMineR::PCA(dat, scale.unit = T, ncp = 4, graph = F) # ncp = 4：降维几个主成分，设置最大即为全部
    49.   #使用ade4包执行PCA
    50.   # center：指定是否对数据进行中心化，默认为 TRUE。中心化意味着将数据减去各自的均值，使得数据在每个维度上的平均值为零。
    51.   # scale：指定是否对数据进行缩放，默认为 TRUE。缩放意味着将数据除以各自的标准差，使得数据在每个维度上的标准差为一。
    52.   # scannf：指定是否计算特征值和特征向量，默认为 FALSE。如果设置为 TRUE，则会计算特征值和特征向量。
    53.   pca1 <- ade4::dudi.pca(dat, center = TRUE, scale = TRUE, scannf = FALSE, nf = 4) # nf= 4：降维几个主成分，设置最大即为全部
    54.   #检测不确定性和显著性
    55.   #执行维数检验
    56.   pc_tested <-testdim(pca1, nrepet = 999)
    57.   
    58.   ###提取bootstrap数据集的PCA结果
    59.   N_PCS <- N_PCS %>% 
    60.     bind_rows(tibble(strap = j, n_pcs = pc_tested$nb.cor)) #第j次运行的PCA
    61.   
    62.   pca_stats <- bind_rows(pca_stats,
    63.                          pca_result$var$contrib %>% # add contributions
    64.                            as_tibble(rownames = "var") %>%
    65.                            pivot_longer(cols = !var, names_to = "PC", values_to = "contrib") %>% 
    66.                            left_join(pca_result$var$coord %>% # add loadings
    67.                                        as_tibble(rownames = "var") %>%
    68.                                        pivot_longer(cols = !var, names_to = "PC", values_to = "loading"),
    69.                                      by = c("var", "PC")
    70.                            ) %>% 
    71.                            mutate(PC = str_sub(PC, start = 5), #提取PC名称中的数字
    72.                                   strap = j) # bootstrap run number
    73.   ) #得到变量贡献和载荷
    74.   
    75.   R2 <- bind_rows(R2,
    76.                   tibble(PC = pca_result[["eig"]]%>% rownames(),
    77.                          exp_var = pca_result[["eig"]][,2],
    78.                          strap = j) %>% 
    79.                     mutate(PC = str_sub(PC, start = 6)) #提取PC名称中的数字
    80.   )
    81. }
    82.  
    83. #保留的PC数量
    84. N_PCS <- N_PCS %>%
    85.   group_by(n_pcs) %>% 
    86.   summarise(n_rep = n()) %>% #对重复值进行计数
    87.   mutate(retained = n_rep / tt * 100) #计算运行次数百分比
    88. pc_ret <- N_PCS %>% filter(retained == max(retained))
    89. #输出结果的摘要
    90. print(paste0("Number of statistical significant components according to Dray method (Dray et al., 2008) was ",pc_ret[1,1], " in ", round(pc_ret[1,3], digits = 1), "% of runs."))
    91.  
    92.  
    93. n_pcs <- NA #保留PC数的初始设置
    94. # n_pcs <- 2 #可以手动设置保留PC数
    95. if (is.na(n_pcs)) {
    96.   n_pcs <- N_PCS %>% 
    97.     filter(retained == max(retained)) %>% 
    98.     select(n_pcs) %>% 
    99.     unlist() %>% unname()
    100. } #按照Dray等人的方法设置保留PC数
    101.  
    102. ##变量贡献和载荷
    103. pca_stats <- pca_stats %>% 
    104.   group_by(PC, var) %>% 
    105.   mutate(
    106.     contrib_mean = mean(contrib),
    107.     contrib_median = median(contrib),
    108.     contrib_std = sd(contrib),
    109.     # contrib_q25 = quantile(contrib, 0.25), contrib_q75 = quantile(contrib, 0.75),
    110.     loading_mean = mean(loading),
    111.     loading_median = median(loading),
    112.     loading_std = sd(loading),
    113.     # loading_q25 = quantile(loading, 0.25), loading_q75 = quantile(loading, 0.75)
    114.   ) %>% 
    115.   ungroup() %>% 
    116.   dplyr::rename(contrib_boot = contrib, loading_boot = loading) #重命名以免后续的匹配过程出现混乱
    117.  
    118. ##修改PC名称
    119. pca_stats <- pca_stats %>%
    120.   mutate(PC_name = paste0("PC", PC))
    121.  
    122. ##解释方差占比
    123. R2 <- R2 %>% 
    124.   group_by(PC) %>% 
    125.   mutate(
    126.     R2_mean = mean(exp_var),
    127.     R2_median = median(exp_var),
    128.     R2_std = sd(exp_var),
    129.     # R2_q25 = quantile(exp_var, 0.25), R2_q75 = quantile(exp_var, 0.75)
    130.   ) %>% 
    131.   ungroup() %>% 
    132.   dplyr::rename(R2_boot = exp_var) #重命名以免后续的匹配过程出现混乱
    133.  
    134. ##添加到pca_stats的表格中
    135. pca_stats <- pca_stats %>% left_join(R2, by = c("PC", "strap"))
    136.  
    137. #对原始数据的PCA
    138. pca_result <- FactoMineR::PCA(iris[,-5], scale.unit = T, ncp = 4, graph = F)
    139.  
    140. #添加原始数据计算得到的实际值
    141. pca_stats <- pca_stats %>% 
    142.   dplyr::left_join( #添加原始数据的R2（不是bootstrapping的均值）
    143.     tibble(PC = pca_result[["eig"]] %>% rownames(),
    144.            R2 = pca_result[["eig"]][,2]
    145.     ) %>% 
    146.       mutate(PC = str_sub(PC, start = 6)), #提取PC数
    147.     by = "PC"
    148.   ) %>% 
    149.   dplyr::left_join( #添加原始数据的变量贡献（不是bootstrapping的均值）
    150.     pca_result$var$contrib %>% #添加贡献
    151.       as_tibble(rownames = "var") %>%
    152.       pivot_longer(cols = !var, names_to = "PC", values_to = "contrib") %>% 
    153.       mutate(PC = str_sub(PC, start = 5)), #提取PC数
    154.     by = c("PC", "var")
    155.   ) %>%
    156.   dplyr::left_join( #添加原始数据的变量载荷（不是bootstrapping的均值）
    157.     pca_result$var$coord %>% #添加载荷
    158.       as_tibble(rownames = "var") %>%
    159.       pivot_longer(cols = !var, names_to = "PC", values_to = "loading") %>% 
    160.       mutate(PC = str_sub(PC, start = 5)), # extract PC numbers
    161.     by = c("PC", "var")
    162.   )
    163.  
    164. ## 添加PC数的保留百分比（在自助法中PC被保留得有多频繁）
    165. pca_stats <- pca_stats %>% 
    166.   dplyr::left_join(N_PCS %>% dplyr::mutate(PC = n_pcs %>% as.character) %>% dplyr::select(PC, retained),
    167.                    by = "PC"
    168.   ) 
    169.  
    170. # -----------------------------------
    171. # 绘制图碎石图
    172. dat_boot <- pca_stats %>%
    173.   dplyr::select(PC_name, PC, R2_boot) %>% unique()%>% #去除重复
    174.   dplyr::mutate(PC = as.character(PC))
    175.  
    176. dat_true <- pca_stats %>%
    177.   dplyr::select(PC_name, PC, R2, R2_median, R2_std) %>% unique() %>% #去除重复
    178.   dplyr::mutate(PC = as.character(PC))
    179.  
    180. p2 <- ggplot(data = dat_true, aes(x = PC_name, y = R2, group = 1)) + # x = PC -> only numbers on axis, x = PC_name -> can give problems with PC10 being ordered before PC2;
    181.   # group 1 是用来避免某些warning/error的
    182.   geom_errorbar(aes(ymin = R2 - R2_std, ymax = R2 + R2_std),
    183.                 color = Three_colorblind[1], linewidth = plot_linewidth, width = 0.4) + # bootstrapping的标准差
    184.   # geom_bar(stat = "identity", position = position_dodge(), fill = Three_colorblind[1], width = 0.61) + #b07a4f, #9c6a5e, #643c3c
    185.   geom_line(color = Three_colorblind[1]) +
    186.   geom_point(color = Three_colorblind[1], size = point_size) + #实际值
    187.   geom_jitter(data = dat_boot, aes(x = PC_name, y = R2_boot, group = 1), alpha = 0.1,
    188.               color = "black", shape = point_shape, size = 0.5, width = 0.1) + #每次自助样本的值
    189.   geom_point(aes(x = PC_name, y = R2_median), color = plot_elements_dark,
    190.              alpha = boot_alpha_main, shape = point_shape, size = point_size) + #添加自助法得到的中位数值
    191.   geom_text(aes(x = PC_name, y = R2 + R2_std + 2, label = paste0(R2 %>% round(digits = 1), "%")),
    192.             nudge_x = 0.33, size = 2) + #添加数值标注
    193.   labs(title = "", x = "", y = "Explained variance") +
    194.   theme_classic() +
    195.   theme(title = element_blank(),
    196.         text = element_text(size = normal_text),
    197.         axis.line = element_line(color = graph_elements_dark),
    198.         axis.ticks.x = element_line(color = graph_elements_dark),
    199.         axis.ticks.y = element_blank(),
    200.         axis.title = element_text(size = title_text, face = "bold"),
    201.         # axis.title.x = element_blank(), #已经在'labs'中指定
    202.         axis.text = element_text(size = normal_text),
    203.         axis.text.y = element_blank(),
    204.         plot.margin = unit(c(0, 1, 0, 1), "cm"),
    205.         legend.position = "none"
    206.   ) +
    207.   NULL
    208. p2
    209.  
    210. # -----------------------------------
    211. # 绘制变量载荷图
    212. dat_boot <- pca_stats %>%
    213.   dplyr::filter(PC <= n_pcs[1]) %>% #去除额外的PC
    214.   dplyr::select(PC_name, var, loading_boot) %>% unique() #去除重复
    215.  
    216. dat_true <- pca_stats %>%
    217.   dplyr::filter(PC <= n_pcs[1]) %>% #去除额外的PC
    218.   dplyr::select(PC_name, var, loading, loading_median, loading_std) %>% unique() #去除重复
    219.  
    220. p3 <- ggplot(data = dat_true, aes(x = var, y = loading)) +
    221.   facet_grid(. ~ PC_name, scales = "free_y") +
    222.   geom_errorbar(aes(ymin = loading - loading_std, ymax = loading + loading_std), # loading_q25, ymax = loading_q75
    223.                 color = Three_colorblind[2], linewidth = plot_linewidth, width = 0.9) + # standard error = std from bootstrapping
    224.   geom_bar(stat = "identity", position = position_dodge(), fill = Three_colorblind[2]) + #b07a4f, #9c6a5e, #643c3c
    225.   geom_hline(yintercept = 0, color = graph_elements_dark) +
    226.   geom_jitter(data = dat_boot, aes(x = var, y = loading_boot), alpha = boot_alpha_small, color = plot_elements_dark,
    227.               shape = point_shape, size = 0.2, width = 0.1) + #每次自助抽样的值
    228.   geom_point(aes(x = var, y = loading_median), alpha = boot_alpha_main, shape = point_shape,
    229.              size = point_size, color = plot_elements_dark) + #添加自助法得到的中位数值
    230.   coord_flip() +  #对调坐标轴以更好地展示图形
    231.   scale_y_continuous(breaks = waiver(), n.breaks = 4) + #修改x轴（对调后，这就是y轴）
    232.   labs(y = "Loadings", x = "", title = "") +
    233.   theme_classic() +
    234.   theme(title = element_text(size = normal_text, face = "bold"),
    235.         text = element_text(size = normal_text),
    236.         axis.line.x = element_line(color = graph_elements_dark),
    237.         axis.line.y = element_blank(),
    238.         axis.ticks.x = element_line(color = graph_elements_dark),
    239.         axis.ticks.y = element_blank(),
    240.         axis.title = element_text(size = title_text),
    241.         axis.text = element_text(size = normal_text),
    242.         axis.text.x = element_text(angle = x_angle, vjust = x_adjust),
    243.         legend.position = "none",
    244.         legend.title = element_text(size = title_text),
    245.         legend.text = element_text(size = subtitle_text),
    246.         legend.key.height = unit(1.0, "mm"),
    247.         legend.key.width = unit(1.0, "mm"),
    248.         plot.margin = unit(c(0, 0, 0, 0), "cm"),
    249.         strip.text = element_text(face = "bold", size = title_text),
    250.         strip.background = element_blank()
    251.   ) +
    252.   NULL
    253. p3
    254.  
    255. # -----------------------------------
    256. # 绘制变量贡献图
    257. dat_boot <- pca_stats %>%
    258.   dplyr::filter(PC <= n_pcs[1]) %>% #去除额外的PC
    259.   dplyr::select(PC_name, var, contrib_boot) %>% unique() #去除重复
    260.  
    261. dat_true <- pca_stats %>%
    262.   dplyr::filter(PC <= n_pcs[1]) %>% # remove additional PCs
    263.   dplyr::select(PC_name, var, contrib, contrib_median, contrib_std) %>% unique() #去除重复
    264.  
    265. p4<- ggplot(data = dat_true, aes(x = var, y = contrib)) +
    266.   facet_grid(. ~ PC_name, scales = "free_y") +
    267.   geom_errorbar(aes(ymin = contrib_median - contrib_std, ymax = contrib_median + contrib_std), # ymin = contrib_q25, ymax = contrib_q75
    268.                 color = Three_colorblind[3], linewidth = plot_linewidth, width = 0.9) + # standard error = standard deviation from bootstrapping
    269.   geom_bar(stat = "identity", position = position_dodge(), fill = Three_colorblind[3]) + #4f85b0, #59918e, #3c6464
    270.   geom_hline(yintercept = 0, color = graph_elements_dark) +
    271.   geom_jitter(data = dat_boot, aes(x = var, y = contrib_boot), alpha = boot_alpha_small, color = plot_elements_dark,
    272.               shape = point_shape, size = 0.2, width = 0.1) + #每次自助抽样的值
    273.   geom_point(aes(x = var, y = contrib_median), alpha = boot_alpha_main, shape = point_shape,
    274.              size = point_size, color = plot_elements_dark) + #添加自助法得到的中位数值
    275.   coord_flip() + #对调坐标轴以更好地展示图形
    276.   scale_y_continuous(breaks = waiver(), n.breaks = 4) + #添加自助法得到的中位数值
    277.   labs(y = "Contribution [%]", x = "", title = "") +
    278.   theme_classic() +
    279.   theme(title = element_text(size = normal_text, face = "bold"),
    280.         text = element_text(size = normal_text),
    281.         axis.line.x = element_line(color = graph_elements_dark),
    282.         axis.line.y = element_blank(),
    283.         axis.ticks.x = element_line(color = graph_elements_dark),
    284.         axis.ticks.y = element_blank(),
    285.         axis.title = element_text(size = title_text),
    286.         axis.text = element_text(size = normal_text),
    287.         axis.text.x = element_text(angle = x_angle, vjust = x_adjust),
    288.         legend.position = "none",
    289.         legend.title = element_text(size = title_text),
    290.         legend.text = element_text(size = subtitle_text),
    291.         legend.key.height = unit(1.0, "mm"),
    292.         legend.key.width = unit(1.0, "mm"),
    293.         # plot.margin = unit(c(0, 0, 0, 0), "cm"),
    294.         strip.text = element_text(face = "bold", size = title_text),
    295.         strip.background = element_blank()
    296.   ) +
    297.   NULL
    298. p4
    299.  
    300. # -----------------------------------
    301. # 拼图
    302. library(patchwork)
    303. p2+p3/p4

    文中用到的数据代码：R语言绘制 PCA 双标图、碎石图、变量载荷图和变量贡献图（self）.zip - 蓝奏云

  • 相关阅读:
    一次不成功的抓包过程
    136. 只出现一次的数字
    【计算机网络】应用层自定义协议
    CPU和GPU有什么区别？
    我们该如何运营Facebook账号呢？
    条件式注解@ConfigurationOnXXX原理详解
    【C++初阶】前言——C++的发展简述及学习方法分享
    2023年9月21日
    基于Matlab使用地面雷达探测和跟踪LEO卫星星座仿真（附源码）
    golang数组和Slice地址
  • 原文地址：https://blog.csdn.net/qq_45100200/article/details/132941674