====== LRT似然比 ====== lmtest 似然比检验Likelihood Ratio Test 比较两个线性模型 似然比检验在统计学中到底是什么❓ 我们先看图中的英文定义，LRT是一种假设检验，用于判断某组样本更可能是来源于两种未知分布中的哪一种。常见的形式是： H0:θ=θ0 vs Ha:θ=θa 其中θ是某未知分布F的参数。θ0，θa是我们的两种估计，而LRT就是帮助我们判断哪一种估计更make sense。通过LRT的命名，我们就能知道他用的核心值是likelihood value，其test-statistic如图中定义。 如果样本更倾向于Fθ0(x)的分布，那λ会偏大还是偏小？ ——偏大！（为什么呢？）回忆一下likelihood value的含义：表示样本出现的可能性。 如果λ偏大，说明分子大分母小，那也就是Fθ0(x)下样本出现的概率 大于 Fθa(x)下样本出现的概率，那也就说明样本更倾向于θ0的参数估计，i.e. H0，所以最后会NOT Reject H0。 General总结来看，当样本量不够大，不能确保使用上述定理时，我们可以通过significance level α的值来反向计算c的取值。 这里需要同学们回忆一下significance level α的概念：“拒真”的概率，即P(Rejection H0|H0 is true)。所以我们可以以Fθ0(x)为真实分布，计算在该分布下的likelihood value，随后计算得到c的取值。 {{:yuxi:lrt.jpeg?400|}} library(lmtest) Loading required package: zoo Attaching package: 'zoo' The following objects are masked from 'package:base': as.Date, as.Date.numeric ## with data from Greene (1993): ## load data and compute lags data("USDistLag") usdl <- na.contiguous(cbind(USDistLag, lag(USDistLag, k = -1))) colnames(usdl) <- c("con", "gnp", "con1", "gnp1") fm1 <- lm(con ~ gnp + gnp1, data = usdl) fm2 <- lm(con ~ gnp + con1 + gnp1, data = usdl) ## various equivalent specifications of the LR test ## 下面4种操作方法是等价的 lrtest(fm2, fm1) Likelihood ratio test Model 1: con ~ gnp + con1 + gnp1 Model 2: con ~ gnp + gnp1 #Df LogLik Df Chisq Pr(>Chisq) 1 5 -56.069 2 4 -65.871 -1 19.605 9.524e-06 *** --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 lrtest(fm2, 2) # Remove variable with idx 2 Likelihood ratio test Model 1: con ~ gnp + con1 + gnp1 Model 2: con ~ gnp + gnp1 #Df LogLik Df Chisq Pr(>Chisq) 1 5 -56.069 2 4 -65.871 -1 19.605 9.524e-06 *** --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 lrtest(fm2, "con1") # Remove variable cond 1 Likelihood ratio test Model 1: con ~ gnp + con1 + gnp1 Model 2: con ~ gnp + gnp1 #Df LogLik Df Chisq Pr(>Chisq) 1 5 -56.069 2 4 -65.871 -1 19.605 9.524e-06 *** --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 lrtest(fm2, . ~ . - con1) # Remove con1 by formula Likelihood ratio test Model 1: con ~ gnp + con1 + gnp1 Model 2: con ~ gnp + gnp1 #Df LogLik Df Chisq Pr(>Chisq) 1 5 -56.069 2 4 -65.871 -1 19.605 9.524e-06 *** --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Survival 模型交互项： library(survival) fit2 = coxph(Surv(time, status) ~ age + sex*factor(ph.ecog), data = lung) lmtest::lrtest(fit2, . ~ . - sex:factor(ph.ecog)) Likelihood ratio test Model 1: Surv(time, status) ~ age + sex * factor(ph.ecog) Model 2: Surv(time, status) ~ age + sex + factor(ph.ecog) #Df LogLik Df Chisq Pr(>Chisq) 1 7 -728.93 2 5 -729.05 -2 0.2263 0.893 summary(fit2) Call: coxph(formula = Surv(time, status) ~ age + sex * factor(ph.ecog), data = lung) n= 227, number of events= 164 (1 observation deleted due to missingness) coef exp(coef) se(coef) z Pr(>|z|) age 0.009667 1.009714 0.009601 1.007 0.3140 sex -0.641345 0.526583 0.385890 -1.662 0.0965 . factor(ph.ecog)1 0.335017 1.397964 0.609291 0.550 0.5824 factor(ph.ecog)2 0.618092 1.855384 0.686960 0.900 0.3683 factor(ph.ecog)3 1.941570 6.969687 1.033409 1.879 0.0603 . sex:factor(ph.ecog)1 0.063350 1.065400 0.451918 0.140 0.8885 sex:factor(ph.ecog)2 0.221865 1.248403 0.507016 0.438 0.6617 sex:factor(ph.ecog)3 NA NA 0.000000 NA NA --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 exp(coef) exp(-coef) lower .95 upper .95 age 1.0097 0.9904 0.9909 1.029 sex 0.5266 1.8990 0.2472 1.122 factor(ph.ecog)1 1.3980 0.7153 0.4235 4.615 factor(ph.ecog)2 1.8554 0.5390 0.4827 7.131 factor(ph.ecog)3 6.9697 0.1435 0.9195 52.827 sex:factor(ph.ecog)1 1.0654 0.9386 0.4394 2.583 sex:factor(ph.ecog)2 1.2484 0.8010 0.4621 3.372 sex:factor(ph.ecog)3 NA NA NA NA Concordance= 0.636 (se = 0.025 ) Likelihood ratio test= 31.09 on 7 df, p=6e-05 Wald test = 30.72 on 7 df, p=7e-05 Score (logrank) test = 33.92 on 7 df, p=2e-05 除了用 likelihood ratio test，还可以使用 anova() 去对比。 但如果数据不一致时，上述模型无法见效。这时，应把不同的对比作为变量构建模型，然后可以采用 aov() 进行方差分析，看变量是否显著以评估是否有组间差异。