R语言作为统计分析和数据科学领域的首选编程语言，其强大的统计函数库是其核心竞争力。本文将从基础统计函数出发，深入探讨高级应用技巧，并针对数据分析中的常见问题提供解决方案。
一、基础统计函数详解
1.1 描述性统计函数
R语言提供了丰富的基础统计函数，这些函数是数据分析的基石。让我们从最基础的描述性统计开始。
基本统计量计算
# 创建示例数据
set.seed(123) # 设置随机种子以确保结果可重现
data <- rnorm(100, mean = 50, sd = 10) # 生成100个正态分布的随机数
# 基础统计函数
mean_val <- mean(data) # 算术平均值
median_val <- median(data) # 中位数
sd_val <- sd(data) # 标准差
var_val <- var(data) # 方差
min_val <- min(data) # 最小值
max_val <- max(data) # 最大值
range_val <- range(data) # 极差（最小值和最大值）
# 输出结果
cat("平均值:", mean_val, "\n")
cat("中位数:", median_val, "\n")
cat("标准差:", sd_val, "\n")
cat("方差:", var_val, "\n")
cat("最小值:", min_val, "\n")
cat("最大值:", max_val, "\n")
cat("极差:", range_val[1], "-", range_val[2], "\n")
代码解释
：
set.seed(123)
：设置随机数生成器的种子，确保每次运行代码时生成相同的随机数，便于结果重现。
rnorm()
：生成正态分布随机数，参数包括样本量、均值和标准差。
mean()
、
median()
、
sd()
、
var()
：分别计算算术平均值、中位数、标准差和方差。
min()
、
max()
、`R语言作为统计分析和数据科学领域的首选编程语言，其强大的统计函数库是其核心竞争力。本文将从基础统计函数出发，深入探讨高级应用技巧，并针对数据分析中的常见问题提供解决方案。
一、基础统计函数详解
1.1 描述性统计函数
R语言提供了丰富的基础统计函数，这些函数是数据分析的基石。让我们从最基础的描述性统计开始。
基本统计量计算
# 创建示例数据
set.seed(123) # 设置随机种子以确保结果可重现
data <- rnorm(100, mean = 50, sd = 10) # 生成100个正态分布的随机数
# 基础统计函数
mean_val <- mean(data) # 算术平均值
median_val <- median(data) # 中位数
sd_val <- sd(data) # 标准差
var_val <- var(data) # 方差
min_val <- min(data) # 最小值
max_val <- max(data) # 最大值
range_val <- range(data) # 极差（最小值和最大值）
# 输出结果
cat("平均值:", mean_val, "\n")
cat("中位数:", median_val, "\1. **数据准备**：首先创建了一个包含100个元素的向量`data`，这些数据服从均值为50、标准差为10的正态分布。
2. **统计量计算**：使用R内置函数分别计算了算术平均值、中位数、标准差、方差、最小值、最大值和极差。
3. **结果输出**：使用`cat()`函数将结果格式化输出，便于阅读。
#### 描述性统计汇总
对于更全面的描述性统计，我们可以使用`summary()`函数或`psych`包中的`describe()`函数：
```r
# 使用summary()函数
summary_stats <- summary(data)
print(summary_stats)
# 使用psych包进行更详细的描述性统计
# install.packages("psych") # 如果未安装，先安装包
library(psych)
detailed_stats <- describe(data)
print(detailed_stats)
输出结果示例
：
vars n mean sd median trimmed mad min max range skew kurtosis se
1 1 100 50.2 9.84 50.15 50.23 9.78 26.81 73.89 47.08 -0.06 -0.23 0.98
关键指标解释
：
vars
：变量编号
n
：样本量
mean
：均值
sd
：标准差
median
：中位数
trimmed
：修剪均值（去除极端值后的均值）
mad
：中位数绝对偏差
min/max
：最小/最大值
range
：极差
skew
：偏度（衡量分布不对称性）
kurtosis
：峰度（衡量分布尖锐程度）
se
：标准误
1.2 频率统计与交叉表
在处理分类数据时，频率统计和交叉表是必不可少的工具。
频率统计
# 创建分类数据
category_data <- sample(c("A", "B", "C", "D"), 100, replace = TRUE,
prob = c(0.2, 0.3, 0.25, 0.25))
# 基础频率统计
freq_table <- table(category_data)
print(freq_table)
# 比例统计
prop_table <- prop.table(freq_table)
print(prop_table)
# 带百分比的频率表
freq_with_pct <- addmargins(prop.table(freq_table), 1)
print(freq_with_pct)
代码解释
：
sample()
：从指定类别中随机抽取样本，
prob
参数指定每个类别的概率。
table()
：计算频数。
prop.table()
：计算比例（默认按列，参数1表示按行）。
addmargins()
：添加汇总行（参数1表示添加行汇总）。
交叉表
# 创建两个分类变量
gender <- sample(c("男", "女"), 100, replace = TRUE)
grade <- sample(c("优秀", "良好", "及格"), 100, replace = TRUE)
# 交叉表
cross_table <- table(gender, grade)
print(cross_table)
# 带百分比的交叉表
# 行百分比
row_pct <- prop.table(cross_table, 1)
print(round(row_pct, 3))
# 列百分比
col_pct <- prop.table(cross_table, 2)
print(round(col_pct, 3))
# 总体百分比
total_pct <- prop.table(cross_table)
print(round(total_pct, 3))
应用场景
：交叉表常用于市场分析、用户行为分析等场景，帮助理解两个分类变量之间的关系。
1.3 基础假设检验
t检验
# 单样本t检验
t.test(data, mu = 50) # 检验均值是否等于50
# 独立样本t检验
group1 <- rnorm(50, mean = 48, sd = 8)
group2 <- rnorm(50,方差分析（ANOVA）是用于比较三个或更多组均值差异的统计方法。R语言提供了多种方差分析函数，最常用的是`aov()`函数。
```r
# 创建多组数据
group_a <- rnorm(30, mean = 50, sd = 8)
group_b <- rnorm(30, mean = 55, sd = 8)
group_c <- rnorm(30, mean = 60, sd = 8)
# 合并数据
data_anova <- data.frame(
value = c(group_a, group_b, group_c),
group = factor(rep(c("A", "B", "C"), each = 30))
)
# 单因素方差分析
anova_result <- aov(value ~ group, data = data_anova)
summary(anova_result)
# 多因素方差分析（示例）
# 假设有两个因素：治疗方法和性别
treatment <- factor(rep(c("DrugA", "DrugB", "Placebo"), each = 20))
gender <- factor(rep(c("M", "F"), each = 10, times = 3))
response <- rnorm(60, mean = 50, sd = 10)
multi_anova <- aov(response ~ treatment * gender)
summary(multi_anova)
代码解释
：
aov()
：拟合方差分析模型，公式格式为
因变量 ~ 自变量
。
summary()
：显示方差分析表，包括F值、p值等。
treatment * gender
：表示主效应和交互效应。
非参数检验
当数据不满足正态分布假设时，可以使用非参数检验：
# Wilcoxon秩和检验（曼-惠特尼U检验）
wilcox.test(group1, group2)
# Kruskal-Wallis检验（多组比较）
kruskal.test(value ~ group, data = data_anova)
# Friedman检验（重复测量）
# 创建重复测量数据
subject <- rep(1:20, each = 3)
time <- rep(c("T1", "T2", "T3"), 20)
score <- rnorm(60, mean = 50, sd = 5)
friedman.test(score ~ time | subject)
二、高级统计函数应用
2.1 回归分析
线性回归
# 创建模拟数据
set.seed(123)
n <- 100
x1 <- rnorm(n, 10, 2)
x2 <- rnorm(n, 5, 1)
error <- rnorm(n, 0, 3)
y <- 20 + 3*x1 + 2*x2 + error
# 简单线性回归（单个自变量）
simple_lm <- lm(y ~ x1)
summary(simple_lm)
# 多元线性回归
multi_lm <- lm(y ~ x1 + x2)
summary(multi_lm)
# 带交互项的回归
interaction_lm <- lm(y ~ x1 * x2)
summary(interaction_lm)
# 回归诊断
par(mfrow = c(2, 2))
plot(multi_lm)
回归结果解读
：
Coefficients
：回归系数估计值
Std. Error
：标准误
t value
：t统计量
Pr(>|t|)
：p值（<0.05通常认为显著）
R-squared
：决定系数，解释模型拟合优度
广义线性模型（GLM）
# 逻辑回归（二分类）
# 创建二分类响应变量
y_binary <- ifelse(y > median(y), 1, 0)
logit_model <- glm(y_binary ~ x1 + x2, family = binomial(link = "logit"))
summary(logit_model)
# 预测概率
predicted_probs <- predict(logit_model, type = "response")
# 泊松回归（计数数据）
count_data <- rpois(n, lambda = 5 + 0.5*x1)
poisson_model <- glm(count_data ~ x1, family = poisson)
summary(poisson_model)
2.2 时间序列分析
# 创建时间序列数据
ts_data <- ts(rnorm(100, mean = 10, sd = 2), start = c(2020, 1), frequency = 12)
# 基本时间序列分析
plot(ts_data, main = "时间序列图")
abline(h = mean(ts_data), col = "red", lty = 2)
# 自相关函数（ACF）
acf(ts_data, main = "自相关函数")
# 偏自相关函数（PACF）
pacf(ts_data, main = "偏自相关函数")
# 时间序列分解
decomposed <- decompose(ts_data)
plot(decomposed)
# ARIMA模型
# 自动选择最佳ARIMA模型
library(forecast)
auto_arima_model <- auto.arima(ts_data)
summary(auto_arima_model)
# 预测未来12个时间点
forecast_result <- forecast(auto_arima_model, h = 12)
plot(forecast_result)
2.3 主成分分析（PCA）
# 创建高维数据
set.seed(123)
high_dim_data <- matrix(rnorm(100*10), nrow = 100, ncol = 10)
colnames(high_dim_data) <- paste0("Var", 1:10)
# 执行PCA
pca_result <- prcomp(high_dim_data, scale. = TRUE)
# 查看结果
summary(pca_result)
# 特征值（方差解释）
pca_result$sdev^2
# 载荷矩阵（变量与主成分的关系）
print(pca_result$rotation)
# 得分矩阵（样本与主成分的关系）
scores <- pca_result$x
# 可视化
biplot(pca_result, scale = 0)
三、常见问题与解决方案
3.1 缺失值处理
缺失值是数据分析中常见的问题，R语言提供了多种处理方法。
识别缺失值
# 创建含缺失值的数据
data_with_na <- c(1, 2, NA, 4, NA, 6, 7)
# 识别缺失值
is.na(data_with_na) # 返回逻辑向量
which(is.na(data_with_na)) # 返回缺失值位置
# 统计缺失值数量
sum(is.na(data_with_na))
处理缺失值
# 1. 删除缺失值
complete_data <- na.omit(data_with_na) # 删除含NA的行/列
# 或
clean_data <- data_with_na[!is.na(data_with_na)]
# 2. 用均值填充
mean_imputed <- data_with_na
mean_imputed[is.na(mean_imputed)] <- mean(mean_imputed, na.rm = TRUE)
# 3. 用中位数填充
median_imputed <- data_with_na
median_imputed[is.na(median_imputed)] <- median(median_imputed, na.rm = TRUE)
# 4. 用众数填充（适用于分类数据）
get_mode <- function(v) {
uniqv <- unique(v)
uniqv[which.max(tabulate(match(v, uniqv)))]
}
mode_imputed <- data_with_na
mode_imputed[is.na(mode_imputed)] <- get_mode(mode_imputed)
# 5. 使用mice包进行多重插补
library(mice)
# 创建含多个变量的数据框
df <- data.frame(
x = c(1, 2, NA, 4, 5),
y = c(NA, 2, 3, 4, 5),
z = c(1, NA, 3, 4, NA)
)
# 执行多重插补
imp <- mice(df, m = 5, maxit = 50, method = 'pmm', seed = 123)
# 查看插补结果
print(imp)
# 提取完整数据
complete_data <- complete(imp, 1)
处理策略选择
：
删除法
：适用于缺失比例很小（%）且随机缺失的情况
填充法
：适用于缺失比例中等（5-20%）的情况
多重插补
：适用于缺失比例较大（>20%）或需要保留数据不确定性的情况
3.2 异常值检测与处理
异常值检测方法
# 创建含异常值的数据
set.seed(123)
normal_data <- rnorm(100, mean = 0, sd = 1)
data_with_outliers <- c(normal_data, 5, -5, 8)
# 1. 基于标准差的方法
mean_val <- mean(data_with_outliers)
sd_val <- sd(data_with_outliers)
outliers_sd <- data_with_outliers[abs(data_with_outliers - mean_val) > 3 * sd_val]
# 2. 基于四分位距（IQR）的方法
Q1 <- quantile(data_with_outliers, 0.25)
Q3 <- quantile(data_with_outliers, 0.75)
IQR <- Q3 - Q1
lower_bound <- Q1 - 1.5 * IQR
upper_bound <- Q3 + 1.5 * IQR
outliers_iqr <- data_with_outliers[data_with_outliers < lower_bound | data_with_outliers > upper_bound]
# 3. 使用箱线图检测
boxplot(data_with_outliers, main = “箱线图检测异常值”)
# 4. 使用LOF算法（局部异常因子）
library(DMwR2)
lof_scores <- lof(data_with_outliers, k = 5)
outliers_lof <- data_with_outliers[lof_scores > 1.5] # 阈值可调整
# 5. 使用孤立森林
library(isotree)
iso_model <- isolation.forest(data_with_outliers)
outliers_iso <- data_with_outliers[iso_model$scores > 0.6] # 阈值可调整
异常值处理
# 1. 删除异常值
clean_data <- data_with_outliers[!data_with_outliers %in% outliers_iqr]
# 2. Winsorization（缩尾处理）
winsorize <- function(x, lower = 0.05, upper = 0.95) {
q <- quantile(x, c(lower, upper), na.rm = TRUE)
x[x < q[1]] <- q[1]
x[x > q[2]] <- q[2]
return(x)
}
winsorized_data <- winsorize(data_with_outliers)
# 3. 使用中位数替换
median_replaced <- data_with_outliers
median_replaced[median_replaced < lower_bound | median_replaced > upper_bound] <- median(data_with_outliers)
3.3 数据标准化与归一化
标准化（Z-score）
# 创建数据
raw_data <- matrix(rnorm(100*5), nrow = 100, ncol = 5)
# 手动标准化
standardized_manual <- scale(raw_data)
# 使用caret包
library(caret)
preprocess_params <- preProcess(raw_data, method = c("center", "scale"))
standardized_caret <- predict(preprocess_params, raw_data)
# 验证
colMeans(standardized_caret) # 均值应为0
apply(standardized_caret, 2, sd) # 标准差应为1
归一化（Min-Max）
# 手动归一化
normalize <- function(x) {
return((x - min(x)) / (max(x) - min(x)))
}
normalized_manual <- apply(raw_data, 2, normalize)
# 使用caret包
preprocess_params <- preProcess(raw_data, method = "range")
normalized_caret <- predict(preprocess_params, raw_data)
3.4 大数据处理技巧
数据处理优化
# 1. 使用data.table包处理大数据
library(data.table)
# 创建大型数据框
large_df <- data.frame(
id = 1:1e6,
group = sample(LETTERS[1:20], 1e6, replace = TRUE),
value = rnorm(1e6)
)
# 转换为data.table
dt <- as.data.table(large_df)
# 高效分组统计
system.time({
result <- dt[, .(mean_value = mean(value),
sd_value = sd(value)),
by = group]
})
# 2. 使用dplyr包
library(dplyr)
system.time({
result <- large_df %>%
group_by(group) %>%
summarise(mean_value = mean(value),
sd_value = sd(value))
})
# 3. 并行计算
library(parallel)
# 创建任务函数
process_chunk <- function(chunk) {
return(mean(chunk$value))
}
# 分割数据并并行处理
chunks <- split(large_df, cut(1:nrow(large_df), 8))
results <- mclapply(chunks, process_chunk, mc.cores = 4)
内存管理技巧
# 1. 查看对象大小
object.size(large_df)
# 2. 删除大对象释放内存
rm(large_df)
gc() # 强制垃圾回收
# 2. 使用ff包处理超出内存的数据
library(ff)
# 创建ffdf对象
ffdf <- as.ffdf(large_df)
# 操作方式与普通data.frame类似，但数据存储在磁盘上
3.5 结果可视化与报告
自定义统计报告函数
# 创建综合统计报告函数
create_statistical_report <- function(data, group_var = NULL) {
# 加载必要包
if (!require(psych)) install.packages("psych")
library(psych)
report <- list()
# 基本描述性统计
if (is.null(group_var)) {
report$descriptive <- describe(data)
} else {
# 分组描述性统计
report$descriptive <- describeBy(data, group_var, mat = TRUE)
}
# 相关性分析
if (ncol(data) > 1) {
report$correlation <- cor(data, use = "complete.obs")
}
# 缺失值统计
report$missing <- sapply(data, function(x) sum(is.na(x)))
# 异常值统计（基于IQR）
find_outliers_iqr <- function(x) {
if (!is.numeric(x)) return(0)
Q1 <- quantile(x, 0.25, na.rm = TRUE)
Q3 <- quantile(x, 0.75, na.rm = TRUE)
IQR <- Q3 - Q1
lower <- Q1 - 1.5 * IQR
upper <- Q3 + 1.5 * IQR
sum(x < lower | x > upper, na.rm = TRUE)
}
report$outliers <- sapply(data, find_outliers_iqr)
# 返回报告
return(report)
}
# 使用示例
sample_data <- data.frame(
age = c(25, 30, 35, NA, 40, 100, 45, 50),
income = c(50000, 60000, 65000, 70000, 75000, 80000, 85000, 90000),
score = c(80, 85, 90, 95, 100, 105, 110, 115)
)
report <- create_statistical_report(sample_data)
print(report)
四、高级技巧与最佳实践
4.1 函数式编程在统计中的应用
R语言是函数式编程语言，利用这一特性可以编写更高效的统计代码。
# 使用lapply进行批量统计计算
calculate_stats <- function(data_list) {
lapply(data_list, function(x) {
if (is.numeric(x)) {
c(mean = mean(x, na.rm = TRUE),
sd = sd(x, na.rm = TRUE),
n = sum(!is.na(x)))
} else {
NULL
}
})
}
# 应用示例
data_list <- list(
var1 = rnorm(100),
var2 = rnorm(100),
var3 = sample(c("A", "B"), 100, replace = TRUE)
)
stats <- calculate_stats(data_list)
4.2 使用purrr包进行高级函数式编程
library(purrr)
# 创建统计函数列表
stat_functions <- list(
mean = ~mean(.x, na.rm = TRUE),
median = ~median(.x, na.rm = TRUE),
sd = ~sd(.x, na.rm = TRUE),
n = ~sum(!is.na(.x))
)
# 批量应用到数据框的每一列
data <- data.frame(
x = rnorm(100),
y = rnorm(100),
z = rnorm(100)
)
# 使用map_dfr将结果合并为数据框
stats_summary <- map_dfr(names(data), function(col) {
map_dfr(stat_functions, function(fn) {
tibble(value = fn(data[[col]]))
}, .id = "stat") %>%
mutate(variable = col)
}) %>%
pivot_wider(names_from = stat, values_from = value)
print(stats_summary)
4.3 可重现研究的最佳实践
# 1. 使用set.seed()确保随机过程可重现
set.seed(123)
# 2. 使用sessionInfo()记录环境信息
sessionInfo()
# 3. 使用renv包管理项目依赖
# install.packages(“renv”)
# renv::init() # 初始化项目环境
# renv::snapshot() # 保存当前包版本
# 4. 使用targets包管理分析流程
library(targets)
# 定义分析流程
list(
tar_target(data, rnorm(100)),
tar_target(stats, list(mean = mean(data), sd = sd(data))),
tar_target(plot, plot(data))
)
# 5. 使用knitr/R Markdown生成报告
# 在R Markdown中嵌入统计分析代码
4.4 性能优化技巧
# 1. 向量化操作优于循环
# 慢的循环方式
slow_function <- function(n) {
result <- numeric(n)
for (i in 1:n) {
result[i] <- i^2
}
return(result)
}
# 快的向量化方式
fast_function <- function(n) {
return((1:n)^2)
}
# 性能对比
system.time(slow_function(1e6))
system.time(fast_function(1e6))
# 2. 使用Rcpp编写C++扩展
library(Rcpp)
cppFunction('
NumericVector fast_sum(NumericVector x) {
int n = x.size();
NumericVector result(n);
double total = 0;
for (int i = 0; i < n; ++i) {
total += x[i];
result[i] = total;
}
return result;
}
')
# 3. 使用data.table进行大数据处理
library(data.table)
dt <- data.table(x = rnorm(1e6), y = rnorm(1语言的统计函数进行深入分析。接下来，我们将探讨一些高级统计方法和实际应用技巧。
## 五、高级统计方法详解
### 5.1 贝叶斯统计分析
贝叶斯方法在现代统计分析中越来越重要，R语言提供了强大的贝叶斯分析工具。
#### 基础贝叶斯推断
```r
# 安装和加载必要的包
if (!requireNamespace("rstan", quietly = TRUE)) {
install.packages("rstan")
}
library(rstan)
# 贝叶斯正态分布参数估计
# 假设我们有以下数据
data <- c(5.2, 5.1, 4.9, 5.0, 5.3, 5.2, 4.8, 5.1, 5.0, 5.2)
# 编写Stan模型
stan_code <- "
data {
int
vector[N] y;
}
parameters {
real mu;
real
}
model {
// 先验分布
mu ~ normal(5, 0.5);
sigma ~ exponential(2);
// 似然函数
y ~ normal(mu, sigma);
}
”
# 准备数据
stan_data <- list(N = length(data), y = data)
# 运行模型
fit <- stan(model_code = stan_code, data = stan_data,
iter = 2000, chains = 4)
# 查看结果
print(fit)
plot(fit)
使用brms包进行贝叶斯回归
# 安装和加载brms
if (!requireNamespace("brms", quietly = TRUE)) {
install.packages("brms")
}
library(brms)
# 贝叶斯线性回归
# 创建数据
set.seed(123)
n <- 100
x <- rnorm(n)
y <- 2 + 3*x + rnorm(n, 0, 1)
# 拟合模型
bayes_lm <- brm(y ~ x, data = data.frame(x, y))
# 查看结果
summary(bayes_lm)
# 可视化后验分布
plot(bayes_lm)
# 预测
predictions <- predict(bayes_lm, newdata = data.frame(x = c(0, 1, 2)))
print(predictions)
5.2 生存分析
生存分析在医学研究、工程可靠性分析中应用广泛。
# 安装和加载survival包
if (!requireNamespace("survival", quietly = TRUE)) {
install.packages("survival")
}
library(survival)
# 创建生存数据
set.seed(123)
n <- 200
time <- rexp(n, rate = 0.1)
status <- rbinom(n, 1, 0.8) # 80%的事件发生率
treatment <- sample(c("A", "B"), n, replace = TRUE)
age <- rnorm(n, 50, 10)
surv_data <- data.frame(time, status, treatment, age)
# 创建生存对象
surv_obj <- Surv(time = surv_data$time, event = surv_data$status)
# Kaplan-Meier生存曲线
km_fit <- survfit(surv_obj ~ treatment, data = surv_data)
# 绘制生存曲线
plot(km_fit, col = c("blue", "red"), lwd = 2,
xlab = "Time", ylab = "Survival Probability",
main = "Kaplan-Meier Survival Curves")
legend("topright", legend = c("Treatment A", "Treatment B"),
col = c("blue", "red"), lwd = 2)
# Log-rank检验
survdiff(surv_obj ~ treatment, data = surv_data)
# Cox比例风险模型
cox_model <- coxph(surv_obj ~ treatment + age, data = surv_data)
summary(cox_model)
# 检验比例风险假设
cox.zph(cox_model)
5.3 混合效应模型
混合效应模型用于处理具有层次结构的数据。
# 安装和加载lme4包
if (!requireNamespace("lme4", quietly = TRUE)) {
install.packages("lme4")
}
library(lme4)
# 创建层次数据
set.seed(123)
n_schools <- 20
n_students_per_school <- 15
school_id <- rep(1:n_schools, each = n_students_per_school)
student_id <- 1:(n_schools * n_students_per_school)
treatment <- rep(sample(0:1, n_schools, replace = TRUE), each = n_students_per_school)
ability <- rnorm(n_schools * n_students_per_school)
# 随机截距：不同学校的基线水平不同
school_effect <- rnorm(n_schools, 0, 2)
score <- 50 + 5*treatment + 2*ability + school_effect[school_id] + rnorm(n_schools * n_students_per_school, 0, 3)
hierarchical_data <- data.frame(
school_id = factor(school_id),
student_id = factor(student_id),
treatment = factor(treatment),
ability = ability,
score = score
)
# 随机截距模型
ri_model <- lmer(score ~ treatment + ability + (1 | school_id),
data = hierarchical_data)
summary(ri_model)
# 随机斜率和截距模型
rs_model <- lmer(score ~ treatment + ability + (treatment | school_id),
data = hierarchical_data)
summary(rs_model)
# 模型比较
anova(ri_model, rs_model)
# 提取固定效应和随机效应
fixef(ri_model)
ranef(ri_model)
5.4 结构方程模型（SEM）
# 安装和加载lavaan包
if (!requireNamespace("lavaan", quietly = TRUE)) {
install.packages("lavaan")
}
library(lavaan)
# 创建模拟数据
set.seed(123)
n <- 500
# 潜变量：感知价值、满意度、忠诚度
perceived_value <- rnorm(n)
satisfaction <- 0.6*perceived_value + rnorm(n, 0, 0.8)
loyalty <- 0.5*satisfaction + 0.3*perceived_value + rnorm(n, 0, 0.7)
# 观测变量
x1 <- 0.7*perceived_value + rnorm(n, 0, 0.7)
x2 <- 0.8*perceived_value + rnorm(n, 0, 0.6)
x3 <- 0.75*perceived_value + rnorm(n, 0, 0.65)
y1 <- 0.8*satisfaction + rnorm(n, 0, 0.6)
y2 <- 0.85*satisfaction + rnorm(n, 0, 0.55)
z1 <- 0.9*loyalty + rnorm(n, 0, 0.45)
z2 <- 0.8*loyalty + rnorm(n, 0, 0.6)
sem_data <- data.frame(x1, x2, x3, y1, y2, z1, z2)
# 定义模型
model <- '
# 测量模型
PV =~ x1 + x2 + x3
SAT =~ y1 + y2
LOY =~ z1 + z2
# 结构模型
SAT ~ PV
LOY ~ SAT + PV
'
# 拟合模型
fit <- sem(model, data = sem_data)
# 查看结果
summary(fit, standardized = TRUE, fit.measures = TRUE)
# 可视化
library(semPlot)
semPaths(fit, what = "std", layout = "tree",
edge.label.cex = 0.8, sizeMan = 5)
六、实际应用案例
6.1 客户流失预测分析
# 模拟电信客户数据
set.seed(123)
n <- 1000
customer_data <- data.frame(
customer_id = 1:n,
tenure = sample(1:72, n, replace = TRUE),
monthly_charges = runif(n, 30, 120),
total_charges = tenure * monthly_charges,
contract = sample(c("Month-to-month", "One year", "Two year"), n,
prob = c(0.6, 0.25, 0.15), replace = TRUE),
internet_service = sample(c("DSL", "Fiber optic", "No"), n,
prob = c(0.3, 0.4, 0.3), replace = TRUE),
payment_method = sample(c("Electronic check", "Mailed check",
"Bank transfer", "Credit card"), n, replace = TRUE),
churn = rbinom(n, 1, 0.267) # 约27%的流失率
)
# 数据预处理
customer_data$contract <- factor(customer_data$contract)
customer_data$internet_service <- factor(customer_data$internet_service)
customer_data$payment_method <- factor(customer_data$payment_method)
# 探索性分析
library(ggplot2)
# 流失率按合同类型
ggplot(customer_data, aes(x = contract, fill = factor(churn))) +
geom_bar(position = "fill") +
labs(title = "流失率按合同类型", y = "比例")
# 逻辑回归模型
churn_model <- glm(churn ~ tenure + monthly_charges + total_charges +
contract + internet_service + payment_method,
data = customer_data, family = binomial)
summary(churn_model)
# 预测概率
customer_data$predicted_prob <- predict(churn_model, type = "response")
# 评估模型
# 计算ROC曲线
library(pROC)
roc_obj <- roc(customer_data$churn, customer_data$predicted_prob)
plot(roc_obj, main = "ROC Curve")
auc(roc_obj) # AUC值
# 混淆矩阵
predicted_class <- ifelse(customer_data$predicted_prob > 0.5, 1, 0)
confusionMatrix <- table(Actual = customer_data$churn, Predicted = predicted_class)
print(confusionMatrix)
# 计算准确率、召回率等指标
accuracy <- sum(diag(confusionMatrix)) / sum(confusionMatrix)
precision <- confusionMatrix[2, 2] / sum(confusionMatrix[, 2])
recall <- confusionMatrix[2, 2] / sum(confusionMatrix[2, ])
cat("Accuracy:", accuracy, "\n")
cat("Precision:", precision, "\n")
cat("Recall:", recall, "\n")
6.2 A/B测试分析
# 模拟A/B测试数据
set.seed(123)
n_A <- 1000
n_B <- 1000
# 控制组A：转化率5%
conversion_A <- rbinom(n_A, 1, 0.05)
# 实验组B：转化率7%
conversion_B <- rbinom(n_B, 1, 0.07)
# 基本统计
conversion_rate_A <- mean(conversion_A)
conversion_rate_B <- mean(conversion_B)
cat("转化率A:", conversion_rate_A, "\n")
cat("转化率B:", conversion_rate_B, "\n")
# 比例检验
prop_test <- prop.test(x = c(sum(conversion_A), sum(conversion_B)),
n = c(n_A, n_B))
print(prop_test)
# 计算统计功效
power.prop.test(n = n_A, p1 = 0.05, p2 = 0.07, sig.level = 0.05)
# 可视化结果
library(ggplot2)
ab_test_results <- data.frame(
group = c("A", "B"),
conversion_rate = c(conversion_rate_A, conversion_rate_B),
lower_ci = c(prop.test(c(sum(conversion_A), sum(conversion_A)), c(n_A, n_A))$conf.int[1],
prop.test(c(sum(conversion_B), sum(conversion_B)), c(n_B, n_B))$conf.int[1]),
upper_ci = c(prop.test(c(sum(conversion_A), sum(conversion_A)), c(n_A, n_A))$conf.int[2],
prop.test(c(sum(conversion_B), sum(conversion_B)), c(n_B, n_B))$conf.int[2])
)
ggplot(ab_test_results, aes(x = group, y = conversion_rate, fill = group)) +
geom_col() +
geom_errorbar(aes(ymin = lower_ci, ymax = upper_ci), width = 0.2) +
labs(title = "A/B测试结果", y = "转化率", x = "组别") +
theme_minimal()
6.3 时间序列预测：销售预测
# 模拟销售数据
set.seed(123)
n <- 120 # 10年月度数据
# 基础趋势
trend <- 100 + 1:n * 0.5
# 季节性
seasonal <- 20 * sin(2 * pi * (1:n) / 12)
# 随机噪声
noise <- rnorm(n, 0, 5)
sales <- trend + seasonal + noise
# 转换为时间序列对象
sales_ts <- ts(sales, start = c(2014, 1), frequency = 12)
# 分解时间序列
decomp <- decompose(sales_ts)
plot(decomp)
# ARIMA预测
library(forecast)
# 自动选择ARIMA模型
auto_arima <- auto.arima(sales_ts)
summary(auto_arima)
# 预测未来12个月
forecast_result <- forecast(auto_arima, h = 12)
# 可视化
plot(forecast_result, main = "未来12个月销售预测", xlab = "时间", ylab = "销售额")
# 模型诊断
checkresiduals(auto_arima)
# 精确评估
# 将数据分为训练集和测试集
train <- window(sales_ts, end = c(2023, 6))
test <- window(sales_ts, start = c(2023, 7))
# 在训练集上拟合
model <- auto.arima(train)
# 预测测试集
pred <- forecast(model, h = length(test))
# 计算误差指标
rmse <- sqrt(mean((test - pred$mean)^2))
mae <- mean(abs(test - pred$mean))
mape <- mean(abs((test - pred$mean) / test)) * 100
cat("RMSE:", rmse, "\n")
cat("MAE:", mae, "\n")
cat("MAPE:", mape, "%\n")
七、性能优化与并行计算
7.1 代码性能分析
# 使用profvis进行性能分析
if (!requireNamespace("profvis", quietly = TRUE)) {
install.packages("profvis")
}
library(profvis)
# 创建一个性能较差的函数
slow_function <- function(n) {
result <- numeric(0)
for (i in 1:n) {
result <- c(result, i^2) # 每次迭代都扩展向量，效率低
}
return(result)
}
# 分析性能
profvis({
slow_function(10000)
})
7.2 并行计算
# 并行计算示例：蒙特卡洛模拟
library(parallel)
# 检测核心数
num_cores <- detectCores() - 1 # 保留一个核心给系统
# 串行版本
monte_carlo_serial <- function(n) {
results <- numeric(n)
for (i in 1:n) {
x <- runif(1)
y <- runif(1)
results[i] <- ifelse(x^2 + y^2 <= 1, 1, 0)
}
return(mean(results) * 4) # 估算π
}
# 并行版本
monte_carlo_parallel <- function(n, cores = num_cores) {
# 将任务分割到多个核心
results <- mclapply(1:cores, function(i) {
sub_n <- ceiling(n / cores)
count <- 0
for (j in 1:sub_n) {
x <- runif(1)
y <- runif(1)
if (x^2 + y^2 <= 1) count <- count + 1
}
return(count)
}, mc.cores = cores)
total_count <- sum(unlist(results))
return((total_count / n) * 4)
}
# 性能对比
n <- 1e6
system.time(serial_result <- monte_carlo_serial(n))
system.time(parallel_result <- monte_carlo_parallel(n))
cat("串行结果:", serial_result, "\n")
cat("并行结果:", parallel_result, "\n")
7.3 使用Rcpp加速计算
# 安装和加载Rcpp
if (!requireNamespace("Rcpp", quietly = TRUE)) {
install.packages("Rcpp")
}
library(Rcpp)
# 编写C++函数
cppFunction('
NumericVector fast_cumsum(NumericVector x) {
int n = x.size();
NumericVector result(n);
double total = 0;
for (int i = 0; i < n; ++i) {
total += x[i];
result[i] = total;
}
return result;
}
')
# 比较性能
x <- rnorm(1e6)
system.time(r_result <- cumsum(x))
system.time(cpp_result <- fast_cumsum(x))
# 验证结果是否相同
all.equal(r_result, cpp_result)
八、总结与建议
8.1 学习路径建议
基础阶段
：掌握基础统计函数、数据操作和可视化
中级阶段
：学习回归分析、方差分析、时间序列等标准统计方法
高级阶段
：掌握贝叶斯统计、混合效应模型、SEM等复杂方法
专家阶段
：精通性能优化、并行计算、可重现研究等高级技能
8.2 最佳实践清单
代码规范
：使用有意义的变量名，添加注释，保持函数单一职责
数据验证
：在分析前检查数据质量（缺失值、异常值）
版本控制
：使用Git管理代码版本
文档化
：使用R Markdown记录分析过程和结果
测试
：编写单元测试确保函数正确性
性能意识
：对大数据集考虑内存使用和计算效率
8.3 推荐资源
书籍
：《R语言实战》、《统计学习导论》、《贝叶斯数据分析》
在线课程
：Coursera上的R语言专项课程、DataCamp
社区
：Stack Overflow、RStudio Community、GitHub
包文档
：CRAN任务视图、包vignettes
通过掌握这些统计函数和技巧，你将能够高效地解决各种数据分析问题，从简单的描述性统计到复杂的模型构建，都能游刃有余。记住，优秀的数据分析师不仅要知道如何使用工具，更要理解统计方法背后的原理，这样才能在实际应用中做出正确的判断。