四分格相关系数的标准差及显著性检验

本文详细解释了如何使用R语言进行四分格相关系数的显著性检验,包括计算标准差、求解相关系数及判断显著性等步骤,并提供完整代码实现。

摘要生成于 C知道 ,由 DeepSeek-R1 满血版支持, 前往体验 >

    http://www.john-uebersax.com/stat/tetra.htm 提到可用两种方法做四分格相关系数的显著性检验,这里用到的是第一种方法,也就是利用标准差的检验.

     <The Tetrachoric Correlation and its Asymptotic Standard Error>一文讲到四分格相关系数的标准差为s:

tetrachoric.se <- function(a, b, c, d) {
  n <- a + b + c + d
  z1 <- qnorm((a + c)/n)
  z2 <- qnorm((a + b)/n)
  r <- psych::tetrachoric(matrix(c(a, b, c, d), 2, 2))$rho
  f1 <- pnorm((z1 - r * z2)/sqrt(1 - r * r)) - 0.5
  f2 <- pnorm((z2 - r * z1)/sqrt(1 - r * r)) - 0.5
  t <- (a + d) * (b + c)/4 + 
    (a + c) * (b + d) * f2 * f2 + 
    (a + b) * (c + d) * f1 * f1 +
    2 * (a * d - b * c) * f1 * f2 +
    (c * d - b * a) * f2 +
    (b * d - a * c) * f1
  s <- n^(-3/2) * sqrt(t)/ 
    mnormt::dmnorm(c(z1, z2), varcov = matrix(c(1,  r,  r, 1), 2, 2))
  c(n, z1, z2, r, f1, f2, t, s)
}


    但是给出的假定相关系数为0时的标准差表达式有误,具体为下列代码中的s1,修正为s2:

tetrachoric.se0 <- function(a, b, c, d) {
  n <- a + b + c + d
  z1 <- qnorm((a + c)/n)
  z2 <- qnorm((a + b)/n)
  r <- 0
  f1 <- pnorm((z1 - r * z2)/sqrt(1 - r * r)) - 0.5
  f2 <- pnorm((z2 - r * z1)/sqrt(1 - r * r)) - 0.5
  t <- (a + d) * (b + c)/4 + 
    (a + c) * (b + d) * f2 * f2 + 
    (a + b) * (c + d) * f1 * f1 +
    2 * (a * d - b * c) * f1 * f2 +
    (c * d - b * a) * f2 +
    (b * d - a * c) * f1
  s <- n^(-3/2) * sqrt(t)/ 
    mnormt::dmnorm(c(z1, z2), varcov = matrix(c(1,  r,  r, 1), 2, 2))
  s1 <- n^(-5/2) * sqrt((a + b) * (a + c) * (b + d) * (b + c))/ 
    mnormt::dmnorm(c(z1, z2), varcov = matrix(c(1,  0,  0, 1), 2, 2))
  s2 <- n^(-5/2) * sqrt((a + b) * (a + c) * (b + d) * (d + c))/ 
    mnormt::dmnorm(c(z1, z2), varcov = matrix(c(1,  0,  0, 1), 2, 2))
  c(n, z1, z2, r, f1, f2, t, s, s1, s2)
}

    其中s为直接将相关系数设为0之后求得的标准差, s1为文章提供的表达式求得(估计印刷错误),s2为修正之后的值,与s的值相同.


    有了相关系数为0的标准差之后,就可以检验相关系数的显著性了,p值为:

2 * pnorm(abs(r)/s2, lower.tail = F)

    写一个完整的求四分格相关系数显著性的R语言代码:

tetrachoric.sig <- function(a, b, c, d) {
  n <- a + b + c + d
  z1 <- qnorm((a + c)/n)
  z2 <- qnorm((a + b)/n)
  r <- psych::tetrachoric(matrix(c(a, b, c, d), 2, 2))$rho
  se <- n^(-5/2) * sqrt((a + b) * (a + c) * (b + d) * (d + c))/ 
    mnormt::dmnorm(c(z1, z2), varcov = matrix(c(1,  0,  0, 1), 2, 2))
  p <- 2 * pnorm(abs(r)/se, lower.tail = F)
  list(rho = r, p.value = p)
}

tetrachoric.sig <- function(m) {
  n <- sum(m)
  z1 <- qnorm(sum(m[1, ])/n)
  z2 <- qnorm(sum(m[, 1])/n)
  r <- psych::tetrachoric(m)$rho
  se <- n^(-5/2) * sqrt(prod(rowSums(m), colSums(m)))/
  mnormt::dmnorm(c(z1, z2), varcov = matrix(c(1,  0,  0, 1), 2, 2))
  p <- 2 * pnorm(abs(r)/se, lower.tail = F)
  list(rho = r, p.value = p)
}


转载于:https://my.oschina.net/u/1791586/blog/336815

<think>我们面对的任务是基于NinaproDb5B类数据库,使用三种机器学习算法(SVM、随机森林、神经网络)对17种手势的sEMG信号进行分类。整个流程包括数据预处理、特征提取、模型训练与比较以及显著性分析。步骤1:数据预处理引用[1]中提到数据处理包括滤波、归一化等步骤。具体到我们的任务:a)滤波:通常使用带通滤波(如20-450Hz)去除噪声和工频干扰。b)活动段提取:由于EMG信号在动作开始和结束时有静息段,我们需要检测动作的开始和结束。可以使用能量阈值法(例如计算信号的均方根(RMS)然后设定阈值)来提取活动段。c)滑动窗口分割:为了获得足够的数据样本,我们使用滑动窗口(例如窗口长度200ms,重叠率50%)将活动段分割成多个样本段。步骤2:特征提取根据任务要求,我们需要提取三个特征:均方根值(RMS)、绝对平均值(MAV)和小波系数(WaveletCoefficients)。a)RMS:计算窗口内信号的均方根值。b)MAV:计算窗口内信号的绝对平均值。c)小波系数:通常使用小波变换(如Daubechies小波)分解信号,然后提取特定尺度的小波系数作为特征。我们可能需要选择几个尺度的系数(例如第3层到第7层的细节系数)并计算它们的能量或标准差等统计量。步骤3:模型训练与比较我们将使用三种算法:支持向量机(SVM)、随机森林(RandomForest)和神经网络(NeuralNetwork)。a)数据划分:将数据集按照一定比例(如70%训练,30%测试)划分,并注意保持类别平衡(使用分层抽样)。b)模型训练:-SVM:选择核函数(如RBF),并通过交叉验证调整C和gamma参数。-随机森林:调整树的数量(n_estimators)和最大深度(max_depth)。-神经网络:设计一个多层感知机(MLP),包括输入层(特征维度)、隐藏层(如两个隐藏层,每层128个神经元)和输出层(17个神经元,使用softmax激活)。使用ReLU激活函数和Adam优化器。c)评估指标:准确率(分类正确率)、混淆矩阵、ROC曲线(由于是多类问题,需要为每个类别绘制ROC曲线,或使用宏平均/微平均ROC曲线)。步骤4:结果可视化a)分类正确率柱状图:绘制三种算法在测试集上的准确率柱状图。b)混淆矩阵:对每种算法绘制混淆矩阵,观察各类别的分类情况。c)ROC曲线:由于类别较多(17类),我们可以选择几个代表性的类别绘制ROC曲线,或者绘制多类ROC曲线的平均(宏平均或微平均)。步骤5:显著性分析使用SPSS软件进行显著性分析(通常使用方差分析ANOVA)来比较三种算法的性能差异是否显著。具体步骤:a)由于我们使用相同的数据集,可以采用重复测量设计(每个样本被三种算法分类,得到三个准确率?注意:这里不是每个样本的准确率,而是整个测试集的准确率。因此,我们需要采用多次重复实验(如10次)得到多个准确率值,然后进行配对t检验或ANOVA)。但实际上,我们通常会在相同的测试集上运行三种算法,得到三个准确率(每个算法一个值),但这样只有三个值无法进行显著性检验。因此,我们需要对每个算法进行多次运行(例如使用不同的随机种子,或者使用交叉验证)来获得多个准确率值(比如10次10折交叉验证,每次得到10个准确率,然后取平均作为该次的结果,这样每个算法有10个准确率值)。然后,我们可以用这三个样本(每个算法10个值)进行单因素方差分析(One-wayANOVA),如果结果显著,则进行事后检验(如TukeyHSD)比较两两之间的差异。但是,任务要求使用SPSS,因此我们可以将数据整理成如下式:算法(分类变量,3水平:SVM、RF、NN)和准确率(连续变量)。然后进行单因素方差分析。然而,由于我们使用相同的数据集,多次运行的结果可能高度相关,因此更合适的方法是使用重复测量ANOVA(将每个算法的多次运行作为重复测量)。或者,我们也可以使用非参数检验(如Friedman检验)来比较多个相关样本。考虑到实际操作的可行性,我们可以在整个流程中采用10次10折交叉验证,这样每个算法得到100个准确率值(10折*10次),然后我们可以对每个算法计算平均准确率。但注意,这100个值并不是独立的,所以我们可以将每次10折交叉验证的平均准确率作为一次运行结果(即每个算法有10个值,对应10次重复),然后进行方差分析。步骤6:代码实现我们将使用Python进行数据处理、特征提取和模型训练。SPSS部分需要将结果导出为CSV,然后在SPSS中进行分析。现在,我们按照步骤详细展开:1.数据预处理-滤波:使用4阶巴特沃斯带通滤波器(20-450Hz)-活动段提取:计算整个记录的RMS,然后设定阈值(例如最大RMS的20%),连续超过阈值一定时间(如300ms)认为是动作开始,低于阈值一定时间认为是动作结束。-滑动窗口:窗口长度200ms,重叠100ms(即50%重叠)。2.特征提取-RMS:$RMS=\sqrt{\frac{1}{N}\sum_{i=1}^{N}x_i^2}$-MAV:$MAV=\frac{1}{N}\sum_{i=1}^{N}|x_i|$-小波系数:使用'db7'小波进行5层分解,然后提取第3层到第5层的细节系数(cD3,cD4,cD5)和近似系数(cA5)。对于每个系数,我们计算其能量(平方和)和标准差,这样每个小波分解层可以得到2个特征(能量和标准差),3层共6个特征。注意:每个通道都会提取这些特征,所以如果有8个通道,则每个窗口的特征数为:(RMS+MAV+6)*8=(2+6)*8=64个特征。3.模型训练-数据标准化:使用StandardScaler(均值为0,方差为1)对特征进行标准化。-SVM:使用sklearn的SVC,核函数选择RBF,使用网搜索调整C和gamma。-随机森林:使用RandomForestClassifier,调整n_estimators(例如100,200,500)和max_depth(例如10,20,None)。-神经网络:使用MLPClassifier,隐藏层设为两个,每层128个神经元,激活函数为ReLU,优化器为Adam,学习率0.001,最大迭代次数500。4.评估-准确率:使用sklearn的accuracy_score。-混淆矩阵:使用confusion_matrix,然后使用seaborn绘制热图。-ROC曲线:由于是多类,我们可以使用OneVsRest策略,然后计算每个类别的ROC曲线,再计算宏平均ROC曲线(将每个类别的ROC曲线平均)并绘制。5.显著性分析-进行10次10折交叉验证,每次记录10个准确率(10折),然后计算这10折的平均值作为该次运行的准确率。这样每个算法有10个准确率值(10次运行)。-将数据整理成CSV,包含三列:SVM的10个值,RF的10个值,NN的10个值。-在SPSS中,使用重复测量ANOVA(或非参数检验)比较三种算法的差异。现在,我们开始编写代码框架(由于实际数据需要下载和处理,这里只提供伪代码和关键步骤):注意:NinaproDB5B类数据库包含8个电极,17种手势动作。伪代码:步骤1:加载数据从NinaproDB5B类数据库中加载数据。每个数据文件包含一个受试者的记录,包括EMG信号和标签。我们需要加载多个受试者的数据?任务没有明确,通常我们可以按受试者独立分析,也可以合并。这里我们假设合并所有受试者。步骤2:预处理滤波:使用scipy.signal.butter设计带通滤波器,然后filtfilt进行零相位滤波。活动段提取:计算整个记录的RMS(使用滑动窗口,比如200ms窗口,100ms步长),然后设定阈值(例如整个记录RMS最大值的20%)。然后找到连续超过阈值的区域。滑动窗口分割:在活动段内,使用200ms窗口,100ms步长分割。步骤3:特征提取对每个窗口,计算每个通道的RMS和MAV。对每个通道进行小波分解(使用pywt.wavedec),得到cA5,cD5,cD4,cD3(5层分解,所以有cA5,cD5,cD4,cD3,cD2,cD1,但我们只需要3-5层,即cD3,cD4,cD5和cA5)。然后计算cD3,cD4,cD5,cA5的能量和标准差(注意:cA5是近似系数,cD是细节系数)。这样每个通道有4个系数,每个系数计算2个特征(能量和标准差),共8个特征,再加上RMS和MAV,每个通道共10个特征。8个通道共80个特征。步骤4:划分数据集将特征和标签合并,然后划分训练集和测试集(70%训练,30%测试)。注意使用分层抽样(stratify)。步骤5:训练模型分别训练SVM、随机森林和神经网络。步骤6:评估在测试集上计算准确率、绘制混淆矩阵和ROC曲线。步骤7:多次运行进行显著性分析为了得到多个准确率值,我们进行10次重复实验(每次重新划分训练集和测试集?注意:这样会引入随机性,但数据划分不同会导致结果不同)。或者,我们使用10次10折交叉验证,这样每次划分都不同,然后记录每次交叉验证的平均准确率(10折的平均),这样每个算法有10个值(10次运行)。步骤8:导出数据到SPSS将三种算法的10个准确率值保存为CSV文件,然后在SPSS中分析。代码实现的关键点:1.数据加载:需要知道NinaproDB5B类数据库的数据式。通常,数据是.mat文件,可以使用scipy.io.loadmat加载。2.预处理和特征提取需要高效,因为数据量可能很大。3.小波分解:注意信号长度。我们的窗口是200ms,采样频率可能是2000Hz,所以窗口有400个点。进行5层小波分解是可行的(因为2^5=32<400)。由于实际代码较长,这里只给出关键步骤的代码片段。我们首先写一个函数用于预处理和特征提取:注意:实际应用中,我们需要考虑计算效率,可以使用并行处理(如joblib)来加速特征提取。由于时间限制,我们只提供关键步骤的代码框架,具体实现需要根据数据调整。下面我们开始写代码:</think>### 基于Ninapro Db5 B类数据库的sEMG手势分类方案 #### 1. 数据预处理 **(1) 滤波处理** 使用4阶巴特沃斯带通滤波器(20-450Hz)去除噪声和工频干扰: ```python from scipy.signal import butter, filtfilt def bandpass_filter(signal, low=20, high=450, fs=2000, order=4): nyq = 0.5 * fs b, a = butter(order, [low/nyq, high/nyq], btype='band') return filtfilt(b, a, signal) ``` **(2) 活动段提取** 基于信号能量阈值法检测动作起止点: ```python import numpy as np def detect_activity(emg, threshold=0.2, min_duration=0.3): rms = np.sqrt(np.mean(emg**2, axis=0)) threshold_val = threshold * np.max(rms) activity = rms > threshold_val # 形态学处理去除短时波动 return activity_segments ``` **(3) 滑动窗口分割** 窗口长度200ms,重叠率50%: ```python def sliding_window(data, window_size=400, overlap=200): # 2000Hz采样率 segments = [] for start in range(0, len(data)-window_size, overlap): segments.append(data[start:start+window_size]) return np.array(segments) ``` #### 2. 特征提取 提取每个通道的时域和小波特征: ```python import pywt def extract_features(segment): features = [] # 时域特征 features.append(np.sqrt(np.mean(segment**2))) # RMS features.append(np.mean(np.abs(segment))) # MAV # 小波特征(db4小波5层分解) coeffs = pywt.wavedec(segment, 'db4', level=5) for c in coeffs[1:3]: # 取D3-D4细节系数 features.append(np.std(c)) # 小波系数标准差 return np.array(features) ``` #### 3. 机器学习模型 **(1) 支持向量机 (SVM)** ```python from sklearn.svm import SVC svm_model = SVC(kernel='rbf', C=10, gamma=0.01, probability=True) ``` **(2) 随机森林** ```python from sklearn.ensemble import RandomForestClassifier rf_model = RandomForestClassifier(n_estimators=200, max_depth=10) ``` **(3) 神经网络** ```python from tensorflow.keras.models import Sequential from tensorflow.keras.layers import Dense nn_model = Sequential([ Dense(128, activation='relu', input_shape=(n_features,)), Dense(64, activation='relu'), Dense(17, activation='softmax') ]) nn_model.compile(optimizer='adam', loss='sparse_categorical_crossentropy') ``` #### 4. 模型评估与可视化 **(1) 分类正确率柱状图** ```python import matplotlib.pyplot as plt accuracies = [svm_acc*100, rf_acc*100, nn_acc*100] plt.bar(['SVM', 'Random Forest', 'Neural Network'], accuracies) plt.ylabel('Accuracy (%)') plt.title('Gesture Classification Performance') ``` **(2) 混淆矩阵** ```python from sklearn.metrics import ConfusionMatrixDisplay ConfusionMatrixDisplay.from_predictions(y_true, y_pred, display_labels=gesture_labels, cmap='Blues') ``` **(3) ROC曲线** ```python from sklearn.metrics import RocCurveDisplay fig, ax = plt.subplots(figsize=(10,8)) for i in range(17): RocCurveDisplay.from_predictions(y_true_binary[i], y_prob[i], name=gesture_labels[i], ax=ax) ``` #### 5. 显著性分析 (SPSS) 1. 导出模型性能指标(准确率、F1-score等) 2. 在SPSS中选择:`分析 → 比较均值 → 重复测量方差分析` - 主体内因子:算法类型(3水平) - 主体内变量:各次交叉验证的准确率 3. 进行球形检验后,若p<0.05则使用Greenhouse-Geisser校正 4. 事后检验:Bonferroni校正的多重比较 #### 6. 关键注意事项 1. **数据划分**:采用分层5折交叉验证,确保每类手势比例一致 2. **特征标准化**:使用`StandardScaler`进行Z-score标准化 3. **类别不平衡**:对少数手势采用SMOTE过采样 4. **参数优化**:使用贝叶斯优化搜索最佳超参数 > **结果解读**:通过显著性分析(如重复测量ANOVA),可确定三种算法在分类性能上是否存在统计学显著差异(p<0.05)。若存在显著差异,进一步的事后检验可明确最优算法[^2]。 --- ### 相关问题 1. 如何选择适合sEMG信号分析的小波基函数? 2. 在Ninapro数据库处理中,如何解决不同受试者间的个体差异问题? 3. 除RMS和MAV外,哪些时频域特征对sEMG手势识别效果显著? 4. 如何设计端到端的深度学习架构直接处理原始sEMG信号? 5. 在显著性分析中,何时选用参数检验(如ANOVA)或非参数检验(如Friedman检验)? [^1]: Ninapro-EMG数据库为开发更加智能、响应性的辅助技术提供了强大的数据支撑。 [^2]: 手势动作分类结果为识别出的手势标签中所占比例最高的标签。
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值