将 mob() 树(partykit 包)与 nls() 模型一起使用

Using mob() trees (partykit package) with nls() model

我正在尝试将基于模型的递归分区 (MOB) 与 mob() 函数(来自 partykit 包)一起使用,以分离使用 nls() 功能。我必须定义我的模型并确定起始值。我一直在尝试查看是否可以将其与 mob() 函数一起使用,但无济于事。

我尝试按照第 7 页上的这个示例进行操作: https://cran.r-project.org/web/packages/partykit/vignettes/mob.pdf 我创建了一个拟合函数来估计起始值并将 return nls() 的估计值等。但在那之后我似乎什么也做不了。我想知道是否可以使用带有系数以及因变量和自变量的自定义模型,并将它们包含在 mob() 中并让它工作。我尝试了 lmtree() 函数,但当然这只会给出一条直线。

我的代码如下。基本上我使用分段线性回归来获取我正在使用的双指数曲线的起始值。这基本上是我得到的最远的。参数估计会给出错误等,即使你超过了它也不会 运行。我只需要知道 mob() 函数是否可以 运行 nls().

我加载了示例数据,但如果可以使用 nls()

    photo.try <- function(y, x,start = NULL, weights = NULL, offset = NULL, estfun = FALSE, object = TRUE) 
        {
            lin.mod1 <- lm(y ~ x)
            segmented.mod.2 <- segmented(lin.mod1, seg.Z = ~x, psi=1)
            segmented.mod1 <- segmented(lin.mod1, seg.Z = ~x, psi =  segmented.mod.2$psi[1,2])
            nls(y ~ (a*exp(-b * x) - c* exp(-d* x)), start = list(a = -1*(intercept(segmented.mod1)[[1]][1,1]) , b = slope(segmented.mod1)[[1]][1,1], 
            c = -1*(intercept(segmented.mod1)[[1]][2,1]), 
            d = -1*slope(segmented.mod1)[[1]][2,1]))

        }

photo_form <- Pn ~ (a*exp(-b * PAR) - c* exp(-d* PAR))| Species


photo_tree <- mob(photo_form, data = eco, fit = (photo.try))

这是我的示例数据:

eco <- structure(list(Species = structure(c(1L, 1L, 1L, 1L, 1L, 1L, 
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 
2L, 2L, 2L, 2L, 2L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 
3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 
3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 
3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 
3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 
3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 
3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 4L, 4L, 4L, 4L, 4L, 4L, 
4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 
4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 
4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 
4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 
4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 
4L, 4L, 4L, 4L, 4L, 4L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 
5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 
5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 
5L, 5L, 5L, 5L, 5L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 
6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 
6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 
6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 
6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 
6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 
6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 
6L, 6L, 6L, 7L, 7L, 7L, 7L, 7L, 7L, 7L, 7L, 7L, 7L, 7L, 7L, 7L, 
7L, 7L, 7L, 7L, 7L, 7L, 7L, 7L, 7L, 7L, 7L, 7L, 7L, 7L, 7L, 7L, 
7L, 7L, 7L, 7L, 7L, 7L, 7L, 7L, 7L, 7L, 7L, 7L, 7L, 7L, 7L, 7L, 
7L, 7L, 7L, 7L, 7L, 7L, 7L, 7L, 7L, 7L, 7L, 7L, 7L, 7L, 7L, 7L, 
7L, 7L, 7L, 7L, 7L, 7L, 7L, 7L, 7L, 7L, 7L, 7L, 7L, 7L, 7L, 7L, 
7L, 7L, 7L, 7L, 7L, 7L, 7L, 7L, 7L, 7L, 7L, 7L, 7L, 7L, 7L, 7L, 
7L, 7L, 7L, 7L, 7L, 7L, 7L, 7L, 7L, 7L, 7L, 7L, 7L, 7L), .Label = c("Bogum", 
"Clethra", "Eugene", "Guarria", "Melo", "Santa", "Sapium"), class = "factor"), 
    PAR = c(0, 58.6, 101.4, 228.6, 462.4, 904.7, 1565.8, 1992.1, 
    2395.9, 0, 72.8, 125.9, 232.8, 411, 841.1, 1669.6, 2394.5, 
    2394.9, 0, 53.5, 122.1, 231.6, 451, 808.5, 1575, 2394.6, 
    2395.1, 0, 70.9, 104.8, 251.1, 474.6, 858.3, 1612.3, 2393.3, 
    2395.1, 0, 63.1, 124.6, 277.1, 417.7, 824.4, 1649.6, 2377.7, 
    2381.9, 0, 31, 46.5, 115.7, 228.1, 424.3, 822.5, 1644.2, 
    2380.7, 2381.2, 0, 50.1, 118.1, 203.3, 413.2, 804.5, 1587.3, 
    2385.3, 0, 28.8, 36.9, 101.2, 211.7, 423.1, 793, 0, 43.6, 
    106.7, 200.8, 468.6, 808.4, 1567, 2367.1, 2376.5, 0.1, 40.4, 
    104.1, 202.2, 447.3, 794.7, 1546, 2391.8, 2393.3, 0.1, 44.1, 
    107.5, 227.4, 429.6, 802.5, 1668.4, 2391, 0, 42.2, 125.3, 
    126.2, 127.3, 240.3, 433.4, 791, 1600, 2396.8, 2397, 2399.3, 
    0, 72.7, 118.1, 236.9, 425, 828.4, 1613.3, 1615.4, 2396.1, 
    2396.5, 2397.2, 2397.5, 0, 62, 116.2, 235.5, 401.7, 879, 
    879.8, 1552.2, 1553.9, 2394.3, 2394.4, 2394.7, 2396.6, 0, 
    84.8, 135, 209.8, 425.3, 859.1, 1597.6, 2377.3, 2379.5, 2385.1, 
    0.1, 62, 106.3, 226.2, 442.9, 822.5, 1462.3, 2389.8, 2392.1, 
    0.1, 0.1, 73.9, 126, 249.8, 428.5, 846.5, 1555.3, 2390.1, 
    2390.7, 2390.8, 0, 68.7, 121.5, 209.7, 426.2, 803, 1525.9, 
    2389.8, 0, 52.8, 96.9, 211.1, 441.3, 787.9, 1566.5, 2415.2, 
    2415.3, 2415.5, 2417.5, 2417.7, 2418.5, 0.1, 46.5, 108.4, 
    233.5, 461.7, 792.3, 1635.7, 2415.1, 2415.6, 2415.6, 2416.5, 
    2416.6, 2417.8, 0.1, 68.3, 110, 239.5, 531.7, 847.2, 1591.4, 
    2387.3, 2387.6, 2389.7, 0, 49.7, 114.6, 230.6, 397.7, 398.2, 
    817.7, 1596.4, 2376.2, 2376.4, 2380.9, 0, 62.9, 65.5, 117, 
    209, 431.2, 854.5, 1611.3, 2387.3, 2388.5, 2390.3, 0, 49.1, 
    108.9, 200.3, 408.8, 842.2, 1630.2, 2386.5, 2386.8, 2388.2, 
    0, 64.8, 122.9, 226, 422.9, 801.6, 1635.7, 2383.6, 2383.6, 
    2384.3, 2386.1, 0, 36.7, 143.2, 213.7, 444.9, 814.9, 816.2, 
    1496.5, 2384.7, 2386.5, 2388.6, 0.1, 45.6, 105.2, 206.7, 
    494.8, 901.2, 1610.9, 2388, 2388.1, 2388.3, 2388.6, 0, 0.1, 
    45.9, 48.5, 100.2, 209.4, 432.4, 778, 1600.3, 2408.8, 2408.8, 
    0, 71.8, 121.6, 216.4, 404.3, 815.2, 1622, 2414.9, 2415.1, 
    2416.1, 2416.1, 0, 36.2, 97.5, 186.7, 417.9, 840.4, 1597.5, 
    2390.7, 2390.9, 2391.2, 2391.2, 2391.5, 2392.1, 2392.5, 0, 
    53.8, 138.2, 227, 403.6, 800.8, 1642.3, 2396.9, 2397.1, 0, 
    57.9, 95.1, 246.6, 466.8, 796.2, 1574.2, 2395.5, 2397.3, 
    0, 54.9, 94.9, 201.7, 408.1, 822.6, 1596, 2384.1, 0, 55.6, 
    131, 202.5, 419.8, 798.5, 1614, 2387.4, 2387.8, 0, 39.1, 
    109.6, 197.1, 403.3, 835.4, 836.9, 1725.9, 1727.4, 1729.3, 
    1730.6, 54.5, 58.6, 125.4, 226.9, 409, 806.8, 1578.8, 2377.2, 
    2380.1, 2388.3, 0, 68, 127.4, 206.9, 510.5, 814.9, 1561, 
    2404.1, 2404.8, 0, 58.4, 95.3, 229.6, 457.2, 781.5, 1634.4, 
    2399.8, 2401, 2403, 0.1, 56.5, 101.9, 221.8, 394.3, 815.1, 
    1655.4, 2411.8, 2411.9, 0, 50.2, 107.3, 220.5, 434.4, 819.8, 
    1630.6, 2412.4, 2412.6, 0, 48.4, 117.7, 195.3, 403.2, 801, 
    1632.7, 2388.9, 2389.3, 2390.7, 0, 50.4, 120.3, 234.7, 460.3, 
    829.1, 1581.7, 2398.5, 2402.3, 0, 60.8, 105.8, 215.8, 466.6, 
    826, 828.3, 1570.8, 2405.6, 2406.1, 2408.8, 0, 52.6, 106.9, 
    206.5, 414.3, 868.4, 1629.9, 1655.1, 2409.1, 2413, 0, 49.5, 
    100.6, 232.9, 389.4, 808.2, 1588.2, 2412.4, 2413.3, 2415.9, 
    0.1, 70.9, 110.5, 208.4, 409, 807.5, 1579.9, 2382.2, 2382.5, 
    2383.6, 2383.8, 0, 61.5, 106.5, 213.9, 473.8, 814.2, 1561.9, 
    2390.7, 2391.9, 2393.1, 0, 59.9, 64, 112, 216, 397.6, 807.4, 
    1625, 2392.3, 2395.1, 0, 74, 108.8, 109.7, 236.1, 433.6, 
    794.7, 1590.3, 2381.9, 2382.5, 0.1, 56.3, 114.5, 254.1, 487.7, 
    864.3, 1593.5, 2369.3, 2369.3, 2372.3, 2373.9, 0.2, 57.1, 
    110, 201.4, 402.7, 807.2, 1572.9, 2392.8, 2393.5, 0.1, 56.4, 
    122.5, 224.5, 420.2, 853.7, 1502.1, 2390.3, 2392.9, 0, 50.5, 
    53.7, 118.2, 230, 462.8, 794.3, 1513.4, 2391.4, 2392.3, 2393.4, 
    2393.4, 2394.1, 0.1, 49.7, 98.3, 208.3, 383.2, 850.7, 1653.5, 
    2395.3, 2396, 2397.1, 0, 48.4, 121.2, 228.8, 423.9, 817, 
    1708.5, 2389.9, 2389.9, 0, 66.4, 129.7, 209.4, 431.5, 794.1, 
    1673.7, 2383.7, 2384.2, 0, 57, 122.6, 215, 434.1, 838.5, 
    1657.5, 2386.4, 0.1, 22.6, 127.8, 220.4, 404.3, 810.9, 1592.3, 
    2386.7, 2388.7, 0, 49.8, 119.7, 200.5, 463.8, 828.7, 1560.7, 
    2384.5, 2385.7, 2391.2, 0, 73.1, 138.2, 226.6, 408.5, 815.3, 
    1627.3, 2390.2, 2395.4, 0, 61.2, 108.8, 233.8, 417.7, 824.5, 
    1502.7, 2395, 2396.2, 0, 56, 101.4, 226.3, 282.1, 412.9, 
    873.8, 1672.6, 2380.4, 2380.9, 2381.5, 0.1, 70.7, 138, 246, 
    444.4, 817.1, 1643.2, 2391.5, 2391.8, 2392), Pn = c(-0.95, 
    0.75, 0.94, 1.27, 1.5, 1.9, 2.14, 2.35, 2.38, 1.48, 3.51, 
    3.7, 3.99, 4.4, 4.32, 4.52, 4.73, 4.72, 1.97, 3.24, 4.23, 
    4.35, 4.41, 4.66, 4.57, 4.68, 4.88, 1.16, 3.64, 4.05, 4.75, 
    5.42, 5.57, 5.55, 5.89, 5.8, 1.48, 3.89, 4.7, 5.34, 5.47, 
    5.62, 5.71, 5.7, 6.08, 1.26, 0.59, 2.96, 4.34, 5, 4.82, 5.22, 
    5.2, 5.33, 5.51, 1.2, 2.95, 3.67, 3.9, 4.06, 4.59, 4.6, 4.62, 
    2.01, 1.92, 2.41, 2.19, 2.22, 2.41, 2.21, 1.6, 3.29, 3.97, 
    4.39, 4.89, 5.12, 4.93, 5.12, 5.1, 2.39, 3.84, 4.45, 4.63, 
    4.43, 4.93, 4.78, 4.73, 5.04, 3.09, 3.74, 4.03, 3.89, 4.52, 
    4.43, 4.24, 4.26, 1.5, 2.73, 2.83, 3.14, 2.89, 3.39, 2.89, 
    2.84, 3.34, 3.11, 3.16, 3.31, 0.1, 1.17, 1.72, 1.61, 1.64, 
    2.06, 2.17, 1.99, 2.31, 2.14, 2.27, 2.08, 0.17, 1.17, 1.32, 
    1.33, 1.4, 1.8, 1.48, 2, 1.81, 1.95, 2.09, 1.73, 1.85, 2.95, 
    4.33, 4.82, 4.98, 4.97, 5.03, 5.08, 5.22, 5.32, 4.88, 2.17, 
    3.08, 3.32, 3.42, 3.45, 3.67, 3.64, 3.71, 3.71, 2.85, 2.33, 
    3.15, 2.81, 3.22, 2.99, 3.16, 3.33, 3.56, 3.61, 3.63, 2.52, 
    3.55, 4.07, 4.1, 4.17, 4.41, 4.53, 4.56, 2.06, 2.57, 2.91, 
    2.61, 3.08, 3.29, 3.99, 6.49, 5.23, 6.08, 5.74, 4.41, 6.5, 
    1.59, 3.22, 3.59, 3.75, 3.84, 4.5, 4.93, 6.87, 6.75, 6.97, 
    6.53, 6.04, 6.82, 1.28, 3.56, 4.39, 5.27, 5.51, 6.38, 7.05, 
    7.46, 7.16, 7.24, 0.87, 2.45, 3.86, 4.32, 4.57, 4.43, 4.68, 
    4.71, 4.86, 4.36, 4.68, 1.06, 2.79, 4.05, 4.86, 5.48, 5.9, 
    6.38, 6.79, 7.46, 7.12, 7.03, 2.76, 3.92, 3.96, 4.07, 4.2, 
    4.5, 4.91, 5.52, 5.49, 5.33, 2.84, 4.78, 4.83, 4.76, 4.74, 
    4.84, 5.19, 5.59, 5.74, 5.7, 5.65, 3.02, 3.61, 4.14, 4.23, 
    4.45, 4.37, 4.5, 4.6, 4.78, 4.79, 4.85, 2.71, 4.26, 5.42, 
    6.24, 6.58, 6.63, 6.55, 7.29, 7.43, 7.24, 7, 3.36, 2.19, 
    2.86, 2.87, 2.37, 3.16, 2.68, 3, 3.4, 3.6, 4.35, 1.28, 2.62, 
    2.92, 3.3, 3.35, 3.58, 3.73, 4.02, 4, 3.7, 3.75, 1.61, 2.26, 
    2.5, 2.52, 2.71, 2.61, 2.75, 3.19, 2.92, 3.99, 4.36, 3.67, 
    4.14, 4.37, -0.28, 1.91, 2.78, 2.84, 2.96, 3.04, 3.24, 3.44, 
    3.58, 1.78, 4.12, 4.58, 4.33, 4.8, 4.7, 5.02, 5.09, 5.22, 
    2.79, 4.71, 4.89, 4.93, 4.87, 4.92, 4.83, 4.81, 1.66, 3, 
    4.04, 4.35, 4.56, 4.75, 4.75, 4.66, 4.89, 1.56, 2.77, 3.86, 
    3.58, 3.7, 3.76, 3.58, 4.55, 4.63, 4.05, 3.73, 1.76, 2.71, 
    2.98, 3.01, 3.06, 3.22, 2.99, 3.15, 3.32, 3.34, 1.58, 3.76, 
    4.97, 5.21, 5.29, 5.5, 5.59, 5.71, 5.74, 1.89, 2.67, 3.01, 
    3.14, 3.39, 3.57, 3.45, 3.91, 4.11, 3.94, 1.15, 2.88, 3.63, 
    4.32, 4.09, 4.43, 4.58, 4.61, 4.63, 1.23, 2.26, 3.15, 3.33, 
    3.3, 3.61, 3.46, 3.65, 3.67, 0.19, 2.23, 3.43, 4.1, 4.85, 
    5.21, 5.8, 6.27, 6.34, 6.08, 1.94, 3.72, 4.88, 5.51, 6.71, 
    6.51, 6.96, 7.01, 7.4, 0.48, 2.29, 2.5, 2.87, 3.18, 3.51, 
    3.13, 3.86, 4.13, 4.34, 4.03, 1.63, 3.64, 5.15, 5.95, 6.43, 
    6.57, 6.61, 6.51, 6.65, 6.56, 1.93, 3.95, 4.63, 5.66, 6.03, 
    6.28, 6.67, 6.69, 6.95, 6.75, 0.93, 3.14, 3.46, 3.9, 4.19, 
    4.27, 4.77, 5.39, 5.36, 5.24, 5.02, 1.71, 3.31, 3.86, 4.02, 
    4.02, 4.29, 4.36, 4.73, 4.88, 4.59, 1.63, 2.65, 2.63, 2.48, 
    2.93, 3.45, 4.01, 4.67, 5.02, 5.08, 1.93, 3.54, 3.8, 3.81, 
    4.04, 4.17, 4.38, 4.55, 4.99, 4.99, 1.29, 2.73, 3.32, 3.66, 
    3.77, 3.79, 4.14, 4.37, 4.22, 4.1, 4.14, 1.06, 2.89, 3.65, 
    4.01, 4.11, 4.19, 4.66, 5.03, 5.12, 0.97, 2.45, 2.99, 3.32, 
    3.34, 3.35, 3.47, 3.12, 3.38, 2.29, 1.72, 4.33, 5.49, 6.44, 
    6.96, 7.91, 7.49, 8.45, 8.21, 8.17, 8.71, 8.35, 0.29, 2.99, 
    3.93, 4.52, 5.69, 6.23, 6.23, 6.81, 6.96, 6.68, 0.99, 3.67, 
    4.62, 5.52, 5.86, 6.23, 5.91, 6.64, 6.29, -0.08, 3.34, 4.89, 
    6.02, 6.37, 6.59, 6.99, 6.95, 7.2, 0.99, 2.28, 2.72, 2.67, 
    2.99, 3.18, 3.55, 3.58, 1.31, 2.18, 5.55, 7.37, 8.42, 9.14, 
    9.44, 9.26, 9.5, 1.23, 3.11, 5.01, 6.21, 7.14, 7.44, 7.79, 
    7.73, 8.1, 7.96, 1.35, 3.33, 5.67, 6.58, 7.05, 7.36, 7.73, 
    7.75, 7.99, 0.4, 2.25, 2.83, 3.31, 3.55, 3.66, 3.96, 3.54, 
    3.77, 1.46, 2.91, 3.51, 3.64, 4.5, 3.83, 3.96, 4.17, 4.66, 
    4.09, 4.44, 2.41, 4.77, 5.49, 6.05, 6.15, 6.28, 6.6, 6.76, 
    6.75, 6.78)), .Names = c("Species", "PAR", "Pn"), class = "data.frame", row.names = c(NA, 
-628L))

是的,我们可以! ;-)

原则上,你试图做正确的事情,但有几个方面不太正确。主要问题是您如何传递数据和公式:由于 mob()nls() 指定其公式的方式一无所知,因此需要使用一个简单的公式 Pn ~ PAR | Species 然后fit 函数需要知道如何处理数据。 mob() 提供的预处理可以设置模型矩阵(截距、dummy/contrast 编码等)或模型框架(其中因子仍然是因子等)。在这种情况下,最简单的方法是使用默认模型矩阵,然后在拟合函数中省略截距。

您的代码的第二个问题是您使用了 fit 函数的扩展规范(带有 estfunobject 参数)但只提供了拟合模型对象。使用该规范 mob() 期望 fit 函数使用 coefficientsobjfun 等设置合适的列表

结合起来,这意味着您的 fit 函数应该如下所示:

photofit <- function(y, x = NULL, start = NULL, weights = NULL, offset = NULL, ...,
  estfun = FALSE, object = FALSE)
{
  ## only use first real regressor (without intercept)
  x <- x[, 2]

  ## obtain starting values if necessary
  if(is.null(start)) {
    aux_lm <- lm(y ~ x)
    aux_seg_2 <- segmented::segmented(aux_lm, seg.Z = ~ x, psi = 1)
    aux_seg_1 <- segmented::segmented(aux_lm, seg.Z = ~ x, psi = aux_seg_2$psi[1, 2])
    start <- list(
      a = -1 * (segmented::intercept(aux_seg_1)[[1]][1, 1]),
      b =           segmented::slope(aux_seg_1)[[1]][1, 1], 
      c = -1 * (segmented::intercept(aux_seg_1)[[1]][2, 1]), 
      d = -1 *      segmented::slope(aux_seg_1)[[1]][2, 1]
    )
  } else {
    start <- as.list(start)
  }

  ## estimate NLS model
  rval <- nls(y ~ (a * exp(-b * x) - c * exp(-d * x)), start = start)

  ## return processed information for mob()
  list(
    coefficients = coef(rval),
    objfun = deviance(rval),
    estfun = if(estfun) sandwich::estfun(rval) else NULL,
    object = if(object) rval else NULL    
  )
}

然后你就可以种植MOB树了。指定 verbose = TRUE 控制选项将在您等待时为您提供一些进度信息:

photomob <- mob(Pn ~ PAR | Species, data = eco, fit = photofit,
  control = mob_control(verbose = TRUE))
coef(photomob)
##           a             b         c             d
## 4  2.967680 -3.216708e-05  1.519680  1.076879e+01
## 5 -1.811596  1.967366e-02 -3.573079 -4.877852e-05
## 6 -2.772783  1.438087e-02 -4.177953 -7.821814e-05
## 8 -2.427253  1.757744e-02 -4.449105 -1.328930e-04
## 9 -4.579248  1.020021e-02 -5.714575 -7.502393e-05

然后您还可以想象这棵树。默认情况下,每个节点中都会显示一个数字摘要,但您也可以轻松显示拟合曲线:

plot(photomob)
plot(photomob, terminal_panel = node_bivplot, tnex = 2)

如您所见,树选择了五个具有不同参数的终端节点。我建议您对适合不同节点的模型进行更多诊断,因为我不确定所有参数的识别效果如何。我对 NLS 不是很熟悉,可能完全错了,但似乎并非所有参数都能可靠地确定。

作为一个示例,我执行以下操作:我从树中提取所有九个合适的 nls 对象。对于来自根节点(节点 1)的模型,我通过对所有观察方向的梯度贡献求和来计算梯度(由 estfun() 方法计算):

photonls <- refit.modelparty(photomob)
library("sandwich")
colSums(estfun(photonls[[1]]))
##             a             b             c             d 
##  2.010552e-05  5.753230e-02 -1.166331e-04  6.771585e+00 

参数 a-c 的梯度相当接近于零,但对于 d 则不是。这也可能会影响 mob() 中基于观察梯度贡献(也称为模型分数或估计函数)的推理。

总之:你想做的,都能做到!但我建议考虑一个更简单的模型。如果你这样做,你只需要相应地修改 photofit() 功能,然后 运行 再次通过 mob() 它。