MATLAB和R及Python病例对照分析

MATLAB病例对照分析

在MATLAB中，病例对照分析通常用于比较暴露组与非暴露组的结果，并计算风险比、优势比（Odds Ratio, OR）等统计量。以下是一个基于2x2列联表的病例对照分析步骤，包括OR计算、置信区间、卡方检验等内容。

1. 设置数据

假设我们有以下数据：

• exposed_cases: 暴露的病例数。
• not_exposed_cases: 未暴露的病例数。
• exposed_controls: 暴露的对照数。
• not_exposed_controls: 未暴露的对照数。

可以在MATLAB中定义如下：

% 定义病例和对照的数据
exposed_cases = 50;
not_exposed_cases = 30;
exposed_controls = 20;
not_exposed_controls = 100;

2. 构建列联表

构建一个2x2列联表来整理数据。

% 构建2x2列联表
data_matrix = [exposed_cases, not_exposed_cases; exposed_controls, not_exposed_controls];
disp('2x2列联表:');
disp(data_matrix);

3. 计算优势比（Odds Ratio）

优势比的计算公式为：

$$\text{OR}=\frac{(\text{exposed\_cases}\times \text{not\_exposed\_controls})}{(\text{not\_exposed\_cases}\times \text{exposed\_controls})}$$

在MATLAB中计算如下：

% 计算优势比（OR）
odds_ratio = (exposed_cases * not_exposed_controls) / (not_exposed_cases * exposed_controls);
disp(['优势比（OR）: ', num2str(odds_ratio)]);

4. 计算优势比的置信区间

优势比的95%置信区间可使用以下公式：

$$\ln (\text{OR})\pm Z\times \sqrt{\frac{1}{\text{exposed\_cases}}+\frac{1}{\text{not\_exposed\_cases}}+\frac{1}{\text{exposed\_controls}}+\frac{1}{\text{not\_exposed\_controls}}}$$

其中， $Z\approx 1.96$ 用于95%的置信水平。

% 计算95%置信区间
Z = 1.96;
ln_OR = log(odds_ratio);
SE_ln_OR = sqrt(1/exposed_cases + 1/not_exposed_cases + 1/exposed_controls + 1/not_exposed_controls);
CI_lower = exp(ln_OR - Z * SE_ln_OR);
CI_upper = exp(ln_OR + Z * SE_ln_OR);disp(['95% 置信区间: [', num2str(CI_lower), ', ', num2str(CI_upper), ']']);

5. 进行卡方检验

卡方检验用于检测暴露与结果之间是否存在显著关联。

MATLAB提供了chi2gof函数可以执行卡方检验，不过对于2x2列联表更直接的方法是手动计算卡方统计量。

% 手动计算卡方统计量
total = sum(data_matrix, 'all');  % 总样本数
row_totals = sum(data_matrix, 2); % 每行合计
col_totals = sum(data_matrix, 1); % 每列合计% 计算期望值
expected_matrix = (row_totals * col_totals) / total;% 卡方统计量
chi2_stat = sum((data_matrix - expected_matrix).^2 ./ expected_matrix, 'all');
p_value = 1 - chi2cdf(chi2_stat, 1); % 自由度为1disp(['卡方统计量: ', num2str(chi2_stat)]);
disp(['p值: ', num2str(p_value)]);

6. 可选：进行Logistic回归

如果想调整其他变量对结果的影响，可以使用Logistic回归模型。

假设有一个包含所有数据的表格data，其中有两个变量：outcome（结果，1代表病例，0代表对照），和exposure（暴露，1代表暴露，0代表未暴露）。

% 示例数据（假设 data 是一个表格）
outcome = [ones(exposed_cases, 1); ones(not_exposed_cases, 1); zeros(exposed_controls, 1); zeros(not_exposed_controls, 1)];
exposure = [ones(exposed_cases, 1); zeros(not_exposed_cases, 1); ones(exposed_controls, 1); zeros(not_exposed_controls, 1)];
data = table(outcome, exposure);% 使用MATLAB的Logistic回归模型
mdl = fitglm(data, 'outcome ~ exposure', 'Distribution', 'binomial');
disp(mdl);% 提取优势比
odds_ratio_logistic = exp(mdl.Coefficients.Estimate(2));
disp(['Logistic回归的优势比: ', num2str(odds_ratio_logistic)]);

总结

在MATLAB中，完成了病例对照分析的以下步骤：

1. 设置数据并构建列联表：用表格整理病例和对照数据。
2. 计算优势比（OR）：测量暴露与结果之间的关联。
3. 计算置信区间：估计OR的统计不确定性。
4. 进行卡方检验：检测是否存在显著关联。
5. Logistic回归（可选）：调整其他变量对结果的影响。

R病例对照分析

在R中进行病例对照分析可以使用简单的函数来计算优势比（Odds Ratio, OR）、置信区间以及卡方检验。下面是具体步骤，包括数据准备、OR和置信区间计算、卡方检验、以及Logistic回归等。

1. 设置数据

假设我们有以下数据：

• exposed_cases: 暴露的病例数。
• not_exposed_cases: 未暴露的病例数。
• exposed_controls: 暴露的对照数。
• not_exposed_controls: 未暴露的对照数。

在R中，可以定义这些数据并创建一个2x2列联表：

# 定义数据
exposed_cases <- 50
not_exposed_cases <- 30
exposed_controls <- 20
not_exposed_controls <- 100# 创建2x2列联表
table_data <- matrix(c(exposed_cases, not_exposed_cases, exposed_controls, not_exposed_controls),nrow = 2, byrow = TRUE,dimnames = list("Exposure" = c("Exposed", "Not Exposed"),"Outcome" = c("Cases", "Controls")))print("2x2列联表:")
print(table_data)

2. 计算优势比（Odds Ratio）

使用epitools包中的oddsratio函数或fmsb包中的oddsratio函数，可以快速计算OR和95%置信区间。

# 安装并加载epitools包（如果未安装）
if(!require(epitools)) install.packages("epitools")
library(epitools)# 计算优势比
or_result <- oddsratio(table_data)
print("优势比（OR）:")
print(or_result$measure)
print("95%置信区间:")
print(or_result$conf.int)

3. 手动计算OR和置信区间

如果不使用第三方包，可以手动计算OR和其95%置信区间：

$$\text{OR}=\frac{(\text{exposed\_cases}\times \text{not\_exposed\_controls})}{(\text{not\_exposed\_cases}\times \text{exposed\_controls})}$$

$$\ln (\text{OR})\pm Z\times \sqrt{\frac{1}{\text{exposed\_cases}}+\frac{1}{\text{not\_exposed\_cases}}+\frac{1}{\text{exposed\_controls}}+\frac{1}{\text{not\_exposed\_controls}}}$$

其中 $Z\approx 1.96$ 对应95%的置信水平。

# 手动计算优势比
odds_ratio <- (exposed_cases * not_exposed_controls) / (not_exposed_cases * exposed_controls)
print(paste("手动计算的优势比（OR）:", odds_ratio))# 计算95%置信区间
ln_or <- log(odds_ratio)
se_ln_or <- sqrt(1/exposed_cases + 1/not_exposed_cases + 1/exposed_controls + 1/not_exposed_controls)
ci_lower <- exp(ln_or - 1.96 * se_ln_or)
ci_upper <- exp(ln_or + 1.96 * se_ln_or)
print(paste("95% 置信区间: [", ci_lower, ", ", ci_upper, "]"))

4. 进行卡方检验

卡方检验可以用来检验暴露与结果之间是否存在显著性关联，可以使用chisq.test函数来实现。

# 卡方检验
chi_test <- chisq.test(table_data, correct = TRUE)  # Yates校正
print("卡方统计量:")
print(chi_test$statistic)
print("p值:")
print(chi_test$p.value)

5. 进行Logistic回归（可选）

如果需要调整其他变量，可以使用Logistic回归分析。假设我们有一个包含所有数据的数据框 data，其中包含两个变量：outcome（1表示病例，0表示对照）和exposure（1表示暴露，0表示未暴露）。

# 构建数据框
outcome <- c(rep(1, exposed_cases + not_exposed_cases), rep(0, exposed_controls + not_exposed_controls))
exposure <- c(rep(1, exposed_cases), rep(0, not_exposed_cases), rep(1, exposed_controls), rep(0, not_exposed_controls))
data <- data.frame(outcome = outcome, exposure = exposure)# Logistic回归
logit_model <- glm(outcome ~ exposure, data = data, family = binomial)
summary(logit_model)# 计算Logistic回归的优势比
or_logit <- exp(coef(logit_model)["exposure"])
print(paste("Logistic回归的优势比:", or_logit))

总结

在R中，我们完成了病例对照分析的以下步骤：

1. 设置数据并构建列联表：将病例和对照数据整理成2x2表格。
2. 计算优势比（OR）和置信区间：使用epitools包或者手动计算OR和其置信区间。
3. 进行卡方检验：用来检验暴露和结果之间的显著性关联。
4. Logistic回归（可选）：用来调整其他潜在的影响因素对结果的影响。

👉更新：亚图跨际