R语言循环执行Tobit回归时censReg未生效问题排查

问题背景 #

我正在对多个数据集循环开展回归分析，分组维度为每种抗生素对应不同浓度。部分数据存在下限为3.912023的左删失，因此这类子集需要调用censReg()替代lm()建模。
但实际运行过程中，即便多个子集包含该删失下限值，所有回归的输出结果都是lm()的返回值，原代码如下：

model.summary<- tibble()

for(j in unique(myData$antibiotic)){
        abx<- filter(myData, antibiotic == j)
        for(i in unique(abx$concentration)){
            if (min(abx$datapoints) == 3.912023){
                conc<-filter(abx, concentration == i)
                model<-censReg(x ~ y, data = conc, left = 3.912023) 
                model.coef<-coef(model)[2] 
                model.list<-c(j,i ,model.coef)
                model.summary<-bind_rows(model.summary, model.list)
} else {
                conc<-filter(abx, concentration == i)
                model<-lm(x ~ y, data = conc)
                model.coef<-coef(model)[2] 
                model.list<-c(j,i, model.coef)
                model.summary<-bind_rows(model.summary, model.list)
    }   

     }
}

问题原因 #

原代码存在两处核心错误导致censReg()分支从未触发：

• 判断对象错误：条件判断使用的是整个抗生素分组abx的datapoints最小值，而非按浓度拆分后的子集最小值，只要同抗生素下有一个浓度分组不存在删失，整个抗生素下的所有浓度分组都会走lm()分支
• 浮点数精度问题：手动输入的近似值3.912023和实际计算得到的log(50)存在存储精度误差，直接用==做相等判断几乎不会触发条件

修复方案 #

调整判断逻辑为按浓度拆分后的子集做判断，同时直接使用log(50)作为判断阈值和删失边界，避免浮点近似误差，修复后代码如下：

model_summary <- tibble ()

for (j in unique(censored2$antibiotic)) {
        abxs <- filter(censored2, antibiotic == j)
        print(j)
        for (i in unique(abxs$MIC)) {
                print(i)
                MICs <- filter(abxs, MIC == i)
                if(min(MICs$logCFU) == log(50) ){
                         print("censored")
                        my_model <- censReg(logCFU ~ hours, data = MICs, left = log(50))
                        my_coef <- coef(my_model)[2]
                        my_list <- c(j,i,my_coef)
                        model_summary <- bind_rows(model_summary, my_list)
                 } else {
                         print("not censored")
                        my_model <- lm(logCFU ~ hours, data = MICs)
                        my_coef <- coef(my_model)[2]
                        my_list <- c(j,i,my_coef)
                        model_summary <- bind_rows(model_summary,my_list)
                 }
        }
}

内容的提问来源于stack exchange，提问作者EJuskewitz