1、基本概念
1.1、PBRTQC
基于患者样本的实时质量控制(patient-based real-time quality control,PBRTQC)系统是一种基于统计学以及数学模型的质量控制方法,是根据患者样本检测结果,利用统计学模型建立的一套以实时监测实验室检测质量的模型或者规则,是医学检验实验室提高质量控制体系的重要发展方向,最早由Hoffmann和Waid于1965年提出
1.2、PBRTQC模型
- 五个关键参数:上截断值(UTL)、下截断值(LTL)、浮动窗口大小(N)、上控制限(UCL)、下控制限(LCL)
- 排除规则:对于偏态分布以及存在极端值的分布数据,应选择一个适宜的截断值区间,将在截断值区间外的数据进行缩尾或者去除
- SPC算法:浮动均值(MA)、浮动分位数(MQ)、指数加权移动平均值(EWMA)用于监测定值误差(CE)和百分比误差(PE);浮动标准差(MovSD)、浮动非正常值患者数(MovSO)用于监测随机误差(RE)
- 控制限:1. 分布法:通过分析历史数据生成的PBRTQC计算值,得出计算值的均值与标准差,根据正态分布或其它分布特征计算出对应的控制限
2. 百分位点法:根据目标假阳性报警率(DFAR),如0.1%,考虑双边的情况下,那么控制限则为历史数据PBRTQC计算值的0.05%和99.95%百分位数点
1.3、模型评估指标
- 假阳性报警率(FAR):指误报次数与样本数的比值,文献中通常要求该值≤0.1%,即每进行1000例患者样本的实时质控,最多允许1次误报,表征模型的特异性,FAR越小,模型特异性越好
- 误差检出前平均患者样本数(ANPed):误差出现时至模型检测出误差所经历的患者样本例数,表征模型的灵敏度,例数越少,模型灵敏度越高
- 二者关系:二者之间存在矛盾关系,FAR越低,ANPed越大,实际应用中应综合考虑,能满足当前项目的需求即可
2、实现简要流程
- 数据类型为信号量,函数模型修改为mq,mq分位数设置为50,转换方式修改为不转换
- 根据条件查询数据集合(仪器+批号+项目+时间)
- 定量定性项目过滤,定量项目过滤掉单位为S/CO,定性项目液体类型为C(校准)或单位为S/CO
- cut-off过滤(不包括校准点数据)
- 根据washStatus判断是否需要排除修改前后信号值不一致的数据
- 得到两份数据。allDataList:经过上述步骤剩余的数据。list:经过上述步骤再过滤掉校准点(液体类型C)数据
- list为空、list大于20000,list大小小于训练数据集样本容量n。直接返回
- 寻找定标线。液体类型为c(校准),连续12个点确定一个定标线(不连续的舍弃)
- 暂存list
- 遍历函数模型集合,一个模型出一个图
- 暂存全量数据list
- 计算控制线
- 每次移动训练数据集个数个数据。计算截断限LTL和UTL(通过移动分位数mq计算)
- 根据LTL和UTL截断sublist数据。数据被截断并且是定标点,将前一个点设置为定标点
- 数据转换sublist,使数据更符合正态性(0-Box-Cox变换、1-对数变换、2-倒数变换、3-平方根变换)
- 记录最优lambda,加权失控报警0.5N和2N最优lambda用N计算出的lambda
- 计算本次控制线(LCL,UCL)
- 判断控制线的合理性:根据控制线偏差和Moving-Slope斜率最大值,最小值判断(斜率判断最新版本已废除)
- 满足合理性判断时,暂存N转换后的数据,用作失控报警中的加权失控计算(0.5N,2N)
- 未找到合适控制线时,直接返回并提示用户
- 使用计算控制线得到的LTL和UTL截断全量数据,得到list
- list数据转换(0-Box-Cox变换、1-对数变换、2-倒数变换、3-平方根变换)
- 暂存转换后的数据,用作失控报警中的加权失控计算
- list统计值计算(0-MA、1-MQ、2-EWMA、3-MovSD、4-MR)
- Moving-Slope计算系统误差点(最新版本废弃)
- 失控报警:失控点计算(0-单点失控,1-连续X点失控,2-加权失控,3-加权连续失控)
- 找到所有定标线
- 响应封装
3、代码实现
3.1、基本流程框架
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128
| public AjaxResult pbrtqcLineData(PbrtqcDTO dto) { PBRTQCResultVO response = new PBRTQCResultVO(); if (dto.getDataType() == 1) { dto.setFuncModel(1); dto.setMqQuantile(50.0); dto.setTransType(null); } List<PbrtqcWashed> allDataList = iPbrtqcWashedService.selectLineData(dto); if (dto.getProjectType() == 0){ allDataList = allDataList.stream().filter(item -> !"S/CO".equals(item.getUnit())).collect(Collectors.toList()); }else { allDataList = allDataList.stream().filter(item -> "C".equals(item.getLiquidType()) || "S/CO".equals(item.getUnit())).collect(Collectors.toList()); } if (1 == dto.getCutOffCondition()){ allDataList = allDataList.stream().filter(item -> "C".equals(item.getLiquidType()) || item.getValue() < dto.getCutOffValue()).collect(Collectors.toList()); }else if (2 == dto.getCutOffCondition()){ allDataList = allDataList.stream().filter(item -> "C".equals(item.getLiquidType()) || BigDecimal.valueOf(item.getValue()).compareTo(BigDecimal.valueOf(dto.getCutOffValue())) >= 0).collect(Collectors.toList()); } if (dto.getWashStatus() == 0){ allDataList = allDataList.stream().filter(item -> "C".equals(item.getLiquidType()) || item.getModifiedSignal().equals(item.getOriginalSignal())).collect(Collectors.toList()); } List<PbrtqcWashed> list = allDataList.stream().filter(item -> !"C".equals(item.getLiquidType())).collect(Collectors.toList()); if (CollectionUtils.isEmpty(list)) { return AjaxResult.error("查询数据为空"); } if (list.size() > 20000) { return AjaxResult.error("数据量过大(超过20000条),请调整查询条件"); } if (dto.getN() > list.size()) { return AjaxResult.error("总数据量为:" + list.size() + ",小于训练数据集样本N,请调整N大小!"); } List<CalibrationPoint> calibrationDetailList = findCalibrationLine(allDataList); List<PbrtqcWashed> tempList = saveDataList(list); String unit = list.get(0).getUnit(); List<PBRTQCData> pbrtqcDataList = new ArrayList<>(); List<Integer> funcModels = Arrays.stream(dto.getFuncModels().split(",")).map(Integer::parseInt).collect(Collectors.toList()); for (int i=0; i< funcModels.size(); i++){ if (funcModels.get(i) == 3){ funcModels.remove(i); funcModels.add(3); } } dto.setIndex(0); for (Integer funcModel : funcModels) { List<PbrtqcWashed> modelDataList = saveDataList(list); dto.setFuncModel(funcModel); if (funcModel == 3) { dto.setLimitQuantile(dto.getSdLimitQuantile()); } else { dto.setLimitQuantile(dto.getMqLimitQuantile()); } ControlLineData controlLine = calculateControlLine(dto, tempList, 1, 0, dto.getIndex()); if ((dto.getIndex() == 0 && controlLine.getUcl() == null) || (dto.getIndex() > 0 && dto.getFindControlLine() == null)){ return AjaxResult.error("未找到合适的控制线,请设置合理的控制限偏差"); } List<PbrtqcWashed> cutDataList = modelDataList; if (dto.getCutType() != 0) { QuantileData quantileData = new QuantileData(dto.getLtl(),dto.getUtl()); cutDataList = cutData(dto, quantileData, modelDataList); } dataTrans(dto, cutDataList, 0); List<PbrtqcWashed> transData = new ArrayList<>(); if (dto.getAlarmMode() == 2 || dto.getAlarmMode() == 3){ transData = saveDataList(cutDataList); } List<PBRTQCLineVO> dataList = funcModelCalculate(dto, cutDataList,1); List<PBRTQCLineVO> systemErrorDataList = systemErrorPoint(dto, dataList); List<PBRTQCLineVO> alarmDataList = outOfControlAlarm(dto, tempList, dataList, controlLine, transData); List<String> lots = dataList.stream().map(PBRTQCLineVO::getName).distinct().collect(Collectors.toList()); lots.add(0, "失控点"); lots.add(1, "系统误差点"); alarmDataList.forEach(item -> dataList.get(item.getId()).setName("失控点")); List<PBRTQCLineVO> calibrationPointList = dataList.stream().filter(item -> item.getIsCalibration() != null && item.getIsCalibration().equals(1)).collect(Collectors.toList()); List<Integer> calibrationList = calibrationPointList.stream().map(PBRTQCLineVO::getId).collect(Collectors.toList()); List<PBRTQCLineVO> systemErrorData = alarmDataList.stream().filter(item -> systemErrorDataList.stream().anyMatch(obj -> obj.getId().equals(item.getId()))).collect(Collectors.toList()); systemErrorData.forEach(item -> dataList.get(item.getId()).setName("系统误差点")); PBRTQCData data = new PBRTQCData(); data.setList(dataList); data.setLots(lots); data.setUnit(unit); data.setMethodName(DictUtils.getDictLabel("ims_pbrtqc_func_model", dto.getFuncModel().toString())); data.setAlarmList(alarmDataList); data.setSystemErrorList(systemErrorData); data.setLcl(controlLine.getLcl()); data.setUcl(controlLine.getUcl()); data.setOriginLcl(controlLine.getOriginLcl()); data.setOriginUcl(controlLine.getOriginUcl()); data.setCalibrationList(calibrationList); pbrtqcDataList.add(data); dto.setIndex(dto.getIndex() + 1); } response.setList(pbrtqcDataList); response.setLtl(dto.getLtl()); response.setUtl(dto.getUtl()); response.setCalibrationPointList(calibrationDetailList); return AjaxResult.success(response); }
|
3.2、控制线计算流程
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88
| private ControlLineData calculateControlLine(PbrtqcDTO dto, List<PbrtqcWashed> list, Integer transFlag, Integer alarmFlag, Integer index) {
int size = list.size(); Integer offset = dto.getLimitOffset(); Integer n = dto.getN(); Double limitQuantile = dto.getLimitQuantile(); ControlLineData controlLineData = new ControlLineData(); if(index > 0){ ControlLineData controlLine = calculateCurrentControlLine(dto.getTransDataList(), dto); Double ucl = controlLine.getUcl(); Double lcl = controlLine.getLcl(); Double originLcl = controlLine.getOriginLcl(); Double originUcl = controlLine.getOriginUcl(); boolean msFlag = true; if (dto.getMinSlope() != null && dto.getMaxSlope() != null){ msFlag = movingSlope(dto, dto.getStatisticList2()); } if ((originUcl - originLcl) / (originLcl + originUcl) <= limitQuantile / 100 && msFlag) { dto.setFindControlLine(1); controlLineData.setLcl(lcl); controlLineData.setUcl(ucl); controlLineData.setOriginUcl(originUcl); controlLineData.setOriginLcl(originLcl); dto.setBestMean2(dto.getUseMean2()); dto.setBestSD2(dto.getUseSD2()); return controlLineData; } } for (int i = 0; i < size - offset; i++) { if (i > 100) { break; } List<PbrtqcWashed> subList = list.stream().skip(offset * i).limit(n).collect(Collectors.toList()); subList = saveDataList(subList); if (subList.size() < n) { return controlLineData; } QuantileData quantileData = new QuantileData(); if (dto.getCutType() != 0) { quantileData = calculateQuantile(subList, dto); subList = cutData(dto, quantileData, subList); } List<PbrtqcWashed> transDatas = dataTrans(dto, subList, transFlag); List<PbrtqcWashed> tempTransDataList = new ArrayList<>(); if (dto.getAlarmMode() == 2 || dto.getAlarmMode() == 3) { tempTransDataList = saveDataList(transDatas); } ControlLineData currentControlLine = calculateCurrentControlLine(transDatas, dto); Double ucl = currentControlLine.getUcl(); Double lcl = currentControlLine.getLcl(); Double originLcl = currentControlLine.getOriginLcl(); Double originUcl = currentControlLine.getOriginUcl(); boolean msFlag = true; if (dto.getMinSlope() != null && dto.getMaxSlope() != null){ msFlag = movingSlope(dto, dto.getStatisticList()); } if ((originUcl - originLcl) / (originLcl + originUcl) <= limitQuantile / 100 && msFlag) { if (index == 0){ dto.setLtl(quantileData.getLtl()); dto.setUtl(quantileData.getUtl()); dto.setTransDataList(transDatas); dto.setBestMean(dto.getUseMean()); dto.setBestSD(dto.getUseSD()); } controlLineData.setLcl(lcl); controlLineData.setUcl(ucl); controlLineData.setOriginUcl(originUcl); controlLineData.setOriginLcl(originLcl); dto.setTransDataN(tempTransDataList); break; } } return controlLineData; }
|
3.3、数据转换box-cox转换
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141
| private ControlLineData calculateControlLine(PbrtqcDTO dto, List<PbrtqcWashed> list, Integer transFlag, Integer alarmFlag, Integer index) {
int size = list.size(); Integer offset = dto.getLimitOffset(); Integer n = dto.getN(); Double limitQuantile = dto.getLimitQuantile(); ControlLineData controlLineData = new ControlLineData(); if(index > 0){ ControlLineData controlLine = calculateCurrentControlLine(dto.getTransDataList(), dto); Double ucl = controlLine.getUcl(); Double lcl = controlLine.getLcl(); Double originLcl = controlLine.getOriginLcl(); Double originUcl = controlLine.getOriginUcl(); boolean msFlag = true; if (dto.getMinSlope() != null && dto.getMaxSlope() != null){ msFlag = movingSlope(dto, dto.getStatisticList2()); } if ((originUcl - originLcl) / (originLcl + originUcl) <= limitQuantile / 100 && msFlag) { dto.setFindControlLine(1); controlLineData.setLcl(lcl); controlLineData.setUcl(ucl); controlLineData.setOriginUcl(originUcl); controlLineData.setOriginLcl(originLcl); dto.setBestMean2(dto.getUseMean2()); dto.setBestSD2(dto.getUseSD2()); return controlLineData; } } for (int i = 0; i < size - offset; i++) { if (i > 100) { break; } List<PbrtqcWashed> subList = list.stream().skip(offset * i).limit(n).collect(Collectors.toList()); subList = saveDataList(subList); if (subList.size() < n) { return controlLineData; } QuantileData quantileData = new QuantileData(); if (dto.getCutType() != 0) { quantileData = calculateQuantile(subList, dto); subList = cutData(dto, quantileData, subList); } List<PbrtqcWashed> transDatas = dataTrans(dto, subList, transFlag); List<PbrtqcWashed> tempTransDataList = new ArrayList<>(); if (dto.getAlarmMode() == 2 || dto.getAlarmMode() == 3) { tempTransDataList = saveDataList(transDatas); } ControlLineData currentControlLine = calculateCurrentControlLine(transDatas, dto); Double ucl = currentControlLine.getUcl(); Double lcl = currentControlLine.getLcl(); Double originLcl = currentControlLine.getOriginLcl(); Double originUcl = currentControlLine.getOriginUcl(); boolean msFlag = true; if (dto.getMinSlope() != null && dto.getMaxSlope() != null){ msFlag = movingSlope(dto, dto.getStatisticList()); } if ((originUcl - originLcl) / (originLcl + originUcl) <= limitQuantile / 100 && msFlag) { if (index == 0){ dto.setLtl(quantileData.getLtl()); dto.setUtl(quantileData.getUtl()); dto.setTransDataList(transDatas); dto.setBestMean(dto.getUseMean()); dto.setBestSD(dto.getUseSD()); } controlLineData.setLcl(lcl); controlLineData.setUcl(ucl); controlLineData.setOriginUcl(originUcl); controlLineData.setOriginLcl(originLcl); dto.setTransDataN(tempTransDataList); break; } } return controlLineData; }
public static double findBestLambda(double[] data) { double lambdaStart = -5.0; double lambdaEnd = 5.0; double lambdaStep = 0.01; double G = 1.0; for (double x : data) { G *= Math.pow(x, 1.0 / data.length); } double minSD = Double.MAX_VALUE; double minLambda = -5; for (double lambda = lambdaStart; lambda <= lambdaEnd; lambda += lambdaStep) { double[] transformedData = boxCoxTransform(data, lambda, G); double avg = calculateAverage(transformedData); double sd = calculateStandardDeviation(transformedData, avg); if (sd < minSD) { minSD = sd; minLambda = lambda; } } return minLambda; }
public static double[] boxCoxTransform(double[] data, double lambda, double G) { double[] transformedData = new double[data.length]; if (lambda == 0) { for (int i = 0; i < data.length; i++) { transformedData[i] = G * Math.log(data[i]); } } else { for (int i = 0; i < data.length; i++) { transformedData[i] = (Math.pow(data[i], lambda) - 1) / (lambda * Math.pow(G, lambda - 1)); } } return transformedData; }
|
3.4、统计值计算之mq
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35
|
public static List<PBRTQCLineVO> MQ(Integer N, Double X, List<PbrtqcWashed> list) { List<PBRTQCLineVO> result = new ArrayList<>(); for (int i = 0; i < list.size() -N + 1; i++) { double val1 = (N + 1) * X / 100; int j = (int) val1; double g = val1 - j; List<Integer> subIndex = new ArrayList<>(); double calValue; List<PbrtqcWashed> subList = list.stream().skip(i).limit(N).collect(Collectors.toList()); List<PbrtqcWashed> collect = subList.stream().sorted(Comparator.comparing(PbrtqcWashed::getValue)).collect(Collectors.toList()); if (g == 0) { calValue = collect.get(j - 1).getValue(); } else { calValue = (1 - g) * collect.get(j - 1).getValue() + g * collect.get(j).getValue(); } PBRTQCLineVO vo = new PBRTQCLineVO(); vo.setId(i); vo.setName(list.get(i+N-1).getReagentBatch()); vo.setValue(calValue); vo.setIsCalibration(list.get(i+N-1).getIsCalibration()); vo.setReadyTime(list.get(i+N-1).getReadyTime()); result.add(vo); } return result; }
|
4、效果图