初始概率小于1的生存分析：Kaplan-Meier曲线能否用于歌曲冲榜概率可视化？

Absolutely—Kaplan-Meier (KM) curves are a perfect fit here, but we need to frame your problem to align with survival analysis principles:

• Event: A song entering the Billboard Hot 100 Top 10.
• Time: Number of weeks the song remains on the full Hot 100 chart.
• Censoring: If a song drops off the chart without ever reaching the Top 10, that’s a right-censored observation (we don’t know if it would have hit the Top 10 if it stayed longer).

Your starting probability of 0.15 translates to the proportion of songs that reach the Top 10 in their first week on the chart. The KM curve will estimate the probability that a song has not yet reached the Top 10 by each week t. To get the probability of having reached the Top 10 by week t, simply calculate 1 - survival probability.

Let’s walk through a step-by-step example with simulated data (you can swap this with your actual dataset once ready):

Step 1: Install and Load the Package #

The survival package is part of R’s core stats toolkit, but load it explicitly:

install.packages("survival")  # Only run if not already installed
library(survival)

Step 2: Prepare Your Dataset #

Your data needs two key columns:

• weeks_on_chart: Number of weeks the song stayed on the Hot 100.
• reached_top10: Binary indicator (1 = song hit Top 10; 0 = song dropped off without hitting Top 10).

Let’s simulate data that matches your 0.15 starting probability:

set.seed(123)  # For reproducible results
n_songs <- 1000

# Simulate weeks on chart (1 to 20 weeks)
weeks_on_chart <- sample(1:20, n_songs, replace = TRUE)

# Simulate Top 10 hits: 15% in week 1, decreasing probability over time
reached_top10 <- ifelse(weeks_on_chart == 1, rbinom(n_songs, 1, 0.15), 
                        rbinom(n_songs, 1, 0.15 / weeks_on_chart))

# Create data frame
song_data <- data.frame(weeks_on_chart, reached_top10)

Step 3: Fit the Kaplan-Meier Model #

Use survfit() to estimate the survival function (probability of not reaching Top 10 by each week):

# Define the survival object: Surv(time, event)
surv_obj <- Surv(time = song_data$weeks_on_chart, event = song_data$reached_top10)

# Fit the KM model (~1 estimates a single curve for all songs)
km_fit <- survfit(surv_obj ~ 1)

Step 4: Plot the Probability of Reaching Top 10 #

Transform the survival probability to get the event probability (reaching Top 10) and plot:

plot(km_fit, 
     main = "Probability of Reaching Billboard Hot 100 Top 10 by Weeks on Chart",
     xlab = "Weeks on Billboard Hot 100",
     ylab = "Probability of Reaching Top 10",
     fun = function(x) 1 - x,  # Convert survival to event probability
     col = "darkblue",
     lwd = 2)

# Add a reference line for your starting probability
abline(h = 0.15, col = "red", lty = 2, lwd = 1.5)
text(x = 5, y = 0.16, labels = "Week 1 Probability (0.15)", col = "red")

Customization Tips #

• To add subgroups (e.g., genre, artist type), modify the model to surv_obj ~ genre to plot multiple curves.
• Add confidence intervals with conf.int = TRUE in survfit().
• For ggplot2-style plots, use the ggsurvfit package (a wrapper for survival that integrates with ggplot).

• Censoring Accuracy: Ensure you correctly mark songs that left the chart without hitting Top 10 as censored—this is critical for reliable KM estimates.
• Starting Probability: If your 0.15 value is hypothetical (not from real data), you may need to adjust your dataset or use parametric survival models (like Weibull) to enforce this starting point.

内容的提问来源于stack exchange，提问作者Sander