R语言学习：实际工作用到的R包，如何增进tidyverse技能，GLM模型，tidymodels包学习

# 我实际工作中常用的R包
# 第一种方式：加载工作所需的R包
# 1 数据导入R包
library(RODBC)
library(readr)
library(readxl)
library(xlsx)
library(openxlsx)
library(haven)
# 2 数据整洁包
library(tidyr)
# 3 数据处理包
library(dplyr)
library(data.table)
# 4 数据可视化包
library(ggplot2)
library(ggthemes)
library(patchwork)
# 5 函数实编程包
library(purrr)
# 6 数据建模包
library(caret)
library(scorecard)
# 7 数据报告包
library(rmarkdown)

# 或者
# 第二种方式：加载工作所需R包
library(pacman)
p_load(
  tidyverse,
  RODBC,
  readxl,
  xlsx,
  openxlsx,
  haven,
  data.table,
  ggthemes,
  patchwork,
  caret,
  scorecard,
  rmarkdown
)

library(tidyverse)
library(palmerpenguins)
library(skimr)

data(package = 'palmerpenguins')

# 数据EDA
skim(penguins)
skim(penguins_raw)
penguins %>% 
  slice_head(n = 5)

# 字体设置
# 使用showtext包
library(showtext) 
font_add_google("Lobster Two", "lobstertwo")
font_add_google("Roboto", "roboto")
font_add_google("Poppins", "poppins")
showtext_auto()

# 主题设置
theme_set(theme_bw())
theme_update(
  legend.text = element_text(size=9, family = "roboto"),
  legend.title = element_text(face="bold", size=12, family = "roboto"),
  legend.position = c(1,0),
  legend.justification = c(1, 0),
  text = element_text(family = "lobstertwo", size = 8, color = "black"),
  plot.title = element_text(family = "lobstertwo", size = 20,
                            face = "bold", color="#2a475e"),
  plot.subtitle = element_text(family = "lobstertwo", size = 15, 
                               face = "bold", color="#1b2838"),
  plot.caption = element_text(size = 10),
  plot.title.position = "plot",
  #plot.caption.position = "plot",
  axis.text = element_text(size = 10, color = "black"),
  axis.title = element_text(size=12),
  axis.ticks = element_blank(),
  axis.line = element_line(colour = "grey50"),
  rect = element_blank(),
  panel.grid = element_line(color = "#b4aea9"),
  panel.grid.minor = element_blank(),
  panel.grid.major.x = element_blank(),
  #panel.grid.major.x = element_line(linetype="dashed"),
  #panel.grid.major.y = element_blank(),
  panel.grid.major.y = element_line(linetype="dashed"),
  plot.background = element_rect(fill = '#fbf9f4', color = '#fbf9f4')
)

# 删除含有缺失值的样例
penguins_comp <- penguins %>%
  drop_na()

# 散点图可视化
penguins_comp %>% 
  ggplot(aes(x=flipper_length_mm, y=bill_length_mm)) +
  geom_point(aes(color=species, shape=species), size=2, alpha=0.8) +
  scale_color_manual(values = c("#386cb0","#fdb462","#7fc97f")) +
  labs(
    title = "Palmer Penguins Data Visualization",
    subtitle = "Scatter plot of flipper lenth vs bill length",
    x = "flip length (mm)",
    y = "bill length (mm)"
  )

# 多图合并
# 使用patchwork包
library(patchwork)
p1 <- penguins_comp %>% 
  ggplot(aes(x=flipper_length_mm, y=bill_length_mm)) +
  geom_point(aes(color=species, shape=species), size=2, alpha=0.8) +
  scale_color_manual(values = c("#386cb0","#fdb462","#7fc97f")) +
  labs(x = "flip length (mm)",
       y = "bill length (mm)")

p2 <- penguins_comp %>% 
  ggplot(aes(x=bill_length_mm, y=bill_depth_mm)) +
  geom_point(aes(color=species, shape=species), size=2, alpha=0.8) +
  scale_color_manual(values = c("#386cb0","#fdb462","#7fc97f")) +
  labs(x = "bill length (mm)",
       y = "bill depth (mm)")

p1 + p2 +
  plot_layout(guides = "collect") +
  plot_annotation(
    title = "Palmer Penguins Data Visualization",
    subtitle = "Scatter plots, left: flip length vs bill lengtt; 
    right: bill length vs bill depth") &
  theme(legend.position = "bottom",
        legend.justification = "center")

p3 <- penguins_comp %>% 
  ggplot(aes(x = sex, fill = species)) +
  geom_bar(alpha = 0.8,width=0.6, show.legend = FALSE) +
  scale_fill_manual(values = c("#386cb0","#fdb462","#7fc97f")) +
  facet_wrap(~species, ncol = 3) +
  theme(strip.text = element_text(size=12, face="bold"))

(p1 + p2) / p3 +
  plot_layout(guides = "collect") +
  plot_annotation(
    title = "Palmer Penguins Data Visualization") &
  theme(legend.position = "bottom",
        legend.justification = "center")

# 1 Logistic Regression
# where F is a binary factor and
# x1-x3 are continuous predictors
fit <- glm(F~x1+x2+x3,data=mydata,family=binomial())
summary(fit) # display results
confint(fit) # 95% CI for the coefficients
exp(coef(fit)) # exponentiated coefficients
exp(confint(fit)) # 95% CI for exponentiated coefficients
predict(fit, type="response") # predicted values
residuals(fit, type="deviance") # residuals

#2 Poisson Regression
# where count is a count and
# x1-x3 are continuous predictors
fit <- glm(count ~ x1+x2+x3, data=mydata, family=poisson())
summary(fit) # display results

# tidymodels包的安装和加载
if(!require("tidymodels")) {
  install.packages("tidymodels")
  library("tidymodels")
}

# 第1步：数据划分 rsample包
set.seed(1234)
# 默认75%做训练，25%做测试
split_iris <- initial_split(iris, prop = 2 / 3)
split_iris
str(split_iris)
# 获取训练集和测试集
training_iris <- training(split_iris)
head(training_iris)
testing_iris <-testing(split_iris)
head(testing_iris)

# 第2步：数据集做预处理
first_recipe <- recipe(training_iris) %>%
  recipes::update_role(Sepal.Length, new_role = "outcome")  %>%
  recipes::update_role(Sepal.Width, new_role = "predictor") %>%
  recipes::update_role(Species, new_role = "predictor")
first_recipe
summary(first_recipe)

# 使用step系列函数集做相应的处理
first_recipe <- first_recipe %>%
  step_dummy(Species, one_hot = TRUE)

first_recipe
summary(first_recipe)

prepped_rec <- prep(first_recipe, verbose = TRUE, retain = TRUE )
summary(prepped_rec)
names(prepped_rec)
prepped_rec$var_info
preproc_train <- recipes::bake(prepped_rec, new_data = NULL)
glimpse(preproc_train)
baked_test_pm <- recipes::bake(prepped_rec, new_data = testing_iris)
glimpse(baked_test_pm)

# 第3步 数据建模
# 设计线性模型
Lin_reg_model <- parsnip::linear_reg()
Lin_reg_model <- 
  Lin_reg_model %>%
  parsnip::set_engine("lm") %>%
  parsnip::set_mode("regression")
Lin_reg_model

iris_reg_wflow <-workflows::workflow() %>%
  workflows::add_recipe(first_recipe) %>%
  workflows::add_model(Lin_reg_model)
iris_reg_wflow

iris_reg_wflow_fit <- parsnip::fit(iris_reg_wflow, data = training_iris)
iris_reg_wflow_fit

# 第4步 模型评价
wf_fit <- iris_reg_wflow_fit %>% 
  pull_workflow_fit()

head(wf_fit$fit$fitted.values)
predict(iris_reg_wflow_fit, new_data = training_iris)
wf_fitted_values <- 
  broom::augment(wf_fit$fit, data = preproc_train) %>% 
  select(Sepal.Length, .fitted:.std.resid)

head(wf_fitted_values)

# 计算rmse指标值
yardstick::rmse(wf_fitted_values, 
                truth = Sepal.Length, 
                estimate = .fitted)

# 测试集上面做评估
predict(iris_reg_wflow_fit, new_data = testing_iris)
test_wf_fitted_values <- data.frame(
  true_label = testing_iris$Sepal.Length,
  pred_result = predict(iris_reg_wflow_fit, new_data = testing_iris)$.pred
)
yardstick::rmse(test_wf_fitted_values, 
                truth = true_label, 
                estimate = pred_result)

# 数据可视化
library(patchwork)
# 训练集
p1 <- wf_fitted_values %>%
  ggplot(aes(x = Sepal.Length, y = .fitted)) +
  geom_point() +
  geom_smooth(method = "lm") +
  labs( x = "True Sepal Length", y = "Predicted Sepal Length") +
  theme_classic() +
  ggtitle("训练集可视化")

# 测试集
p2 <- test_wf_fitted_values %>%
  ggplot(aes(x = true_label, y = pred_result)) +
  geom_point() +
  geom_smooth(method = "lm") +
  labs( x = "True Sepal Length", y = "Predicted Sepal Length") +
  theme_classic() +
  ggtitle("测试集可视化")

p1 + p2