Ggpredict 过度拟合导致趋势线从奇怪的地方开始
Ggpredict over fitting causing trend line to begin in an odd place
我不完全确定我在做什么 wrong/what 来查找。我的objective是用ggpredict和ggplot来展示时间和燃烧年数比例的关系。我猜这与对数转换的时间变量有关?
library(lme4);library(ggplot2);library(ggeffects);library(dplyr)
data = read.csv('FtmpAllyrs10kmC.csv')
数据是这样的:
structure(list(Observ = c(5208, 2828, 1664, 578, 18, 1644, 4741,
751, 689, 3813, 1464, 438, 1553, 4752, 4960, 376, 2482, 1811,
5682, 5441, 4505, 2281, 2103, 2993, 562, 4297, 3592, 5148, 3793,
1621, 1912, 1627, 1737, 4976, 2173, 5132, 5758, 2756, 1789, 5666,
2628, 2593, 794, 5779, 5158, 3123, 4986, 676, 4200, 2442, 2751,
4330, 1802, 2020, 2500, 1056, 959, 3290, 4303, 247, 5586, 922,
1049, 2432, 2076, 2560, 1369, 3636, 3722, 4137, 1561, 4915, 2515,
3034, 5547, 1491, 1247, 4116, 455, 4687, 1697, 5329, 21, 5724,
3701, 5697, 2938, 1721, 61, 998, 4304, 5798, 651, 910, 2689,
3986, 2908, 5753, 2574, 2345, 1940, 4317, 4588, 2179, 665, 4133,
749, 3977, 3134, 4190, 3985, 4937, 2473, 3238, 4987, 3915, 4261,
3521, 2736, 3665, 1797, 5692, 5578, 4087, 2011, 903, 889, 1523,
3396, 2291, 5269, 3644, 3403, 4814, 4618, 16, 77, 5385, 2842,
5816, 2015, 1443, 3183, 3331, 4977, 5380, 989, 4918, 740, 4637,
887, 1557, 4295, 4673, 1918, 5662, 4167, 1384, 3441, 614, 2360,
780, 661, 1267, 2018, 1906, 3402, 677, 5218, 2830, 4979, 3984,
4924, 1125, 2640, 986, 1885, 2573, 5300, 2398, 4832, 4816, 3738,
3276, 3830, 2425, 2054, 4273, 5607, 1678, 378, 1158, 510, 2210,
2399, 1952, 2909, 4945, 2659, 2642), yrblock15 = c(2015, 2010,
2007, 2005, 2004, 2007, 2014, 2005, 2005, 2012, 2007, 2004, 2007,
2014, 2015, 2004, 2009, 2008, 2016, 2016, 2014, 2009, 2008, 2010,
2005, 2013, 2011, 2015, 2012, 2007, 2008, 2007, 2007, 2015, 2008,
2015, 2016, 2010, 2007, 2016, 2009, 2009, 2005, 2016, 2015, 2010,
2015, 2005, 2013, 2009, 2010, 2013, 2008, 2008, 2009, 2006, 2006,
2011, 2013, 2004, 2016, 2006, 2006, 2009, 2008, 2009, 2007, 2012,
2012, 2013, 2007, 2014, 2009, 2010, 2016, 2007, 2006, 2013, 2005,
2014, 2007, 2015, 2004, 2016, 2012, 2016, 2010, 2007, 2004, 2006,
2013, 2016, 2005, 2006, 2009, 2012, 2010, 2016, 2009, 2009, 2008,
2013, 2014, 2008, 2005, 2013, 2005, 2012, 2010, 2013, 2012, 2014,
2009, 2011, 2015, 2012, 2013, 2011, 2010, 2012, 2007, 2016, 2016,
2013, 2008, 2006, 2005, 2007, 2011, 2009, 2015, 2012, 2011, 2014,
2014, 2004, 2004, 2015, 2010, 2016, 2008, 2007, 2011, 2011, 2015,
2015, 2006, 2014, 2005, 2014, 2005, 2007, 2013, 2014, 2008, 2016,
2013, 2007, 2011, 2005, 2009, 2005, 2005, 2006, 2008, 2008, 2011,
2005, 2015, 2010, 2015, 2012, 2014, 2006, 2009, 2006, 2008, 2009,
2015, 2009, 2014, 2014, 2012, 2011, 2012, 2009, 2008, 2013, 2016,
2007, 2004, 2006, 2005, 2008, 2009, 2008, 2010, 2014, 2009, 2009
), circleID = c(258, 128, 314, 128, 18, 294, 241, 301, 239, 213,
114, 438, 203, 252, 10, 376, 232, 11, 282, 41, 5, 31, 303, 293,
112, 247, 442, 198, 193, 271, 112, 277, 387, 26, 373, 182, 358,
56, 439, 266, 378, 343, 344, 379, 208, 423, 36, 226, 150, 192,
51, 280, 2, 220, 250, 156, 59, 140, 253, 247, 186, 22, 149, 182,
276, 310, 19, 36, 122, 87, 211, 415, 265, 334, 147, 141, 347,
66, 5, 187, 347, 379, 21, 324, 101, 297, 238, 371, 61, 98, 254,
398, 201, 10, 439, 386, 208, 353, 324, 95, 140, 267, 88, 379,
215, 83, 299, 377, 434, 140, 385, 437, 223, 88, 37, 315, 211,
371, 36, 65, 447, 292, 178, 37, 211, 3, 439, 173, 246, 41, 319,
44, 253, 314, 118, 16, 77, 435, 142, 416, 215, 93, 33, 181, 27,
430, 89, 418, 290, 137, 437, 207, 245, 173, 118, 262, 117, 34,
291, 164, 110, 330, 211, 367, 218, 106, 252, 227, 268, 130, 29,
384, 424, 225, 390, 86, 85, 323, 350, 148, 332, 316, 138, 126,
230, 175, 254, 223, 207, 328, 378, 258, 60, 410, 149, 152, 209,
445, 409, 392), rain15 = c(347.83, 394.12, 382.2, 382.41, 395.7,
386.08, 383.79, 352.65, 354.31, 366.48, 416.79, 335.17, 409.24,
373, 390.76, 341.35, 387.25, 452.18, 329.14, 365.74, 432.58,
443.36, 375.57, 359.75, 379.14, 386.41, 361.47, 366.1, 382.57,
383.32, 409.56, 390.92, 380.38, 394.94, 366.72, 347.44, 336.88,
410.94, 370.83, 335.88, 368.53, 370.42, 344.56, 323.41, 348.34,
351.07, 382.75, 362.64, 402.7, 396.11, 418.01, 389.14, 462.76,
391.05, 369.47, 399.78, 419.32, 392.97, 389.15, 345.37, 336.22,
405.73, 378.45, 394.7, 388.29, 379.56, 437.29, 415.95, 388.91,
402.43, 397.09, 368.84, 378.54, 361.92, 355.22, 416.46, 361.24,
417.12, 420.92, 386.48, 375.04, 335.03, 385.23, 342.51, 401.27,
341.21, 362.81, 372.85, 396.48, 390.72, 385.06, 343.64, 365.25,
440.76, 364.68, 354.45, 368.7, 324.44, 366.4, 408.43, 405.71,
390.8, 401.09, 364.07, 360.68, 399.39, 348.38, 344.2, 345.23,
401.29, 356.48, 364.21, 376.12, 403.37, 384.1, 355.71, 389.53,
363.28, 417.76, 403.16, 362.28, 333.91, 337.46, 419.51, 389.22,
448.08, 338.46, 397.52, 372.25, 424.25, 349.25, 408.19, 376.68,
375.87, 403.78, 398.73, 386.92, 340.39, 391.58, 335.03, 390.25,
422.05, 423.79, 386.49, 392.97, 334.07, 403.85, 369.54, 348.84,
392.33, 336.68, 399.56, 386.84, 395.97, 409.93, 337.08, 410.27,
450.48, 364.93, 369.08, 413.31, 341.93, 360.06, 362.28, 395.8,
423.56, 376.67, 366.19, 358.88, 390.74, 390.84, 362.84, 370.21,
360.84, 371.9, 410.36, 421.59, 367.48, 355.62, 389.61, 370.81,
374.37, 382.61, 401.78, 373.7, 382.72, 387.56, 388.53, 329.06,
383.78, 336.97, 376.68, 398.57, 370.46, 388.88, 421.66, 369.29,
371.58, 369.01, 369.22), YearsBurnt = c(6, 6, 3.5, 5, 3, 2, 3.5,
2.5, 2, 1.5, 10.5, 3.5, 2.5, 3.5, 4.5, 3, 2, 2.5, 1.5, 3.5, 3.5,
4, 4, 3, 3.5, 2.5, 6, 4.5, 4, 2.5, 3.5, 2, 7, 3, 2.5, 3.5, 13,
3, 3.5, 3.5, 4.5, 3, 1.5, 2, 4, 2, 4.5, 4, 3.5, 2.5, 2, 2, 3,
1, 5, 2.5, 4, 12.5, 2.5, 1.5, 3.5, 1.5, 2.5, 4, 4.5, 10, 3, 3.5,
4.5, 10.5, 1, 4.5, 2, 13.5, 8.5, 10, 1, 4, 3, 3.5, 1.5, 3, 2.5,
2.5, 2.5, 4.5, 4, 1.5, 3, 3.5, 4.5, 1.5, 3, 2.5, 3.5, 8.5, 4,
7, 2.5, 5, 11, 3.5, 11.5, 3, 1.5, 3, 0.5, 4.5, 3.5, 13.5, 7.5,
3.5, 2, 12, 4, 5, 2, 1.5, 3.5, 4.5, 2, 3.5, 3, 4, 1.5, 2, 2.5,
6, 2, 5, 3.5, 4.5, 2, 3.5, 5, 4.5, 3, 4, 14, 3, 1.5, 3.5, 5.5,
3, 4, 3, 7, 4.5, 2.5, 3, 3, 3.5, 3, 9, 5, 6.5, 5, 4, 4, 3.5,
3, 8.5, 1, 4.5, 1.5, 5.5, 3, 2, 2.5, 2.5, 3, 8.5, 2.5, 1, 3.5,
5.5, 5, 1.5, 2, 4.5, 5, 4, 1.5, 3.5, 4.5, 6, 4.5, 3.5, 3, 6.5,
3, 6.5, 3.5, 4.5, 2.5, 2.5, 4, 4, 4, 4.5), YearsNotBurnt = c(9,
9, 11.5, 10, 12, 13, 11.5, 12.5, 13, 13.5, 4.5, 11.5, 12.5, 11.5,
10.5, 12, 13, 12.5, 13.5, 11.5, 11.5, 11, 11, 12, 11.5, 12.5,
9, 10.5, 11, 12.5, 11.5, 13, 8, 12, 12.5, 11.5, 2, 12, 11.5,
11.5, 10.5, 12, 13.5, 13, 11, 13, 10.5, 11, 11.5, 12.5, 13, 13,
12, 14, 10, 12.5, 11, 2.5, 12.5, 13.5, 11.5, 13.5, 12.5, 11,
10.5, 5, 12, 11.5, 10.5, 4.5, 14, 10.5, 13, 1.5, 6.5, 5, 14,
11, 12, 11.5, 13.5, 12, 12.5, 12.5, 12.5, 10.5, 11, 13.5, 12,
11.5, 10.5, 13.5, 12, 12.5, 11.5, 6.5, 11, 8, 12.5, 10, 4, 11.5,
3.5, 12, 13.5, 12, 14.5, 10.5, 11.5, 1.5, 7.5, 11.5, 13, 3, 11,
10, 13, 13.5, 11.5, 10.5, 13, 11.5, 12, 11, 13.5, 13, 12.5, 9,
13, 10, 11.5, 10.5, 13, 11.5, 10, 10.5, 12, 11, 1, 12, 13.5,
11.5, 9.5, 12, 11, 12, 8, 10.5, 12.5, 12, 12, 11.5, 12, 6, 10,
8.5, 10, 11, 11, 11.5, 12, 6.5, 14, 10.5, 13.5, 9.5, 12, 13,
12.5, 12.5, 12, 6.5, 12.5, 14, 11.5, 9.5, 10, 13.5, 13, 10.5,
10, 11, 13.5, 11.5, 10.5, 9, 10.5, 11.5, 12, 8.5, 12, 8.5, 11.5,
10.5, 12.5, 12.5, 11, 11, 11, 10.5), time = c(1.96, 4.94, 3.46,
4.94, 2.73, 6.22, 4.5, 2.67, 4.66, 3.83, 0.38, 2.6, 3.97, 4.18,
3.77, 3.44, 2.9, 3.93, 2.16, 3.51, 2.91, 3.19, 2.73, 6.36, 1.74,
4.39, 4.1, 2.26, 2.36, 5.32, 1.74, 3.66, 1.26, 5.61, 9.04, 4.61,
0.46, 3.98, 2.63, 5.5, 2.56, 5.92, 6.39, 2.26, 3.27, 7.95, 2.93,
4.93, 2.97, 2.43, 5.91, 3.07, 4.27, 3.21, 4.12, 4.72, 1.93, 0.69,
3.51, 4.39, 4.02, 3.18, 2.61, 4.61, 3.67, 0.54, 2.33, 2.93, 2.12,
1.06, 3.95, 2.31, 5.44, 0.17, 1.42, 0.55, 8.35, 2.53, 2.91, 3.26,
8.35, 2.26, 2.23, 7.18, 6.59, 6.36, 4.38, 7.67, 1.93, 3.34, 2.91,
8.54, 5.75, 3.77, 2.63, 0.97, 3.27, 1.58, 7.18, 2.08, 0.69, 5.43,
0.85, 2.26, 3.69, 3.18, 6.18, 2.93, 2.68, 0.69, 0.92, 2.34, 3.26,
0.85, 2.91, 4.3, 3.95, 7.67, 2.93, 2.1, 6.54, 6.31, 3.87, 2.91,
3.95, 3.35, 2.63, 1.49, 4.32, 3.51, 7.06, 2.67, 3.51, 3.46, 1.56,
4.33, 5.64, 2.73, 0.57, 2.87, 3.69, 2.56, 2.33, 4.27, 4.73, 4.02,
0.82, 4.11, 4.88, 2.29, 2.34, 3.72, 4.21, 1.49, 1.56, 3.03, 1.24,
2.65, 5.71, 1.67, 2.71, 1.49, 3.95, 4.51, 3.36, 5.21, 4.18, 4.54,
5.36, 4.25, 3.71, 0.95, 8.92, 3.12, 2.73, 1.36, 1.85, 7.24, 8.11,
2.2, 0.95, 5.16, 1.3, 6.54, 3.01, 1.97, 2.91, 3.26, 3.72, 1.79,
2.56, 1.96, 1.89, 1.89, 2.61, 5.25, 3.25, 5.26, 1.74, 3.73),
claylake = c(0, 0, 0, 0, 0, 17.53, 0.1, 0.59, 0, 9.13, 36.93,
12.75, 0, 0, 0, 0, 0, 0, 0, 0.09, 0.01, 0, 0, 9.43, 74.71,
26.42, 0.23, 0, 0, 35.27, 74.71, 0, 0, 0, 0, 0, 0, 0, 20.81,
9.46, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
1.14, 0, 26.42, 3.62, 0, 0, 0, 0.21, 0, 0, 0, 0.03, 10.43,
0.99, 3.6, 5.32, 0, 0.36, 0, 0, 0.25, 0.01, 0.22, 0, 0, 6.45,
0, 0, 0, 0, 0, 1.71, 0, 0, 0, 0, 0, 20.81, 0, 0, 0.18, 0,
0, 1.14, 0.03, 1.2, 0, 8.97, 0, 0, 0, 0, 1.14, 0, 1.56, 0.22,
1.2, 0, 0, 0.99, 0, 0, 0, 0, 4.14, 0, 0, 0.99, 0, 20.81,
0, 33.61, 0.09, 14.94, 0, 0, 0, 0, 0.41, 0, 2.7, 0, 0.61,
8.97, 0, 0, 0, 0, 1.7, 2.67, 7.71, 0.2, 8.63, 1.56, 0, 0.49,
0, 0, 0, 0, 0, 11.9, 33.08, 0, 0, 0.99, 2.13, 0, 0, 0, 0,
0.03, 0, 0, 0, 0, 0, 0, 2.86, 1.65, 0, 0, 0, 0, 0, 60.14,
0, 0, 0, 0, 0.22, 0, 0, 0, 0, 0.3, 0, 0, 0, 0, 0, 0, 5.57
), spinsandplain = c(81.94, 34.29, 89.55, 34.29, 80.86, 75.92,
81.55, 43.53, 97.3, 87.84, 60.62, 80.81, 11.73, 5.11, 98.67,
79.52, 60.73, 91.65, 2.82, 97.31, 73.65, 72.78, 96.51, 74.02,
25.09, 50.74, 96.62, 88.77, 98.8, 54.04, 25.09, 95.1, 69.85,
99.4, 78.79, 78.77, 48.16, 80.68, 75.79, 66.33, 68.3, 79.11,
91.89, 82.49, 98.33, 90.82, 91.24, 65.01, 69.24, 99.94, 99.75,
18.57, 90.39, 95.56, 71.07, 67.85, 92.37, 85.85, 17.89, 50.74,
79.65, 68.82, 74.05, 78.77, 87.67, 41.11, 91.74, 91.24, 44.8,
86.24, 97.7, 94.17, 85.59, 33.53, 85.23, 94.55, 78.52, 95.49,
73.65, 95.04, 78.52, 82.49, 77.26, 83.4, 98.29, 85.24, 98.78,
87.09, 81.36, 96.62, 3.4, 94.65, 28.6, 98.67, 75.79, 73.34,
98.33, 74.88, 83.4, 88.24, 85.85, 52.44, 95.84, 82.49, 62.11,
98.74, 70.32, 86.18, 95.67, 85.85, 11.42, 85.96, 75.53, 95.84,
95.46, 93.68, 97.7, 87.09, 91.24, 80.03, 87.77, 68.71, 17.51,
95.46, 97.7, 50.7, 75.79, 70.43, 61.06, 97.31, 74.63, 99,
17.89, 89.55, 99.25, 98.08, 97.61, 93.36, 99.03, 38.1, 62.11,
96.9, 88.87, 40.48, 90.21, 73.79, 95.2, 66.53, 96.67, 82.89,
85.96, 97.08, 75.74, 70.43, 99.25, 96.4, 98.88, 98.13, 85.32,
54.19, 99.2, 81.42, 97.7, 82.25, 97.42, 98.1, 5.11, 12.06,
66.14, 52.39, 52.72, 12.32, 87.32, 98.95, 71.55, 90.58, 97.9,
80.62, 93.32, 76, 86.48, 86.42, 39.54, 68.65, 6.05, 86.02,
3.4, 75.53, 97.08, 32.47, 68.3, 81.94, 89.64, 57.4, 74.05,
0.47, 96.76, 86.7, 78.46, 84.81)), row.names = c(5208L, 2828L,
1664L, 578L, 18L, 1644L, 4741L, 751L, 689L, 3813L, 1464L, 438L,
1553L, 4752L, 4960L, 376L, 2482L, 1811L, 5682L, 5441L, 4505L,
2281L, 2103L, 2993L, 562L, 4297L, 3592L, 5148L, 3793L, 1621L,
1912L, 1627L, 1737L, 4976L, 2173L, 5132L, 5758L, 2756L, 1789L,
5666L, 2628L, 2593L, 794L, 5779L, 5158L, 3123L, 4986L, 676L,
4200L, 2442L, 2751L, 4330L, 1802L, 2020L, 2500L, 1056L, 959L,
3290L, 4303L, 247L, 5586L, 922L, 1049L, 2432L, 2076L, 2560L,
1369L, 3636L, 3722L, 4137L, 1561L, 4915L, 2515L, 3034L, 5547L,
1491L, 1247L, 4116L, 455L, 4687L, 1697L, 5329L, 21L, 5724L, 3701L,
5697L, 2938L, 1721L, 61L, 998L, 4304L, 5798L, 651L, 910L, 2689L,
3986L, 2908L, 5753L, 2574L, 2345L, 1940L, 4317L, 4588L, 2179L,
665L, 4133L, 749L, 3977L, 3134L, 4190L, 3985L, 4937L, 2473L,
3238L, 4987L, 3915L, 4261L, 3521L, 2736L, 3665L, 1797L, 5692L,
5578L, 4087L, 2011L, 903L, 889L, 1523L, 3396L, 2291L, 5269L,
3644L, 3403L, 4814L, 4618L, 16L, 77L, 5385L, 2842L, 5816L, 2015L,
1443L, 3183L, 3331L, 4977L, 5380L, 989L, 4918L, 740L, 4637L,
887L, 1557L, 4295L, 4673L, 1918L, 5662L, 4167L, 1384L, 3441L,
614L, 2360L, 780L, 661L, 1267L, 2018L, 1906L, 3402L, 677L, 5218L,
2830L, 4979L, 3984L, 4924L, 1125L, 2640L, 986L, 1885L, 2573L,
5300L, 2398L, 4832L, 4816L, 3738L, 3276L, 3830L, 2425L, 2054L,
4273L, 5607L, 1678L, 378L, 1158L, 510L, 2210L, 2399L, 1952L,
2909L, 4945L, 2659L, 2642L), class = "data.frame")
我创建了一个新变量,因为燃烧年数的比例超出了 15 年(即二项式)
data$fireprop = cbind(data$YearsBurnt,data$YearsNotBurnt)
型号:
mfireprop = glmer(fireprop~log(time)+spinsandplain+rain15+claylake+rain15*log(time)+(1|circleID),na.action=na.fail, family=binomial, data=data)
趋势线代码:
d = ggpredict(mfireprop, terms = "time[exp]")
d = rename(d, "time" = x, "fireprop" = predicted)
ggplot(d, aes(time, fireprop)) +
geom_ribbon(aes(ymin = conf.low, ymax = conf.high), alpha = .1) +
geom_line(size = 2, colour = "black") +
theme_bw()
趋势线如下所示:
为什么x轴在数据停止的10小时不停止?为什么会涨到20,000?为什么 y 轴只达到 0.4?当部分比例为1时?
当我限制 x 轴和 y 轴时,它最终看起来像这样:
但是当我查看上面的原始数据时,趋势线似乎从一个非常奇怪的地方开始。
我不确定我做错了什么?
好的,我已经找到这里的主要问题了。在 ggpredict()
函数的文档中有一个名为 back.transform
的参数,默认为 TRUE。这意味着经过对数转换的数据将自动转换回原始响应范围。这就是为什么如果您检查 ggpredict 对象 d
,您会发现 time
变量实际上确实超过了该对象中的 8000。所以因为你没有标记 back.transform=FALSE
,而且还指定了 time[exp]
,发生的事情是函数自动对你的值求幂,然后你又做了一次。
如果我们查看记录的值:
summary(log(data$time))
Min. 1st Qu. Median Mean 3rd Qu. Max.
-1.7720 0.8154 1.1802 1.0904 1.4793 2.2017
然后我们对最大值取幂,得到之前的最大值:
exp(2.2017) # Exponentiated to get back to years
[1] 9.040369
summary(df$time) # The original variable
Min. 1st Qu. Median Mean 3rd Qu. Max.
0.170 2.260 3.255 3.519 4.390 9.040
如果我们再次取幂,我们最终会得到超过 8000 的最大时间。
exp(9.040369)
[1] 8436.89
所以,为了得到你想要的情节,你只需要在 ggpredict()
:
调用时间后省略 [exp]
d = ggpredict(mfireprop, terms = "time")
d = rename(d, "time" = x, "fireprop" = predicted)
ggplot(d, aes(time, fireprop)) +
geom_ribbon(aes(ymin = conf.low, ymax = conf.high), alpha = .1) +
geom_line(size = 2, colour = "black") +
theme_bw()
时间被截断,因为在时间 0 没有变化。 YearsNotBurnt
始终为 0。因此,如果您从 ggpredict 查看对象 d
,您将在时间 0 的所有列中看到 NaN
。如果将模型简化为以下内容:
mfireprop2 = glmer(fireprop~
log(time) +
(1|circleID),
na.action=na.fail,
family=binomial,
data=df)
您将能够得到该图,但由于变化很小,置信区间将从 1 跨越到 0。我相信这是一个与分离有关的问题,基本上这意味着如果没有变化,或者如果某些东西可以完美地预测结果,二项式模型就不能适合常客模型。
我想我唯一想提的是你在评论中有一个关于“二项式 glm 中的非整数计数!”的问题。这是因为它期望因变量是试验的一部分,不应有小数。您的数据中有一些点似乎是半年的时间间隔。我对您的数据还不够熟悉,无法确定更好的选择是什么,但创建一个比例并在 weights=
参数中给出观察次数可能是一个选项。
我不完全确定我在做什么 wrong/what 来查找。我的objective是用ggpredict和ggplot来展示时间和燃烧年数比例的关系。我猜这与对数转换的时间变量有关?
library(lme4);library(ggplot2);library(ggeffects);library(dplyr)
data = read.csv('FtmpAllyrs10kmC.csv')
数据是这样的:
structure(list(Observ = c(5208, 2828, 1664, 578, 18, 1644, 4741,
751, 689, 3813, 1464, 438, 1553, 4752, 4960, 376, 2482, 1811,
5682, 5441, 4505, 2281, 2103, 2993, 562, 4297, 3592, 5148, 3793,
1621, 1912, 1627, 1737, 4976, 2173, 5132, 5758, 2756, 1789, 5666,
2628, 2593, 794, 5779, 5158, 3123, 4986, 676, 4200, 2442, 2751,
4330, 1802, 2020, 2500, 1056, 959, 3290, 4303, 247, 5586, 922,
1049, 2432, 2076, 2560, 1369, 3636, 3722, 4137, 1561, 4915, 2515,
3034, 5547, 1491, 1247, 4116, 455, 4687, 1697, 5329, 21, 5724,
3701, 5697, 2938, 1721, 61, 998, 4304, 5798, 651, 910, 2689,
3986, 2908, 5753, 2574, 2345, 1940, 4317, 4588, 2179, 665, 4133,
749, 3977, 3134, 4190, 3985, 4937, 2473, 3238, 4987, 3915, 4261,
3521, 2736, 3665, 1797, 5692, 5578, 4087, 2011, 903, 889, 1523,
3396, 2291, 5269, 3644, 3403, 4814, 4618, 16, 77, 5385, 2842,
5816, 2015, 1443, 3183, 3331, 4977, 5380, 989, 4918, 740, 4637,
887, 1557, 4295, 4673, 1918, 5662, 4167, 1384, 3441, 614, 2360,
780, 661, 1267, 2018, 1906, 3402, 677, 5218, 2830, 4979, 3984,
4924, 1125, 2640, 986, 1885, 2573, 5300, 2398, 4832, 4816, 3738,
3276, 3830, 2425, 2054, 4273, 5607, 1678, 378, 1158, 510, 2210,
2399, 1952, 2909, 4945, 2659, 2642), yrblock15 = c(2015, 2010,
2007, 2005, 2004, 2007, 2014, 2005, 2005, 2012, 2007, 2004, 2007,
2014, 2015, 2004, 2009, 2008, 2016, 2016, 2014, 2009, 2008, 2010,
2005, 2013, 2011, 2015, 2012, 2007, 2008, 2007, 2007, 2015, 2008,
2015, 2016, 2010, 2007, 2016, 2009, 2009, 2005, 2016, 2015, 2010,
2015, 2005, 2013, 2009, 2010, 2013, 2008, 2008, 2009, 2006, 2006,
2011, 2013, 2004, 2016, 2006, 2006, 2009, 2008, 2009, 2007, 2012,
2012, 2013, 2007, 2014, 2009, 2010, 2016, 2007, 2006, 2013, 2005,
2014, 2007, 2015, 2004, 2016, 2012, 2016, 2010, 2007, 2004, 2006,
2013, 2016, 2005, 2006, 2009, 2012, 2010, 2016, 2009, 2009, 2008,
2013, 2014, 2008, 2005, 2013, 2005, 2012, 2010, 2013, 2012, 2014,
2009, 2011, 2015, 2012, 2013, 2011, 2010, 2012, 2007, 2016, 2016,
2013, 2008, 2006, 2005, 2007, 2011, 2009, 2015, 2012, 2011, 2014,
2014, 2004, 2004, 2015, 2010, 2016, 2008, 2007, 2011, 2011, 2015,
2015, 2006, 2014, 2005, 2014, 2005, 2007, 2013, 2014, 2008, 2016,
2013, 2007, 2011, 2005, 2009, 2005, 2005, 2006, 2008, 2008, 2011,
2005, 2015, 2010, 2015, 2012, 2014, 2006, 2009, 2006, 2008, 2009,
2015, 2009, 2014, 2014, 2012, 2011, 2012, 2009, 2008, 2013, 2016,
2007, 2004, 2006, 2005, 2008, 2009, 2008, 2010, 2014, 2009, 2009
), circleID = c(258, 128, 314, 128, 18, 294, 241, 301, 239, 213,
114, 438, 203, 252, 10, 376, 232, 11, 282, 41, 5, 31, 303, 293,
112, 247, 442, 198, 193, 271, 112, 277, 387, 26, 373, 182, 358,
56, 439, 266, 378, 343, 344, 379, 208, 423, 36, 226, 150, 192,
51, 280, 2, 220, 250, 156, 59, 140, 253, 247, 186, 22, 149, 182,
276, 310, 19, 36, 122, 87, 211, 415, 265, 334, 147, 141, 347,
66, 5, 187, 347, 379, 21, 324, 101, 297, 238, 371, 61, 98, 254,
398, 201, 10, 439, 386, 208, 353, 324, 95, 140, 267, 88, 379,
215, 83, 299, 377, 434, 140, 385, 437, 223, 88, 37, 315, 211,
371, 36, 65, 447, 292, 178, 37, 211, 3, 439, 173, 246, 41, 319,
44, 253, 314, 118, 16, 77, 435, 142, 416, 215, 93, 33, 181, 27,
430, 89, 418, 290, 137, 437, 207, 245, 173, 118, 262, 117, 34,
291, 164, 110, 330, 211, 367, 218, 106, 252, 227, 268, 130, 29,
384, 424, 225, 390, 86, 85, 323, 350, 148, 332, 316, 138, 126,
230, 175, 254, 223, 207, 328, 378, 258, 60, 410, 149, 152, 209,
445, 409, 392), rain15 = c(347.83, 394.12, 382.2, 382.41, 395.7,
386.08, 383.79, 352.65, 354.31, 366.48, 416.79, 335.17, 409.24,
373, 390.76, 341.35, 387.25, 452.18, 329.14, 365.74, 432.58,
443.36, 375.57, 359.75, 379.14, 386.41, 361.47, 366.1, 382.57,
383.32, 409.56, 390.92, 380.38, 394.94, 366.72, 347.44, 336.88,
410.94, 370.83, 335.88, 368.53, 370.42, 344.56, 323.41, 348.34,
351.07, 382.75, 362.64, 402.7, 396.11, 418.01, 389.14, 462.76,
391.05, 369.47, 399.78, 419.32, 392.97, 389.15, 345.37, 336.22,
405.73, 378.45, 394.7, 388.29, 379.56, 437.29, 415.95, 388.91,
402.43, 397.09, 368.84, 378.54, 361.92, 355.22, 416.46, 361.24,
417.12, 420.92, 386.48, 375.04, 335.03, 385.23, 342.51, 401.27,
341.21, 362.81, 372.85, 396.48, 390.72, 385.06, 343.64, 365.25,
440.76, 364.68, 354.45, 368.7, 324.44, 366.4, 408.43, 405.71,
390.8, 401.09, 364.07, 360.68, 399.39, 348.38, 344.2, 345.23,
401.29, 356.48, 364.21, 376.12, 403.37, 384.1, 355.71, 389.53,
363.28, 417.76, 403.16, 362.28, 333.91, 337.46, 419.51, 389.22,
448.08, 338.46, 397.52, 372.25, 424.25, 349.25, 408.19, 376.68,
375.87, 403.78, 398.73, 386.92, 340.39, 391.58, 335.03, 390.25,
422.05, 423.79, 386.49, 392.97, 334.07, 403.85, 369.54, 348.84,
392.33, 336.68, 399.56, 386.84, 395.97, 409.93, 337.08, 410.27,
450.48, 364.93, 369.08, 413.31, 341.93, 360.06, 362.28, 395.8,
423.56, 376.67, 366.19, 358.88, 390.74, 390.84, 362.84, 370.21,
360.84, 371.9, 410.36, 421.59, 367.48, 355.62, 389.61, 370.81,
374.37, 382.61, 401.78, 373.7, 382.72, 387.56, 388.53, 329.06,
383.78, 336.97, 376.68, 398.57, 370.46, 388.88, 421.66, 369.29,
371.58, 369.01, 369.22), YearsBurnt = c(6, 6, 3.5, 5, 3, 2, 3.5,
2.5, 2, 1.5, 10.5, 3.5, 2.5, 3.5, 4.5, 3, 2, 2.5, 1.5, 3.5, 3.5,
4, 4, 3, 3.5, 2.5, 6, 4.5, 4, 2.5, 3.5, 2, 7, 3, 2.5, 3.5, 13,
3, 3.5, 3.5, 4.5, 3, 1.5, 2, 4, 2, 4.5, 4, 3.5, 2.5, 2, 2, 3,
1, 5, 2.5, 4, 12.5, 2.5, 1.5, 3.5, 1.5, 2.5, 4, 4.5, 10, 3, 3.5,
4.5, 10.5, 1, 4.5, 2, 13.5, 8.5, 10, 1, 4, 3, 3.5, 1.5, 3, 2.5,
2.5, 2.5, 4.5, 4, 1.5, 3, 3.5, 4.5, 1.5, 3, 2.5, 3.5, 8.5, 4,
7, 2.5, 5, 11, 3.5, 11.5, 3, 1.5, 3, 0.5, 4.5, 3.5, 13.5, 7.5,
3.5, 2, 12, 4, 5, 2, 1.5, 3.5, 4.5, 2, 3.5, 3, 4, 1.5, 2, 2.5,
6, 2, 5, 3.5, 4.5, 2, 3.5, 5, 4.5, 3, 4, 14, 3, 1.5, 3.5, 5.5,
3, 4, 3, 7, 4.5, 2.5, 3, 3, 3.5, 3, 9, 5, 6.5, 5, 4, 4, 3.5,
3, 8.5, 1, 4.5, 1.5, 5.5, 3, 2, 2.5, 2.5, 3, 8.5, 2.5, 1, 3.5,
5.5, 5, 1.5, 2, 4.5, 5, 4, 1.5, 3.5, 4.5, 6, 4.5, 3.5, 3, 6.5,
3, 6.5, 3.5, 4.5, 2.5, 2.5, 4, 4, 4, 4.5), YearsNotBurnt = c(9,
9, 11.5, 10, 12, 13, 11.5, 12.5, 13, 13.5, 4.5, 11.5, 12.5, 11.5,
10.5, 12, 13, 12.5, 13.5, 11.5, 11.5, 11, 11, 12, 11.5, 12.5,
9, 10.5, 11, 12.5, 11.5, 13, 8, 12, 12.5, 11.5, 2, 12, 11.5,
11.5, 10.5, 12, 13.5, 13, 11, 13, 10.5, 11, 11.5, 12.5, 13, 13,
12, 14, 10, 12.5, 11, 2.5, 12.5, 13.5, 11.5, 13.5, 12.5, 11,
10.5, 5, 12, 11.5, 10.5, 4.5, 14, 10.5, 13, 1.5, 6.5, 5, 14,
11, 12, 11.5, 13.5, 12, 12.5, 12.5, 12.5, 10.5, 11, 13.5, 12,
11.5, 10.5, 13.5, 12, 12.5, 11.5, 6.5, 11, 8, 12.5, 10, 4, 11.5,
3.5, 12, 13.5, 12, 14.5, 10.5, 11.5, 1.5, 7.5, 11.5, 13, 3, 11,
10, 13, 13.5, 11.5, 10.5, 13, 11.5, 12, 11, 13.5, 13, 12.5, 9,
13, 10, 11.5, 10.5, 13, 11.5, 10, 10.5, 12, 11, 1, 12, 13.5,
11.5, 9.5, 12, 11, 12, 8, 10.5, 12.5, 12, 12, 11.5, 12, 6, 10,
8.5, 10, 11, 11, 11.5, 12, 6.5, 14, 10.5, 13.5, 9.5, 12, 13,
12.5, 12.5, 12, 6.5, 12.5, 14, 11.5, 9.5, 10, 13.5, 13, 10.5,
10, 11, 13.5, 11.5, 10.5, 9, 10.5, 11.5, 12, 8.5, 12, 8.5, 11.5,
10.5, 12.5, 12.5, 11, 11, 11, 10.5), time = c(1.96, 4.94, 3.46,
4.94, 2.73, 6.22, 4.5, 2.67, 4.66, 3.83, 0.38, 2.6, 3.97, 4.18,
3.77, 3.44, 2.9, 3.93, 2.16, 3.51, 2.91, 3.19, 2.73, 6.36, 1.74,
4.39, 4.1, 2.26, 2.36, 5.32, 1.74, 3.66, 1.26, 5.61, 9.04, 4.61,
0.46, 3.98, 2.63, 5.5, 2.56, 5.92, 6.39, 2.26, 3.27, 7.95, 2.93,
4.93, 2.97, 2.43, 5.91, 3.07, 4.27, 3.21, 4.12, 4.72, 1.93, 0.69,
3.51, 4.39, 4.02, 3.18, 2.61, 4.61, 3.67, 0.54, 2.33, 2.93, 2.12,
1.06, 3.95, 2.31, 5.44, 0.17, 1.42, 0.55, 8.35, 2.53, 2.91, 3.26,
8.35, 2.26, 2.23, 7.18, 6.59, 6.36, 4.38, 7.67, 1.93, 3.34, 2.91,
8.54, 5.75, 3.77, 2.63, 0.97, 3.27, 1.58, 7.18, 2.08, 0.69, 5.43,
0.85, 2.26, 3.69, 3.18, 6.18, 2.93, 2.68, 0.69, 0.92, 2.34, 3.26,
0.85, 2.91, 4.3, 3.95, 7.67, 2.93, 2.1, 6.54, 6.31, 3.87, 2.91,
3.95, 3.35, 2.63, 1.49, 4.32, 3.51, 7.06, 2.67, 3.51, 3.46, 1.56,
4.33, 5.64, 2.73, 0.57, 2.87, 3.69, 2.56, 2.33, 4.27, 4.73, 4.02,
0.82, 4.11, 4.88, 2.29, 2.34, 3.72, 4.21, 1.49, 1.56, 3.03, 1.24,
2.65, 5.71, 1.67, 2.71, 1.49, 3.95, 4.51, 3.36, 5.21, 4.18, 4.54,
5.36, 4.25, 3.71, 0.95, 8.92, 3.12, 2.73, 1.36, 1.85, 7.24, 8.11,
2.2, 0.95, 5.16, 1.3, 6.54, 3.01, 1.97, 2.91, 3.26, 3.72, 1.79,
2.56, 1.96, 1.89, 1.89, 2.61, 5.25, 3.25, 5.26, 1.74, 3.73),
claylake = c(0, 0, 0, 0, 0, 17.53, 0.1, 0.59, 0, 9.13, 36.93,
12.75, 0, 0, 0, 0, 0, 0, 0, 0.09, 0.01, 0, 0, 9.43, 74.71,
26.42, 0.23, 0, 0, 35.27, 74.71, 0, 0, 0, 0, 0, 0, 0, 20.81,
9.46, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
1.14, 0, 26.42, 3.62, 0, 0, 0, 0.21, 0, 0, 0, 0.03, 10.43,
0.99, 3.6, 5.32, 0, 0.36, 0, 0, 0.25, 0.01, 0.22, 0, 0, 6.45,
0, 0, 0, 0, 0, 1.71, 0, 0, 0, 0, 0, 20.81, 0, 0, 0.18, 0,
0, 1.14, 0.03, 1.2, 0, 8.97, 0, 0, 0, 0, 1.14, 0, 1.56, 0.22,
1.2, 0, 0, 0.99, 0, 0, 0, 0, 4.14, 0, 0, 0.99, 0, 20.81,
0, 33.61, 0.09, 14.94, 0, 0, 0, 0, 0.41, 0, 2.7, 0, 0.61,
8.97, 0, 0, 0, 0, 1.7, 2.67, 7.71, 0.2, 8.63, 1.56, 0, 0.49,
0, 0, 0, 0, 0, 11.9, 33.08, 0, 0, 0.99, 2.13, 0, 0, 0, 0,
0.03, 0, 0, 0, 0, 0, 0, 2.86, 1.65, 0, 0, 0, 0, 0, 60.14,
0, 0, 0, 0, 0.22, 0, 0, 0, 0, 0.3, 0, 0, 0, 0, 0, 0, 5.57
), spinsandplain = c(81.94, 34.29, 89.55, 34.29, 80.86, 75.92,
81.55, 43.53, 97.3, 87.84, 60.62, 80.81, 11.73, 5.11, 98.67,
79.52, 60.73, 91.65, 2.82, 97.31, 73.65, 72.78, 96.51, 74.02,
25.09, 50.74, 96.62, 88.77, 98.8, 54.04, 25.09, 95.1, 69.85,
99.4, 78.79, 78.77, 48.16, 80.68, 75.79, 66.33, 68.3, 79.11,
91.89, 82.49, 98.33, 90.82, 91.24, 65.01, 69.24, 99.94, 99.75,
18.57, 90.39, 95.56, 71.07, 67.85, 92.37, 85.85, 17.89, 50.74,
79.65, 68.82, 74.05, 78.77, 87.67, 41.11, 91.74, 91.24, 44.8,
86.24, 97.7, 94.17, 85.59, 33.53, 85.23, 94.55, 78.52, 95.49,
73.65, 95.04, 78.52, 82.49, 77.26, 83.4, 98.29, 85.24, 98.78,
87.09, 81.36, 96.62, 3.4, 94.65, 28.6, 98.67, 75.79, 73.34,
98.33, 74.88, 83.4, 88.24, 85.85, 52.44, 95.84, 82.49, 62.11,
98.74, 70.32, 86.18, 95.67, 85.85, 11.42, 85.96, 75.53, 95.84,
95.46, 93.68, 97.7, 87.09, 91.24, 80.03, 87.77, 68.71, 17.51,
95.46, 97.7, 50.7, 75.79, 70.43, 61.06, 97.31, 74.63, 99,
17.89, 89.55, 99.25, 98.08, 97.61, 93.36, 99.03, 38.1, 62.11,
96.9, 88.87, 40.48, 90.21, 73.79, 95.2, 66.53, 96.67, 82.89,
85.96, 97.08, 75.74, 70.43, 99.25, 96.4, 98.88, 98.13, 85.32,
54.19, 99.2, 81.42, 97.7, 82.25, 97.42, 98.1, 5.11, 12.06,
66.14, 52.39, 52.72, 12.32, 87.32, 98.95, 71.55, 90.58, 97.9,
80.62, 93.32, 76, 86.48, 86.42, 39.54, 68.65, 6.05, 86.02,
3.4, 75.53, 97.08, 32.47, 68.3, 81.94, 89.64, 57.4, 74.05,
0.47, 96.76, 86.7, 78.46, 84.81)), row.names = c(5208L, 2828L,
1664L, 578L, 18L, 1644L, 4741L, 751L, 689L, 3813L, 1464L, 438L,
1553L, 4752L, 4960L, 376L, 2482L, 1811L, 5682L, 5441L, 4505L,
2281L, 2103L, 2993L, 562L, 4297L, 3592L, 5148L, 3793L, 1621L,
1912L, 1627L, 1737L, 4976L, 2173L, 5132L, 5758L, 2756L, 1789L,
5666L, 2628L, 2593L, 794L, 5779L, 5158L, 3123L, 4986L, 676L,
4200L, 2442L, 2751L, 4330L, 1802L, 2020L, 2500L, 1056L, 959L,
3290L, 4303L, 247L, 5586L, 922L, 1049L, 2432L, 2076L, 2560L,
1369L, 3636L, 3722L, 4137L, 1561L, 4915L, 2515L, 3034L, 5547L,
1491L, 1247L, 4116L, 455L, 4687L, 1697L, 5329L, 21L, 5724L, 3701L,
5697L, 2938L, 1721L, 61L, 998L, 4304L, 5798L, 651L, 910L, 2689L,
3986L, 2908L, 5753L, 2574L, 2345L, 1940L, 4317L, 4588L, 2179L,
665L, 4133L, 749L, 3977L, 3134L, 4190L, 3985L, 4937L, 2473L,
3238L, 4987L, 3915L, 4261L, 3521L, 2736L, 3665L, 1797L, 5692L,
5578L, 4087L, 2011L, 903L, 889L, 1523L, 3396L, 2291L, 5269L,
3644L, 3403L, 4814L, 4618L, 16L, 77L, 5385L, 2842L, 5816L, 2015L,
1443L, 3183L, 3331L, 4977L, 5380L, 989L, 4918L, 740L, 4637L,
887L, 1557L, 4295L, 4673L, 1918L, 5662L, 4167L, 1384L, 3441L,
614L, 2360L, 780L, 661L, 1267L, 2018L, 1906L, 3402L, 677L, 5218L,
2830L, 4979L, 3984L, 4924L, 1125L, 2640L, 986L, 1885L, 2573L,
5300L, 2398L, 4832L, 4816L, 3738L, 3276L, 3830L, 2425L, 2054L,
4273L, 5607L, 1678L, 378L, 1158L, 510L, 2210L, 2399L, 1952L,
2909L, 4945L, 2659L, 2642L), class = "data.frame")
我创建了一个新变量,因为燃烧年数的比例超出了 15 年(即二项式)
data$fireprop = cbind(data$YearsBurnt,data$YearsNotBurnt)
型号:
mfireprop = glmer(fireprop~log(time)+spinsandplain+rain15+claylake+rain15*log(time)+(1|circleID),na.action=na.fail, family=binomial, data=data)
趋势线代码:
d = ggpredict(mfireprop, terms = "time[exp]")
d = rename(d, "time" = x, "fireprop" = predicted)
ggplot(d, aes(time, fireprop)) +
geom_ribbon(aes(ymin = conf.low, ymax = conf.high), alpha = .1) +
geom_line(size = 2, colour = "black") +
theme_bw()
趋势线如下所示:
为什么x轴在数据停止的10小时不停止?为什么会涨到20,000?为什么 y 轴只达到 0.4?当部分比例为1时?
当我限制 x 轴和 y 轴时,它最终看起来像这样:
但是当我查看上面的原始数据时,趋势线似乎从一个非常奇怪的地方开始。
我不确定我做错了什么?
好的,我已经找到这里的主要问题了。在 ggpredict()
函数的文档中有一个名为 back.transform
的参数,默认为 TRUE。这意味着经过对数转换的数据将自动转换回原始响应范围。这就是为什么如果您检查 ggpredict 对象 d
,您会发现 time
变量实际上确实超过了该对象中的 8000。所以因为你没有标记 back.transform=FALSE
,而且还指定了 time[exp]
,发生的事情是函数自动对你的值求幂,然后你又做了一次。
如果我们查看记录的值:
summary(log(data$time))
Min. 1st Qu. Median Mean 3rd Qu. Max.
-1.7720 0.8154 1.1802 1.0904 1.4793 2.2017
然后我们对最大值取幂,得到之前的最大值:
exp(2.2017) # Exponentiated to get back to years
[1] 9.040369
summary(df$time) # The original variable
Min. 1st Qu. Median Mean 3rd Qu. Max.
0.170 2.260 3.255 3.519 4.390 9.040
如果我们再次取幂,我们最终会得到超过 8000 的最大时间。
exp(9.040369)
[1] 8436.89
所以,为了得到你想要的情节,你只需要在 ggpredict()
:
d = ggpredict(mfireprop, terms = "time")
d = rename(d, "time" = x, "fireprop" = predicted)
ggplot(d, aes(time, fireprop)) +
geom_ribbon(aes(ymin = conf.low, ymax = conf.high), alpha = .1) +
geom_line(size = 2, colour = "black") +
theme_bw()
时间被截断,因为在时间 0 没有变化。 YearsNotBurnt
始终为 0。因此,如果您从 ggpredict 查看对象 d
,您将在时间 0 的所有列中看到 NaN
。如果将模型简化为以下内容:
mfireprop2 = glmer(fireprop~
log(time) +
(1|circleID),
na.action=na.fail,
family=binomial,
data=df)
您将能够得到该图,但由于变化很小,置信区间将从 1 跨越到 0。我相信这是一个与分离有关的问题,基本上这意味着如果没有变化,或者如果某些东西可以完美地预测结果,二项式模型就不能适合常客模型。
我想我唯一想提的是你在评论中有一个关于“二项式 glm 中的非整数计数!”的问题。这是因为它期望因变量是试验的一部分,不应有小数。您的数据中有一些点似乎是半年的时间间隔。我对您的数据还不够熟悉,无法确定更好的选择是什么,但创建一个比例并在 weights=
参数中给出观察次数可能是一个选项。