R语言绘制经验分布函数以及生存函数

关键词：R语言，经验分布函数，生存函数

经验分布函数的数学定义：

R自定义函数--经验分布函数：

scores <- c(63, 72, 74, 79, 82, 82, 87, 89, 90, 90)


# R的自定义函数，函数名为cdf_table
cdf_table <- function(x){


x <- sort(x) #将x按照从小到大的顺序排列，不删掉重复的值。
n <- length(x) #x的长度
tab <- unname(c(table(x))) #x每个数出现的次数，即频数。
pct = tab / n #计算频率


d <- data.frame( #形成一个数据框
x = sort(unique(x)), #删掉x中重复的值，确保x中每个数都是唯一的，然后对这些唯一的数，按照从小到大的顺序进行排列
freq = tab, #变量名为freq，该列中给出x中数值出现的频数(出现的次数)
pct = pct, #变量名为pct，该列给出x中数值对应的频率
cumfreq = cumsum(tab), #变量名为cumfreq，该列中给出频数的累加和，即，该行的第二个值=tab的第一个值，加上 tab的第二个值。
cumpct = cumsum(pct)) #变量名为cumpct，该列中给出频率的累加和。


d #该函数返回d（数据框）
}


knitr::kable(cdf_table(scores)) #生成表格

得到如下表格：

其中

• x表示原始数据（1）先unique，即删掉重复的数值（2）再将剩下的数值按照从小到大的顺序进行排列。
• freq表示x对应的频数；
• pct表示x对应的频率；
• cumfreq表示x对应的频数累加和；
• cumpct表示x对应的频率累加和。

绘制经验分布函数：

d <- data.frame(
scores = c(63, 72, 74, 79, 82, 82, 87, 89, 90, 90)) #原始数据，数据类型为数据框；


dfreq <- cdf_table(d$scores) #返回数据对应的数据框，包含原始数据的取值，每个值对应的频数，频率，频数累加和，频率累加和

p <- ggplot(data = dfreq, mapping = aes(
x = x, y = cumpct )) #基于作图的数据为dfreq（数据框）, mapping给出了数据要映射的x轴和y轴数据

#对上图做进一步细化设定
p + stat_ecdf(geom = "step") +
scale_y_continuous( 
limits = c(-.01, 1.01),  #调整y轴坐标轴的显示范围
expand = c(0, 0), #expand = c(0,0)不放缩，这个命令自己可以尝试修改不同的数字，看下图的效果
name = " 累计百分比") #y轴显示的标签名称

 绘制上图对应的生存函数：

surv_table <- function(x){
  x <- sort(x)
  n <- length(x)
  tab <- unname(c(table(x)))
  pct = tab / n
  d <- data.frame(
    x = sort(unique(x)),
    freq = tab,
    pct = pct,
    cumfreq = cumsum(tab),
    cumpct = cumsum(pct),
    cumsurv = 1-cumsum(pct))
  d
}
knitr::kable(surv_table(scores))


  scores = c(63, 72, 74, 79, 82, 82, 87, 89, 90, 90)


surv_req <-surv_table(scores)
f1 <- c(-Inf,0,0,0,0,1) # 手动增加-Inf处的取值 
f2 <- c(Inf,0,0,0,1,0)  #手动增加Inf处的取值
surv_req <- rbind(f1,surv_req,f2)

p <- ggplot(data = surv_req, mapping = aes(
  x = x, y = cumsurv ))
p + geom_step() +
  scale_y_continuous(
    limits = c(-.01, 1.01),
    expand = c(0, 0),
    name = " 累计百分比")

在绘制生存函数的时候，这里使用了geom_step()函数，需要补充上，阶梯函数在-Inf和Inf处的取值

其他一些可以参考的知识：

关于ggplot函数，可以参考北京大学李东风老师编写的Rbook。下图是李东风老师关于ggplot作图的介绍：

R-base包中的table()函数：

table(x)返回x的每个不同值的频率（出现次数），table函数用于创建列联表，用于统计x中每个不同取值的频率。返回的是一个table类。

unname()函数，用于删除R对象中，变量的名称。

## Stat_ecdf 原函数
stat_ecdf <- function(mapping = NULL, 
                      data = NULL,
                      geom = "step",
                      position = "identity",
                      ...,
                      n = NULL,
                      pad = TRUE,
                      na.rm = FALSE,
                      show.legend = NA,
                      inherit.aes = TRUE) {
  layer(
    data = data,
    mapping = mapping,
    stat = StatEcdf,
    geom = geom,
    position = position,
    show.legend = show.legend,
    inherit.aes = inherit.aes,
    params = list2(
      n = n,
      pad = pad,
      na.rm = na.rm,
      ...
                   )
          )
                                            }


#' @rdname ggplot2-ggproto
#' @format NULL
#' @usage NULL
#' @export
        StatEcdf <- ggproto("StatEcdf", Stat,
          required_aes = c("x|y"),

          default_aes = aes(x = after_stat(ecdf), y = after_stat(ecdf)),

  setup_params = function(self, data, params) {
        params$flipped_aes <- has_flipped_aes(data, params, main_is_orthogonal = FALSE, main_is_continuous = TRUE)

    has_x <- !(is.null(data$x) && is.null(params$x))
    has_y <- !(is.null(data$y) && is.null(params$y))
    if (!has_x && !has_y) {
      cli::cli_abort("{.fn {snake_class(self)}} requires an {.field x} or {.field y} aesthetic.")
    }

    params
  },

  compute_group = function(data, scales, n = NULL, pad = TRUE, flipped_aes = FALSE) {
    data <- flip_data(data, flipped_aes)
    # If n is NULL, use raw values; otherwise interpolate
    if (is.null(n)) {
      x <- unique0(data$x)
    } else {
      x <- seq(min(data$x), max(data$x), length.out = n)
    }

    if (pad) {
      x <- c(-Inf, x, Inf)
    }
    data_ecdf <- ecdf(data$x)(x)

    df_ecdf <- data_frame0(
      x = x,
      y = data_ecdf,
      ecdf = data_ecdf,
      .size = length(x)
    )
    df_ecdf$flipped_aes <- flipped_aes
    flip_data(df_ecdf, flipped_aes)
  }
)

推荐常用R语言包及其函数检索网址：R Package Documentation

李东风老师R书，R语言教程 (pku.edu.cn)