使用R語言做邏輯回歸的時候,當自變量中有分類變量(大於兩個)的時候,對於回歸模型的結果有一點困惑,搜索相關知識發現不少人也有相同的疑問,通過查閱資料這里給出自己的理解。
首先看一個實例(數據下載自:http://freakonometrics.free.fr/db.txt)
> db <- read.table("db.txt",header=TRUE,sep=";")
> head(db)
Y X1 X2 X3 1 1 3.297569 16.25411 B 2 1 6.418031 18.45130 D 3 1 5.279068 16.61806 B 4 1 5.539834 19.72158 C 5 1 4.123464 18.38634 C 6 1 7.778443 19.58338 C
> summary(db)
Y X1 X2 X3 Min. :0.000 Min. :-1.229 Min. :10.93 A:197 1st Qu.:1.000 1st Qu.: 4.545 1st Qu.:17.98 B:206 Median :1.000 Median : 5.982 Median :20.00 C:196 Mean :0.921 Mean : 5.958 Mean :19.94 D:197 3rd Qu.:1.000 3rd Qu.: 7.358 3rd Qu.:21.89 E:204 Max. :1.000 Max. :11.966 Max. :28.71
> reg <- glm(Y~X1+X2+X3,family=binomial,data=db)
> summary(reg)
Call: glm(formula = Y ~ X1 + X2 + X3, family = binomial, data = db) Deviance Residuals: Min 1Q Median 3Q Max -2.98017 0.09327 0.19106 0.37000 1.50646 Coefficients: Estimate Std. Error z value Pr(>|z|) (Intercept) -4.45885 1.04646 -4.261 2.04e-05 *** X1 0.51664 0.11178 4.622 3.80e-06 *** X2 0.21008 0.07247 2.899 0.003745 ** X3B 1.74496 0.49952 3.493 0.000477 *** X3C -0.03470 0.35691 -0.097 0.922543 X3D 0.08004 0.34916 0.229 0.818672 X3E 2.21966 0.56475 3.930 8.48e-05 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 (Dispersion parameter for binomial family taken to be 1) Null deviance: 552.64 on 999 degrees of freedom Residual deviance: 397.69 on 993 degrees of freedom AIC: 411.69 Number of Fisher Scoring iterations: 7
該數據集三個自變量中 X1, X2為連續型變量,X3為分類變量(A,B,C,D,E)。 獲取邏輯回歸結果時發現X3變量的表示形式和X1,X2不一樣,並且分別產生了X3B, X3C, X3D, X3E四個新的變量,但是又沒有X3A變量。后來查閱相關資料才明白原來邏輯回歸中處理分類變量和連續型變量是不一樣的。
當分類自變量的類別大於兩個的時候,需要建立一組虛擬變量(啞變量)來代表變量的歸屬性質。一般虛擬變量的數目比分類變量的數目少一個,少掉的那個就作為參照類(reference category)。例如本例中,A就是參照類,X3B, X3C, X3D, X3E就是四個虛擬變量。參照類的選取是隨意的,R語言邏輯回歸默認將分類變量的第一個factor設置為虛擬變量。此時的回歸模型如下:
四個虛擬變量的取值為1或0,即當觀測值中的分類變量屬於某一組時,該組的虛擬變量值為1,剩下的虛擬變量值為0。
例如,當一組觀測值(X1,X2,X3,Y)中 X3 的值為B時,虛擬變量X3B = 1, X3C, X3D, X3E 都為0,此時:
而當一組觀測值(X1,X2,X3,Y)中 X3 的值為A時, 因為A為參照類,所以此時X3B, X3C, X3D, X3E都為0,此時:
因此在控制變量條件下,即假設兩組觀測值中,X1, X2相同,而X3分別為A和B, 由上面兩式相減可得:
此處odds(B/A)為變量B對變量A的發生比率,即變量B的發生比與變量A的發生比的比值。大於1的發生比率表明事件發生的可能性會提高,或者說自變量對事件發生的概率有正的作用。例如,假如說odds(B/A)的數值大於1,那么說明在X1,X2不變的條件下,X3取值B比X3取值A有更大的概率使Y的值為1。(王濟川,郭志剛. Logistic 回歸模型 —— 方法與應用[M]. 北京:高等教育出版社)
回到開頭的例子,根據結果我們得以得出這樣的結論,變量X3取值A,C,D對Y的影響差不多,而變量X3取值B,E會使得Y取值為1的概率比去A,C,D顯著增大。簡單看一下:
> db_a <- db[db$X3 == "A",] > db_b <- db[db$X3 == "B",] > db_c <- db[db$X3 == "C",] > db_d <- db[db$X3 == "D",] > db_e <- db[db$X3 == "E",] > table(db_a$Y) 0 1 25 172 > table(db_b$Y) 0 1 6 200 > table(db_c$Y) 0 1 21 175 > table(db_d$Y) 0 1 22 175 > table(db_e$Y) 0 1 5 199
大致從結果看出確實變量B,E組的Y值為1的比例要高於A,C,D組。
我們也可以自己定義虛擬變量:
> levels(db$X3) [1] "A" "B" "C" "D" "E" > db$X3 <- relevel(db$X3, "B") > levels(db$X3) [1] "B" "A" "C" "D" "E"
同上面的回歸模型:
> reg <- glm(Y~X1+X2+X3,family=binomial,data=db) > summary(reg) Call: glm(formula = Y ~ X1 + X2 + X3, family = binomial, data = db) Deviance Residuals: Min 1Q Median 3Q Max -2.98017 0.09327 0.19106 0.37000 1.50646 Coefficients: Estimate Std. Error z value Pr(>|z|) (Intercept) -2.71389 1.07274 -2.530 0.011410 * X1 0.51664 0.11178 4.622 3.8e-06 *** X2 0.21008 0.07247 2.899 0.003745 ** X3A -1.74496 0.49952 -3.493 0.000477 *** X3C -1.77966 0.51002 -3.489 0.000484 *** X3D -1.66492 0.50365 -3.306 0.000947 *** X3E 0.47470 0.66354 0.715 0.474364 --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 (Dispersion parameter for binomial family taken to be 1) Null deviance: 552.64 on 999 degrees of freedom Residual deviance: 397.69 on 993 degrees of freedom AIC: 411.69 Number of Fisher Scoring iterations: 7
主要內容就這么多,如果想要更詳細的了解可以參考:王濟川,郭志剛. Logistic 回歸模型 —— 方法與應用[M]. 北京:高等教育出版社
以及鏈接:https://www.r-bloggers.com/logistic-regression-and-categorical-covariates/
版權聲明:本文為博主原創文章,博客地址:http://www.cnblogs.com/Demo1589/p/8973731.html,轉載請注明出處。