有谁知道如何使用 mclogit 中的 mblogit 函数从多项逻辑回归模型创建分类表。
我在网上找不到任何例子,感到非常沮丧。我可以毫无问题地从 multinom 函数获取分类表
这是我的模型的代码
library(mclogit)
multinom_model <- mblogit(Spawn_ID2 ~ average_Ca + average_Mn + average_Zn, data = train, random = ~1 | ID)
数据
train <- structure(list(Spawn_ID2 = structure(c(1L, 1L, 2L, 3L, 3L, 1L,
1L, 1L, 1L, 2L, 3L, 3L, 3L, 1L, 1L, 2L, 2L, 3L, 1L, 1L, 1L, 2L,
2L, 3L, 3L, 1L, 3L, 2L, 3L, 1L, 1L, 2L, 2L, 2L, 3L, 3L, 1L, 2L,
3L, 3L, 3L, 1L, 1L, 2L, 2L, 2L, 3L, 3L, 3L, 2L, 3L, 2L, 3L, 3L,
1L, 1L, 3L, 3L, 1L, 2L, 3L, 3L, 1L, 1L, 3L, 1L, 1L, 1L, 2L, 3L,
3L, 1L, 2L, 3L, 3L, 1L, 1L, 1L, 2L, 2L, 3L, 3L, 3L), .Label = c("1",
"2", "3"), class = "factor"), average_Ca = c(2365097.923, 2435299,
2447003.947, 2345150.667, 2442557.556, 3273328.407, 3270609.56,
3351219.8, 3114987.875, 3165870.444, 3291476.448, 3359619.591,
3016741.214, 1106471.379, 1118023.857, 1137698.28, 1148620.348,
1109930.897, 1760035.014, 1773782.642, 1787187.905, 1691737.88,
1807261.653, 1657850.143, 1798622.771, 1525681.789, 1586493.148,
1429283.944, 1606884.929, 1946063.84, 1841356.667, 2050657.087,
2002965.429, 1911144.875, 1921007.467, 1901236, 2247586.417,
2206708.758, 1840172.944, 1999647.143, 2274722.219, 1167092.839,
1428112.038, 1332175.833, 1290003.867, 1105272.952, 1306802.75,
1173046.273, 1436822.324, 1227869.143, 1183286.154, 1273745.257,
1282425.769, 1339544.737, 4149021.725, 4459955.917, 3787393.672,
4454623.645, 2786507.158, 2755766.667, 2835629.517, 2715092.1,
3039458.864, 3223676.619, 3246355.824, 2585769, 2460529.316,
2548687.2, 2560059.95, 2637233.556, 2568379.765, 1633590.579,
1560962.273, 1654643.316, 1578421, 1414345.261, 1516850.591,
1560330.4, 1461838.893, 1565908.448, 1412577.647, 1448740, 1528982.063
), average_Mn = c(48.80424538, 51.06187909, 57.30724947, 48.09621286,
52.51352389, 36.08437333, 39.8565816, 87.17658667, 99.6882075,
29.64281944, 29.3972869, 80.81569273, 84.66107357, 20.19708207,
17.41069214, 20.387258, 22.29570696, 18.87749483, 43.43359357,
34.87588094, 54.75980762, 34.6280406, 42.9865002, 37.78165262,
42.765985, 34.34720632, 15.02555489, 44.72356, 28.360175, 33.924748,
58.23363333, 39.03274957, 58.37289476, 53.66621875, 60.42372933,
57.034178, 142.1726167, 122.792637, 62.11736722, 75.77778107,
146.6532081, 18.52475484, 13.38198708, 13.06642123, 22.77531778,
15.03726, 12.33745129, 22.65816045, 13.89419776, 36.11867286,
30.50315769, 49.28838743, 36.76186077, 64.12796211, 94.43138575,
95.7703775, 24.0587618, 94.9781371, 90.87561211, 64.18755143,
43.22726966, 70.2781455, 23.74958045, 26.45033048, 26.65918529,
27.59556313, 48.42888105, 78.189086, 26.519214, 26.35578444,
89.83558706, 39.10770263, 38.60200364, 38.73780158, 35.15071813,
23.48386348, 58.94505136, 82.256964, 29.4669725, 65.96776966,
21.99724235, 42.78055176, 76.0168925), average_Zn = c(245.3482731,
305.0549818, 225.5911947, 214.8873857, 306.4080778, 170.2043667,
194.158624, 289.9295867, 246.5235125, 165.7625722, 152.5834483,
289.1114, 230.3927929, 97.88223655, 81.73639143, 107.9251376,
146.4331, 95.72906069, 122.5505496, 114.468913, 88.85797238,
107.9937798, 99.42765429, 111.497753, 93.50748042, 153.0925632,
103.3344659, 187.9697611, 172.0182571, 131.744836, 346.7786333,
192.9989435, 321.3208238, 336.208675, 308.7409067, 335.676, 163.0427083,
221.0996212, 100.6110142, 132.68165, 202.0647875, 151.1675742,
100.6805067, 98.782744, 130.9494524, 114.2988481, 93.53012708,
178.1473955, 109.4961284, 183.6133381, 190.1452077, 180.7220914,
173.6134154, 261.4983868, 178.378095, 212.3042, 141.4086787,
203.6961742, 298.6220947, 295.2357476, 272.2689172, 275.089395,
302.9338182, 263.3047762, 340.78345, 220.1297063, 111.2658379,
133.30894, 199.98872, 209.7802, 159.7065306, 205.9236526, 205.7157364,
214.2047526, 285.9215938, 143.694913, 254.8733864, 334.40914,
199.9216964, 295.7109379, 133.1236882, 200.1539529, 336.4642813
), ID = c("BLA14", "BLA14", "BLA14", "BLA14", "BLA14", "BLA2",
"BLA2", "BLA2", "BLA2", "BLA2", "BLA2", "BLA2", "BLA2", "BLA20",
"BLA20", "BLA20", "BLA20", "BLA20", "BLA209", "BLA209", "BLA209",
"BLA209", "BLA209", "BLA209", "BLA209", "BLA21", "BLA21", "BLA211",
"BLA211", "BLA238", "BLA238", "BLA238", "BLA238", "BLA238", "BLA238",
"BLA238", "BLA24", "BLA24", "BLA24", "BLA24", "BLA24", "BLA25",
"BLA25", "BLA25", "BLA25", "BLA25", "BLA25", "BLA25", "BLA25",
"BLA283", "BLA283", "BLA307", "BLA307", "BLA307", "BLA31", "BLA31",
"BLA31", "BLA31", "BLA42", "BLA42", "BLA42", "BLA42", "BLA47",
"BLA47", "BLA47", "BLA5", "BLA5", "BLA5", "BLA5", "BLA5", "BLA5",
"BLA79", "BLA79", "BLA79", "BLA79", "BLA80", "BLA80", "BLA80",
"BLA80", "BLA80", "BLA80", "BLA80", "BLA80")), row.names = c(NA,
-83L), class = c("tbl_df", "tbl", "data.frame"))
您可能想要一个混淆矩阵,您可以通过
table
计算实际与预测来获得它。在下面的 predict
离子中,模型会选择概率最高的离子。
preds <- predict(multinom_model, type='response')
head(preds)
# 1 2 3
# [1,] 0.3552168 0.2104946 0.4342886
# [2,] 0.3273301 0.2206145 0.4520554
# [3,] 0.3866728 0.1867448 0.4265824
# [4,] 0.3722244 0.2027590 0.4250166
# [5,] 0.3286972 0.2190235 0.4522793
# [6,] 0.4981363 0.1022896 0.3995741
所以你想要
apply
which.max
到每一行。
table(actual=train$Spawn_ID2, pred=apply(preds, MARGIN=1, which.max))
# pred
# actual 1 2 3
# 1 10 3 16
# 2 2 4 15
# 3 8 3 22
您可以进一步将结果通过管道传输到
proportions
。
table(actual=train$Spawn_ID2, pred=apply(preds, 1, which.max)) |> proportions(margin=1)
# pred
# actual 1 2 3
# 1 0.34482759 0.10344828 0.55172414
# 2 0.09523810 0.19047619 0.71428571
# 3 0.24242424 0.09090909 0.66666667
多项式混合模型方程
library(mclogit)
multinom_model <- mblogit(Spawn_ID2 ~ average_Ca + average_Mn + average_Zn, data = train, random = ~1 | ID)
获取混合多项式模型的参数
coefs <- summary(multinom_model)$coefficients
LLs <- coefs[,1] + qnorm(.025)*coefs[,2]
ULs <- coefs[,1] + qnorm(.975)*coefs[,2]
OR <- exp(coefs[,1])
ORLL <- exp(LLs)
ORUL <- exp(ULs)
HHES <- coefs[,1]/1.81 # Hasselblad and Hedges Effect Size
获取混合多项式模型的系数
round(cbind(coefs, LLs, ULs), 3)
获取混合多项模型的优势比
round(cbind(OR, ORLL, ORUL, HHES), 3)