自动诊断模型、出图表和解释结果且可视化的神奇R包-easystats

“

easystats是个生态，最有用的是report + parameter + performance！精力有限，更强大的功能大家自己去开发吧。

# install.packages("easystats")
library("easystats")
#更新使用easystats::install_latest()
#install.packages("deepdep")
library(deepdep)
plot_dependencies("easystats", depth = 2, show_stamp = FALSE)
#----??1.reposrt包----
#---1.1 汇报data----
report(iris)

#分组数据
library(tidyverse)
iris %>%
  select(-starts_with("Sepal")) %>%
  group_by(Species) %>%
  report() %>%
  summary()


#---1.2 汇报模型，可文字可table----
library(report)
model <- lm(Sepal.Length ~ Species, data = iris)
report(model)

#--自动格式化测试t检验 ，结果为文字！
report(t.test(mtcars$mpg ~ mtcars$am))
#--自动格式化测试t检验 ，整理为表格！
report(t.test(mtcars$mpg ~ mtcars$am))%>%
  as.data.frame()
#--测试相关系数检验，结果输出表格
cor.test(iris$Sepal.Length, iris$Sepal.Width) %>%
  report() %>%
  as.data.frame()
#--测试F分析
aov(Sepal.Length ~ Species, data = iris) %>%
  report()%>%
  as.data.frame()
###----glm-logistic----
model <- glm(vs ~ mpg * drat, data = mtcars, family = "binomial")
#文字版
report(model)
#表格版
report(model) %>% as.data.frame()
#对于复杂的报告,可以直接查阅报告的各个部分:
report_model(model)
report_performance(model)
report_statistics(model)

#--混合效应模型
library(lme4)
model <- lme4::lmer(Sepal.Length ~ Petal.Length + (1 | Species), data = iris)
#文字版
report(model)
#表格版
report(model) %>% as.data.frame()
#--bayes模型
library(report)
library(rstanarm)
model <- stan_glm(mpg ~ qsec + wt, data = mtcars)
report(model)

#报告编写工具版本
report(sessionInfo())

#----??2.correlation包----
library("correlation")
library("dplyr")
results <- correlation(iris)
results
summary(results)
#--可视化
library(see)
results %>%
  summary(redundant = TRUE) %>%
  plot()
plot(cor_test(iris, "Sepal.Width", "Sepal.Length"))

##----2.1 stratified correlations----
iris %>%
  select(Species, Sepal.Length, Sepal.Width, Petal.Width) %>%
  group_by(Species) %>%
  correlation()
##----2.2 Bayesian Correlations----
correlation(iris, bayesian = TRUE)

correlation(iris, include_factors = TRUE, method = "auto")
#偏相关
iris %>%
  correlation(partial = TRUE) %>%
  summary()
##----2.3 Gaussian Graphical Models ----
#高斯图形模型
library(see) # for plotting
library(ggraph) # needs to be loaded
plot(correlation(mtcars, partial = TRUE)) +
  scale_edge_color_continuous(low = "#000004FF", high = "orange")
##---2.4 Multilevel Correlations----
iris %>%
  correlation(partial = F, multilevel = TRUE) %>%
  summary()


#----3.datawizard包，不够强----
#install.packages("datawizard")
library(datawizard)
library(dplyr)
#过滤条件的行vs = 0, am = 1
data_match(mtcars, data.frame(vs = 0, am = 1))
#或者
data_filter(mtcars, vs == 0 & am == 1)
#查找列或检索选定列
data_select(iris, starts_with("Sepal"))  %>%  head()
data_extract(mtcars, "gear")

#----4.effectsize ----
#install.packages("effectsize")
library(effectsize)
options(es.use_symbols = TRUE)
##----4.1 标准化差异----
cohens_d(mpg ~ am, data = mtcars)
hedges_g(mpg ~ am, data = mtcars)
glass_delta(mpg ~ am, data = mtcars)
# Dependence
phi(mtcars$am, mtcars$vs)
cramers_v(mtcars$am, mtcars$cyl)
fei(table(mtcars$cyl), p = c(0.1, 0.3, 0.6))
#eta_squared
model <- aov(mpg ~ factor(gear), data = mtcars)
eta_squared(model)
omega_squared(model)
epsilon_squared(model)

d_to_r(d = 0.2)
#OR与RR转换
oddsratio_to_riskratio(2.6, p0 = 0.4)
#
F_to_d(15, df = 1, df_error = 60)
F_to_r(15, df = 1, df_error = 60)
F_to_eta2(15, df = 1, df_error = 60)

##---4.2 指标解释----
interpret_r(r = 0.3)
#不同rule解释不一样
interpret_cohens_d(d = 0.45, rules = "cohen1988")
interpret_cohens_d(d = 0.45, rules = "gignac2016")

#----5.insight ----
#install.packages("insight")
library(insight)
m <- lm(Sepal.Length ~ Species + Petal.Width + Sepal.Width, data = iris )
#get_data
dat <- get_data(m)
pred <- find_predictors(m, flatten = TRUE)
model_info(m)

l <- lapply(pred, function(x) {
  if (is.numeric(dat[[x]])) {
    mean(dat[[x]])
  } else {
    unique(dat[[x]])
  }
})
names(l) <- pred
l <- as.data.frame(l)
cbind(l, predictions = predict(m, newdata = l))

print_params <- function(model) {
  paste0(
    "My parameters are ",
    toString(row.names(summary(model)$coefficients)),
    ", thank you for your attention!"
  )
}

m1 <- lm(Sepal.Length ~ Petal.Width, data = iris)
print_params(m1)
m2 <- mgcv::gam(Sepal.Length ~ Petal.Width + s(Petal.Length), data = iris)
print_params(m2)

#----6.modelbased ----
#install.packages("modelbased")
##----6.1 智能网格来表示复杂的交互作用----
library(ggplot2)
library(see)
library(modelbased)
# Fit model
model <- lm(Sepal.Length ~ Petal.Length * Petal.Width, data = iris)
#visualisation matrix with expected Z-score values
vizdata <- modelbased::visualisation_matrix(model, by = c("Petal.Length", "Petal.Width = c(-1, 0, 1)"))
# 将预期的标准差值还原为实际值
vizdata <- unstandardize(vizdata, select = "Petal.Width")
#添加模型的预测关系
vizdata <- modelbased::estimate_expectation(vizdata)
#将Petal.Width表示为z得分（“-1 SD”、“+2 SD”等
vizdata$Petal.Width <- effectsize::format_standardize(vizdata$Petal.Width, reference = iris$Petal.Width)
#绘图
ggplot(iris, aes(x = Petal.Length, y = Sepal.Length)) +
  # 添加来自原始数据集的点
  geom_point2(aes(fill = Petal.Width), shape = 21, size = 5) +
  # 添加关系线
  geom_line(data = vizdata, aes(y = Predicted, color = Petal.Width), linewidth = 1) +
  # 改进颜色/主题
  scale_color_viridis_d(direction = -1) +
  scale_fill_viridis_c(guide = "none") +
  theme_modern()

##----6.2 估计边际均值，marginal means----
model <- lm(Sepal.Width ~ Species, data = iris)
means <- estimate_means(model)
ggplot(iris, aes(x = Species, y = Sepal.Width)) +
  # Add base data
  geom_violin(aes(fill = Species), color = "white") +
  geom_jitter2(width = 0.05, alpha = 0.5) +
  # Add pointrange and line from means
  geom_line(data = means, aes(y = Mean, group = 1), linewidth = 1) +
  geom_pointrange(
    data = means,
    aes(y = Mean, ymin = CI_low, ymax = CI_high),
    size = 1,
    color = "white"
  ) +
  # Improve colors
  scale_fill_material() +
  theme_modern()
model <- lm(Sepal.Width ~ Species, data = iris)

##----6.3 对比分析----
library(ggplot2)
library(see)
library(modelbased)
model <- lm(Sepal.Width ~ Species, data = iris)
#边际均值
means <- estimate_means(model)
ggplot(means, aes(x = Species, y = Mean)) +
  geom_line(aes(group = 1)) +
  geom_pointrange(aes(color = Species, ymin = CI_low, ymax = CI_high)) +
  theme_modern()
estimate_contrasts(model)

#复杂一些
model <- lm(Sepal.Width ~ Species * Petal.Width, data = iris)
contrasts <- estimate_contrasts(model)
contrasts
##----6.3a 对比分析灯塔图方法-----
library(see)
plot(contrasts, estimate_means(model)) +
  theme_modern()

##----6.4 差异的改变-----
##有问题！
#model <- lm(Sepal.Width ~ Species * Petal.Length, data = iris)
#estimate_contrasts(model, by = "Petal.Length", length = 3)
#contrasts <- estimate_contrasts(model, by = "Petal.Length", length = 20)
#contrasts$Contrast <- paste(contrasts$Level1, "-", contrasts$Level2)
## Plot
#ggplot(contrasts, aes(x = Petal.Length, y = Difference, )) +
#  # Add line and CI band
#  geom_line(aes(color = Contrast)) +
#  geom_ribbon(aes(ymin = CI_low, ymax = CI_high, fill = Contrast), alpha = 0.2) +
#  # Add line at 0, indicating no difference
#  geom_hline(yintercept = 0, linetype = "dashed") +
#  # Colors
#  theme_modern()

 
##----6.5 pre与原始概率比较-----
library(ggplot2)
model1 <- lm(Petal.Length ~ Sepal.Length, data = iris)
pred1 <- estimate_expectation(model1)
model2 <- lm(Petal.Length ~ Sepal.Length * Species, data = iris)
pred2 <- estimate_expectation(model2)
pred2$Petal.Length <- iris$Petal.Length
# Initialize plot for model 1
ggplot(data = pred1, aes(x = pred2$Petal.Length, y = Predicted)) +
  # with identity line (diagonal) representing perfect predictions
  geom_abline(linetype = "dashed") +
  # Add the actual predicted points of the models
  geom_point(aes(color = "Model 1")) +
  geom_point(data = pred2, aes(color = "Model 2")) +
  # Aesthetics changes
  labs(y = "Petal.Length (predicted)", color = NULL) +
  scale_color_manual(values = c("Model 1" = "blue", "Model 2" = "red")) +
  theme_modern()

#----7.performance----
#install.packages("performance")
library(performance)
#计算许多不同模型的 r2，lm glm polr bayes lme4
##----7.1 r2----
model <- lm(mpg ~ wt + cyl, data = mtcars)
r2(model)

model <- glm(am ~ wt + cyl, data = mtcars, family = binomial)
r2(model)

library(MASS)
data(housing)
model <- polr(Sat ~ Infl + Type + Cont, weights = Freq, data = housing)
r2(model)

set.seed(123)
library(rstanarm)
model <- stan_glmer(
  Petal.Length ~ Petal.Width + (1 | Species),
  data = iris,
  cores = 4
)
r2(model)

library(lme4)
model <- lmer(Reaction ~ Days + (1 + Days | Subject), data = sleepstudy)
r2(model)

##----7.2 ICC----
#ICC 提供有关解释方差的信息，可以解释为总体中分组解释的方差比例
library(lme4)
model <- lmer(Reaction ~ Days + (1 + Days | Subject), data = sleepstudy)
icc(model)

#方法2 卡
#library(brms)
#set.seed(123)
#model <- brm(mpg ~ wt + (1 | cyl) + (1 + wt | gear), data = mtcars)
#icc(model)


##----7.3 模型诊断----
###----7.3.1 过离散-----
#观察到的方差高于模型假设的预期方差时，就会发生过度离散
library(glmmTMB)
data(Salamanders)
model <- glm(count ~ spp + mined, family = poisson, data = Salamanders)
check_overdispersion(model)

#观测到的0的数量大于预测的0时,就表示了零通货膨胀
###----7.3.2 0膨胀-----
model <- glm(count ~ spp + mined, family = poisson, data = Salamanders)
check_zeroinflation(model)

###----7.3.3 奇异模型-----
#一个"奇异"模型拟合意味着变量-协方差矩阵的某些维度估计精确为零
library(lme4)
data(sleepstudy)
# prepare data
set.seed(123)
sleepstudy$mygrp <- sample(1:5, size = 180, replace = TRUE)
sleepstudy$mysubgrp <- NA
for (i in 1:5) {
  filter_group <- sleepstudy$mygrp == i
  sleepstudy$mysubgrp[filter_group] <-
    sample(1:30, size = sum(filter_group), replace = TRUE)
}
# fit strange model
model <- lmer(
  Reaction ~ Days + (1 | mygrp / mysubgrp) + (1 | Subject),
  data = sleepstudy
)
check_singularity(model)

###----7.3.4 检查异质性-----
#线性模型假定常误差方差(同调性)
data(cars)
model <- lm(dist ~ speed, data = cars)
check_heteroscedasticity(model)

###----7.3.5 ??全面可视化check_model-----
# defining a model
model <- lm(mpg ~ wt + am + gear + vs * cyl, data = mtcars)
# checking model assumptions
check_model(model)

##----7.4 提取performance与可视化----
m1 <- lm(mpg ~ wt + cyl, data = mtcars)
model_performance(m1)
m2 <- glm(vs ~ wt + mpg, data = mtcars, family = "binomial")
model_performance(m2)
library(lme4)
m3 <- lmer(Reaction ~ Days + (1 + Days | Subject), data = sleepstudy)
model_performance(m3)
counts <- c(18, 17, 15, 20, 10, 20, 25, 13, 12)
outcome <- gl(3, 1, 9)
treatment <- gl(3, 3)
m4 <- glm(counts ~ outcome + treatment, family = poisson())

###----7.4.1 模型比较----
#compare_performance(m1, m2, m3, m4, verbose = FALSE)
compare_performance(m1, m2, m3, m4, rank = TRUE, verbose = FALSE)
###----7.4.2 ??雷达图可视化----
plot(compare_performance(m1, m2, m4, rank = TRUE, verbose = FALSE))

#----??8.performance----
library(parameters)
#parameters包支持的模型
insight::supported_models()

##----8.1 经典模型-----
model <- lm(Sepal.Width ~ Petal.Length * Species + Petal.Width, data = iris)
# 常规
model_parameters(model)
# 参数标准化
model_parameters(model, standardize = "refit")
# heteroscedasticity-consitent SE and CI
model_parameters(model, vcov = "HC3")

##----8.2 混合模型-----
library(lme4)
model <- lmer(Sepal.Width ~ Petal.Length + (1 | Species), data = iris)
# model parameters with CI, df and p-values based on Wald approximation
model_parameters(model)
# model parameters with CI, df and p-values based on Kenward-Roger approximation
model_parameters(model, ci_method = "kenward", effects = "fixed")

##----8.3 结构方程模型-----
library(psych)
model <- psych::fa(attitude, nfactors = 3)
model_parameters(model)

#select_parameters()特征选择，快速选择和保留最相关的预测
lm(disp ~ ., data = mtcars) |>
  select_parameters() |>
  model_parameters()

#----9.see----

library(ggplot2)
# step-1
model <- lm(mpg ~ factor(cyl) * wt, data = mtcars)
# step-2
results <- fortify(model)
# step-3
ggplot(results) +
  geom_point(aes(x = wt, y = mpg, color = `factor(cyl)`)) +
  geom_line(aes(x = wt, y = .fitted, color = `factor(cyl)`))

#β 系数展示
library(parameters)
library(see)
model <- lm(wt ~ am * cyl, data = mtcars)
plot(parameters(model))
 
#贝叶斯回归模型展示
library(bayestestR)
library(rstanarm)
library(see)
set.seed(123)
model <- stan_glm(wt ~ mpg, data = mtcars, refresh = 0)
result <- hdi(model, ci = c(0.5, 0.75, 0.89, 0.95))
plot(result)

#正态检验
library(performance)
library(see)
model <- lm(wt ~ mpg, data = mtcars)
check <- check_normality(model)
plot(check, type = "qq")

#omega2
library(effectsize)
library(see)
model <- aov(wt ~ am * cyl, data = mtcars)
plot(omega_squared(model))

#以模型为基础的 包计算基于模型的估计和预测
library(modelbased)
library(see)
data(mtcars)
mtcars$gear <- as.factor(mtcars$gear)
model <- lm(mpg ~ wt * gear, data = mtcars)
predicted <- estimate_expectation(model, data = "grid")
plot(predicted)
means <- estimate_means(model)
plot(means)

##---9.1 拼图----
p1 <- ggplot(iris, aes(x = Species, y = Sepal.Length, fill = Species)) +
  geom_boxplot() +
  theme_modern(axis.text.angle = 45) +
  scale_fill_material_d()

p2 <- ggplot(iris, aes(x = Species, y = Sepal.Length, fill = Species)) +
  geom_violin() +
  theme_modern(axis.text.angle = 45) +
  scale_fill_material_d(palette = "ice")

p3 <- ggplot(iris, aes(x = Petal.Length, y = Petal.Width, color = Sepal.Length)) +
  geom_point2() +
  theme_modern() +
  scale_color_material_c(palette = "rainbow")

plots(p1, p2, p3, n_columns = 2)
plots(p1, p2, p3,
      n_columns = 2,
      tags = paste("Fig. ", 1:3))

#泡泡图
normal <- ggplot(iris, aes(x = Petal.Width, y = Sepal.Length)) +
  geom_point(size = 8, alpha = 0.3) +
  theme_modern()

new <- ggplot(iris, aes(x = Petal.Width, y = Sepal.Length)) +
  geom_point2(size = 8, alpha = 0.3) +
  theme_modern()

plots(normal, new, n_columns = 2)

#半个小提琴
ggplot(iris, aes(x = Species, y = Sepal.Length, fill = Species)) +
  geom_violindot(fill_dots = "black") +
  theme_modern() +
  scale_fill_material_d()


#------10.仪表盘菜单------
library("easystats")
mymodel <- lm(Sepal.Length ~ Sepal.Width*Species, data = iris)
model_dashboard(mymodel)