Spotify 流量歌曲的音乐特征与受众趋势分析

R
EDA
音乐
可视化
抽样
对 2023 年 Spotify 播放量前列的歌曲进行探索性分析，揭示音乐流行趋势。
作者

Pinn Xu

发布于

2026年2月25日

10 min read · 302 EN words · 2374 中文字
环境与版本 
R==4.5.2
plotly==4.12.0
sampling==2.11
tidyverse==2.0.0

1 数据集概览

本项目使用来自 Kaggle 的数据集 Most Streamed Spotify Songs 2023，包含 2023 年 Spotify 最受欢迎的 953 首歌曲，属性涵盖曲目名称、艺术家姓名、发行日期、Spotify 播放列表与榜单、播放统计、Apple Music 收录情况、Deezer 收录情况、Shazam 榜单及若干音频特征。

本项目使用的核心数值变量包括：streams（播放量）、BPM（每分钟拍数）、valence index（情绪效价指数）以及在 Spotify 和 Deezer 播放列表中的收录次数。streams 表示一首歌曲在 2023 年全年的播放次数，应为正整数；收录次数同样为计数变量（非负整数）。BPM 决定音乐的节奏与速度。情绪效价指数表示歌曲音乐内容的积极程度，编码范围为 0 到 100（0 = 完全消极，100 = 完全积极）。

类别变量包括音乐的调式（mode）与调性（key）。调式分为大调（major）和小调（minor）两类，是两种通用的和弦类别；调性涵盖 C、C#、D、D# 等，决定了歌曲所使用的基础音符与和弦组合。

总体而言，由于变量的维度与类型各异，本数据集被认为适合后续分析。

2 分析目标

本项目旨在分析 2023 年 Spotify 最受欢迎歌曲的音乐特征的显著规律与趋势，并探究各指标之间的潜在关系。此外，将通过不同抽样方法对数据进行抽样，以验证中心极限定理（CLT）的适用性，并比较各抽样方法对总体参数的预测能力。

3 数据准备

在正式分析前，数据集首先从 Kaggle 网站下载并存储至本地文件夹，随后使用 read.csv() 函数导入 R。以下展示了导入后数据集作为 DataFrame 对象的结构。

代码 
# 导入数据集
spotify_data <- read.csv("../../../src/data/spotify-2023.csv")
str(spotify_data)
'data.frame':   953 obs. of  24 variables:
 $ track_name          : chr  "Seven (feat. Latto) (Explicit Ver.)" "LALA" "vampire" "Cruel Summer" ...
 $ artist.s._name      : chr  "Latto, Jung Kook" "Myke Towers" "Olivia Rodrigo" "Taylor Swift" ...
 $ artist_count        : int  2 1 1 1 1 2 2 1 1 2 ...
 $ released_year       : int  2023 2023 2023 2019 2023 2023 2023 2023 2023 2023 ...
 $ released_month      : int  7 3 6 8 5 6 3 7 5 3 ...
 $ released_day        : int  14 23 30 23 18 1 16 7 15 17 ...
 $ in_spotify_playlists: int  553 1474 1397 7858 3133 2186 3090 714 1096 2953 ...
 $ in_spotify_charts   : int  147 48 113 100 50 91 50 43 83 44 ...
 $ streams             : chr  "141381703" "133716286" "140003974" "800840817" ...
 $ in_apple_playlists  : int  43 48 94 116 84 67 34 25 60 49 ...
 $ in_apple_charts     : int  263 126 207 207 133 213 222 89 210 110 ...
 $ in_deezer_playlists : chr  "45" "58" "91" "125" ...
 $ in_deezer_charts    : int  10 14 14 12 15 17 13 13 11 13 ...
 $ in_shazam_charts    : chr  "826" "382" "949" "548" ...
 $ bpm                 : int  125 92 138 170 144 141 148 100 130 170 ...
 $ key                 : chr  "B" "C#" "F" "A" ...
 $ mode                : chr  "Major" "Major" "Major" "Major" ...
 $ danceability_.      : int  80 71 51 55 65 92 67 67 85 81 ...
 $ valence_.           : int  89 61 32 58 23 66 83 26 22 56 ...
 $ energy_.            : int  83 74 53 72 80 58 76 71 62 48 ...
 $ acousticness_.      : int  31 7 17 11 14 19 48 37 12 21 ...
 $ instrumentalness_.  : int  0 0 0 0 63 0 0 0 0 0 ...
 $ liveness_.          : int  8 10 31 11 11 8 8 11 28 8 ...
 $ speechiness_.       : int  4 4 6 15 6 24 3 4 9 33 ...

随后，导入了本项目所使用的主要包（包括 tidyverse、plotly 和 sampling），并使用 tibble() 函数将原始数据集转换为 tibble（df），以便后续数据处理。

以下展示了数据集的初步概览，包括样本量（N = 953）及是否存在缺失值（使用 anyNA() 函数获取缺失值信息）。

代码 
# 加载包
library(readr)
library(tidyverse)
library(sampling)
library(plotly)

# 将 dataframe 转换为 tibble
df <- tibble(spotify_data)

# 样本量
cat("Sample Size:", nrow(df), '\n')

# 检查缺失值
cat("Whether there is NA:", anyNA(df), "\n")
Sample Size: 953 
Whether there is NA: FALSE

3.1 Streams（播放量）

如下所示，streams 变量的数据类型为字符型，而播放量应为数值型而非字符串，因此需要进行类型转换。

代码 
# 查看结构
str(df["streams"])
tibble [953 × 1] (S3: tbl_df/tbl/data.frame)
 $ streams: chr [1:953] "141381703" "133716286" "140003974" "800840817" ...

在将数据类型转换为数值型之前，首先检查是否存在无法识别为数值的非法字符串，使用正则表达式进行数据筛选。结果显示，有 1 条观测的播放量值无效。

代码 
# 查看无效值
cat("Invalid Streams: \n\n")
df |>
  select(streams) |>
  filter(!str_detect(streams, "^(0|[1-9][0-9]*)$"))
Invalid Streams: 

# A tibble: 1 × 1
  streams                                                                       
  <chr>                                                                         
1 BPM110KeyAModeMajorDanceability53Valence75Energy69Acousticness7Instrumentalne…

随后，仅保留播放量有效的观测，并将数据类型转换为数值型。以下展示了转换后缺失值的情况及变量结构。

代码 
# 仅保留有效值
df <- df |>
  filter(str_detect(streams, "^[1-9][0-9]*"))

# 将数据类型强制转换为数值型
df <- df |>
  mutate(streams = as.numeric(streams))

# 检查类型转换后 streams 中是否存在 NA
cat("Whether there is NA in streams:", anyNA(df$streams), "\n\n")

# 新结构
str(df["streams"])
Whether there is NA in streams: FALSE 

tibble [952 × 1] (S3: tbl_df/tbl/data.frame)
 $ streams: num [1:952] 1.41e+08 1.34e+08 1.40e+08 8.01e+08 3.03e+08 ...

3.2 Deezer 播放列表收录次数

同样，Deezer 播放列表收录次数变量的原始类型也为字符型。再次使用 str_detect() 检查是否存在格式无效的数值。结果显示，存在若干逗号分隔格式的数值，且无其他格式的无效数值，因此在类型转换前需进行进一步处理。

代码 
# 查看不是 0 或纯数字的值
cat("Non-Number Values:\n\n")
df |>
  select(in_deezer_playlists) |>
  filter(!str_detect(in_deezer_playlists, "^(0|[1-9][0-9]*)$"))
Non-Number Values:

# A tibble: 79 × 1
   in_deezer_playlists
   <chr>              
 1 2,445              
 2 3,394              
 3 3,421              
 4 4,053              
 5 1,056              
 6 4,095              
 7 1,003              
 8 1,800              
 9 2,703              
10 1,632              
# ℹ 69 more rows
代码 
# 应用新模式：逗号分隔的数字
cat("Neither Pure Number Nor Comma-Separated:\n\n")
df |>
  select(in_deezer_playlists) |>
  filter(!str_detect(
    in_deezer_playlists, "^(0|[1-9][0-9]*|[1-9][0-9]*,[0-9]*)$")
  )
Neither Pure Number Nor Comma-Separated:

# A tibble: 0 × 1
# ℹ 1 variable: in_deezer_playlists <chr>

将所有逗号替换为空字符串后，对变量进行数值类型转换，并检查转换后是否存在缺失值。结果显示无缺失值。

代码 
# 将逗号替换为空字符串
df <- df |>
  mutate(in_deezer_playlists = str_replace_all(in_deezer_playlists, ",", ""))

# 将数据类型强制转换为数值型
df <- df |>
  mutate(in_deezer_playlists = as.numeric(in_deezer_playlists))

# 检查类型转换后是否存在 NA
cat("Whether there is NA in in_deezer_playlists:",
    anyNA(df$in_deezer_playlists), "\n"
    )
Whether there is NA in in_deezer_playlists: FALSE

3.3 调性（Key）

对于调性变量，首先检查其所有类别。结果显示存在无调性信息的空字符串，因此删除这些无调性信息的歌曲。删除后再次检查类别以确保数据有效。

代码 
# 检查 key 变量的类别
cat("Categories of Key:\n")
unique(df$key)

# 删除 key 为空的行
df <- df[df$key != "", ]

# 删除后再次检查
cat("\nCategories After Removal:\n")
unique(df$key)
Categories of Key:
 [1] "B"  "C#" "F"  "A"  "D"  "F#" ""   "G#" "G"  "E"  "A#" "D#"

Categories After Removal:
 [1] "B"  "C#" "F"  "A"  "D"  "F#" "G#" "G"  "E"  "A#" "D#"

4 单变量分析

4.1 调式分布

如下饼图所示，本数据集中两种音乐调式的占比相对均衡。具体而言，大调占比 55.3%，小调紧随其后，占比 44.7%。10.6% 的差距表明，通常与明亮情感特质相关的大调歌曲更为普遍，但小调歌曲仍占据相当比例。

代码 
# 饼图
plot_ly(
  as.data.frame(table(df["mode"])),
  labels = ~mode, values = ~Freq, type = "pie",
  textinfo = "label+percent",
  marker = list(
    colors = c("khaki", "paleturquoise"),
    line = list(width = 2))
) |>
  layout(
    title = list(y = 0.9, text = "音乐调式分布"), showlegend = FALSE,
    margin = list(t = 80, l = 60, r = 60)
  )

4.2 BPM 分布

在 BPM（每分钟拍数）方面，如下直方图所示，其分布在 110–130 BPM 附近达到峰值。大多数歌曲的 BPM 介于 90 到 150 之间，表明大众偏好适中节奏的音乐。

代码 
# BPM 分布
plot_ly(
  df, x = ~bpm, type = "histogram",
  marker = list(
    color = "pink", line = list(width = 1.5, color = "white")
  )
) |>
  layout(
    margin = list(t = 80),
    title = list(y = 0.9, text = "BPM 分布"),
    xaxis = list(
      title = "每分钟拍数（BPM）",
      linewidth = 2,
      dtick = 20,
      ticks = "outside",
      ticklen = 5,
      tickwidth = 1
    ),
    yaxis = list(title = "频数")
  )

5 多变量分析

5.1 播放量最高的前五调性

如下条形图所示，在播放量方面，C# 为最受欢迎的调性，平均播放量最高（约 6 亿次）。E 位居第二，D# 与 A# 均值相近紧随其后。D 在前五中最低。总体而言，各调性之间的差异较为适中，表明这些调性在流行度上较为均衡。

代码 
# 找出平均播放量最高的前五个调性
top5_key <- df |>
  group_by(key) |>
  summarize(mean_streams = mean(streams), .groups = "drop") |>
  arrange(desc(mean_streams)) |>
  head(5) |>
  mutate(key = reorder(key, -mean_streams))  # 将 key 转换为因子
代码 
# 条形图
plot_ly(
  top5_key, x = ~key, y = ~mean_streams, type = "bar",
  marker = list(color = "mediumaquamarine")
) |>
  layout(
    title = list(y = 0.9, text = "最受欢迎的前五个调性"),
    xaxis = list(title = "调性", linewidth = 2, range = c(-0.8, 4.8)),
    yaxis = list(title = list(text = "平均播放量", standoff = 10)),
    margin = list(t = 80, l = 80, r = 20)
  )

5.2 不同调式的情绪效价分布

箱线图展示了不同调式类别下情绪效价（valence）的分布。小调歌曲的情绪效价中位数（54）略高于大调歌曲（48），两者的分布范围较宽，变异程度相似，四分位距存在重叠。尽管小调的平均情绪效价略高，但大量重叠表明单凭调式并不能强烈预测情绪效价。

代码 
# 箱线图
plot_ly(
  df, x = ~mode, y = ~valence_., type = "box",
  fillcolor = "thistle",
  line = list(width = 2, color = "darkslateblue")
) |>

  layout(
    margin = list(t = 80),  # 顶部边距
    title = list(
      y = 0.9,
      text = "不同调式的情绪效价分布"
    ),
    xaxis = list(title = "音乐调式"),
    yaxis = list(title = "情绪效价")
  )

5.3 Spotify 与 Deezer 播放列表收录次数

如下散点图展示了歌曲在 Spotify 与 Deezer 播放列表中收录次数的关系，呈现出正相关趋势。随着歌曲在一个平台上播放量的增加，在另一个平台上的出现频率也随之提升。总体而言，数据最集中于 Spotify 收录次数 1 万以下、Deezer 收录次数 2 千以下的区间，少数高表现的异常值在 Spotify 超过 5 万次、Deezer 超过 1.2 万次。

代码 
# 散点图
plot_ly(df, x = ~in_spotify_playlists) |>
  add_markers(
    y = ~in_deezer_playlists,
    marker = list(color = "salmon", opacity = 0.6)
  ) |>
  layout(
    margin = list(t = 80),
    title = list(
      y = 0.9,
      text = "Spotify 与 Deezer 播放列表收录次数的关系"
    ),
    xaxis = list(title = "Spotify 播放列表收录次数"),
    yaxis = list(title = "Deezer 播放列表收录次数")
  )

6 中心极限定理

中心极限定理指出，随着样本量的增加，样本均值的分布将越来越接近正态分布（对称钟形），无论总体分布的形状如何。

本节从 df 中的 bpm 变量中抽取不同样本量（n = 10、20、30、40）的样本，并重复抽样 1000 次以获取样本均值分布。下方的直方图矩阵在相同的 x 轴和 y 轴范围下展示了四种样本量的样本均值分布，以便比较其形状。矩阵之后是所有四种样本均值分布与总体统计量的对比输出。

结果表明，随着样本量的增加，样本均值的分布显著更加集中，标准差更小。同时，所有四种样本均值的均值均接近总体均值，体现了样本均值分布对总体的预测能力。此外，所有结果中的样本均值分布均呈对称形态，与正态分布形状相似。

综上，本节结果与中心极限定理的预测一致。

代码 
# 定义变量
sizes <- c(10, 20, 30, 40) # 样本量
sample_means <- numeric(1000) # 样本均值
means_sample_means <- numeric(length(sizes)) # 样本均值的均值
sds_sample_means <- numeric(length(sizes)) # 样本均值的标准差
p <- list() # 图形

# 重复抽样 1000 次
set.seed(123)
for (i in 1:length(sizes)) {

  # 抽取样本并计算均值
  sample_means <- replicate(
    1000, mean(sample(df$bpm, size = sizes[i], replace = TRUE), na.rm = TRUE)
  )

  # 存储样本均值的均值与标准差
  means_sample_means[i] <- mean(sample_means, na.rm = TRUE)
  sds_sample_means[i] <- sd(sample_means, na.rm = TRUE)

  # 存储图形信息
  p[[i]] <- plot_ly(
    x = sample_means, type = "histogram",
    xbins = list(size = 2),
    marker = list(color = "pink")
  ) |>
    layout(
      xaxis = list(
        title = paste("n =", sizes[i]),
        range = c(90, 150),
        dtick = 15,
        ticks = "outside",
        ticklen = 5,
        tickwidth = 1
      ),
      yaxis = list(range = c(0, 180), dtick = 60)
    )
}

# 显示图形
subplot(p, nrows = 2, margin = 0.1, titleX = TRUE) |>
  layout(
    title = list(y = 0.9, text = "BPM 样本均值分布"),
    showlegend = FALSE,
    margin = list(t = 80)
  )
代码 
# 总体均值与标准差
pop_mean <- mean(df$bpm, na.rm = TRUE)
pop_sd <- sd(df$bpm, na.rm = TRUE)

# 显示结果
cat(paste("Population Mean:", round(pop_mean, 2)), "\n")
cat(paste("Population SD:", round(pop_sd, 2)), "\n\n")
cat("Means of Sample Means:\n")
cat(paste0(
    round(means_sample_means, 2), " (n=", seq(10, 40, 10), ")",
    collapse = "\n"
  ))
cat("\n\n")
cat("SDs of Sample Means:\n")
cat(paste0(
    round(sds_sample_means, 2), " (n=", seq(10, 40, 10), ")",
    collapse = "\n"
  ))
Population Mean: 122.84 
Population SD: 28.2 

Means of Sample Means:
122.99 (n=10)
123.27 (n=20)
123.09 (n=30)
122.65 (n=40)

SDs of Sample Means:
8.85 (n=10)
6.41 (n=20)
5.25 (n=30)
4.34 (n=40)

7 抽样方法

7.1 抽样

本节使用并比较三种抽样方法：有放回的简单随机抽样、以 bpm 为概率的不等概率系统抽样，以及按 mode 分层的无放回分层抽样。样本量设定为 100，并使用 sampling 包中的函数进行抽样。

代码 
# 有放回简单随机抽样
n <- 100
set.seed(123)
s <- srswr(n, nrow(df))
rows <- (1:nrow(df))[s!=0]
rows <- rep(rows, s[s!=0])
sample1 <- df[rows, ]


# 不等概率系统抽样
set.seed(123)
pik <- inclusionprobabilities(df$bpm, n) # 以 bpm 计算抽中概率
s <- UPsystematic(pik)
sample2 <- df[s!= 0, ]


# 无放回分层抽样
# 按 mode 设置层
order_index <- order(df$mode)
data <- df[order_index, ]
prop_table <- table(data$mode) / nrow(data)
sizes <- round(prop_table * n)
# 抽样
set.seed(32)
st <- sampling::strata(
  data, stratanames = "mode", size = sizes, method = "srswor"
)
sample3 <- sampling::getdata(data, st)

7.2 方法比较

以下展示了从每个样本和总体中计算出的调式类别比例。如结果所示，分层抽样方法在预测总体调式比例方面表现最佳，因为它在抽样前根据调式类别划定层次。相比之下，简单随机抽样的表现最差，这可能是由于其随机性较强、缺乏分层抽样那样的结构性规律。此外，系统抽样的预测效果居中，优于简单随机抽样但不及分层抽样。这可能源于系统抽样按固定间隔抽取样本，样本在数据集中分布相对均匀，有助于提高类别属性比例的准确性；与此同时，以 bpm 为基础的不等概率选取也可能影响 mode 比例的分布规律。

代码 
# 调式比例
cat("Proportion of Modes in Population:")
(pct0 <- round(table(df$mode) / nrow(df), 2)); cat("\n")

cat("Proportion of Modes in Sample 1 (Simple Random):")
(pct1 <- table(sample1$mode) / nrow(sample1)); cat("\n")
cat("Proportion of Modes in Sample 2 (Systematic):")
(pct2 <- table(sample2$mode) / nrow(sample2)); cat("\n")
cat("Proportion of Modes in Sample 3 (Stratified):")
(pct3 <- table(sample3$mode) / nrow(sample3))
Proportion of Modes in Population:
Major Minor 
 0.55  0.45 

Proportion of Modes in Sample 1 (Simple Random):
Major Minor 
 0.63  0.37 

Proportion of Modes in Sample 2 (Systematic):
Major Minor 
 0.46  0.54 

Proportion of Modes in Sample 3 (Stratified):
Major Minor 
 0.55  0.45
回到顶部