Получите пиковую скорость и возраст для конкретного объекта при пиковых значениях скорости из линейных смешанных сплайновых моделей

Я подгоняю линейную модель смешанных эффектов с естественной функцией сплайна для возраста, чтобы описать нелинейную траекторию для повторяющегося результата y (содержание минералов в костях в граммах) во времени (age в годах).

Как я могу решить производные сплайна, чтобы получить кривую скорости и оценить для каждого человека их пиковую скорость (граммы / год) и возраст при максимальной скорости (годы) по этой модели?

Это пример данных

dat <- structure(list(id = c(1001L, 1001L, 1001L, 1001L, 1001L, 1002L, 
    1003L, 1004L, 1004L, 1004L, 1004L, 1004L, 1004L, 1004L, 1005L, 
    1005L, 1005L, 1005L, 1005L, 1006L, 1006L, 1006L, 1006L, 1006L, 
    1007L, 1007L, 1008L, 1008L, 1008L, 1008L, 1008L, 1009L, 1009L, 
    1009L, 1010L, 1010L, 1010L, 1011L, 1012L, 1012L, 1012L, 1013L, 
    1013L, 1014L, 1015L, 1015L, 1015L, 1016L, 1016L, 1016L, 1016L, 
    1016L, 1017L, 1017L, 1018L, 1020L, 1020L, 1021L, 1021L, 1021L, 
    1021L, 1022L, 1022L, 1023L, 1023L, 1023L, 1023L, 1023L, 1023L, 
    1023L, 1023L, 1023L, 1023L, 1024L, 1024L, 1024L, 1024L, 1024L, 
    1025L, 1025L, 1025L, 1026L, 1026L, 1026L, 1026L, 1027L, 1027L, 
    1028L, 1028L, 1028L, 1028L, 1028L, 1028L, 1028L, 1029L, 1029L, 
    1029L, 1029L, 1029L, 1029L, 1030L, 1030L, 1030L, 1030L, 1030L, 
    1030L, 1030L, 1030L, 1031L, 1031L, 1031L, 1031L, 1032L, 1032L, 
    1032L, 1032L, 1032L, 1033L, 1033L, 1033L, 1033L, 1034L, 1034L, 
    1034L, 1034L, 1034L, 1035L, 1035L, 1036L, 1037L, 1037L, 1037L, 
    1037L, 1039L, 1039L, 1040L, 1040L, 1040L, 1040L, 1040L, 1040L, 
    1041L, 1041L, 1041L, 1041L, 1041L, 1041L, 1042L, 1042L, 1042L, 
    1042L, 1042L, 1042L, 1042L, 1043L, 1043L, 1043L, 1043L, 1044L, 
    1044L, 1044L, 1045L, 1045L, 1045L, 1045L, 1045L, 1045L, 1047L, 
    1048L, 1048L, 1049L, 1049L, 1049L, 1049L, 1051L, 1051L, 1052L, 
    1052L, 1052L, 1052L, 1052L, 1053L, 1053L, 1053L, 1053L, 1053L, 
    1054L, 1054L, 1054L, 1054L, 1054L, 1054L, 1054L, 1054L, 1056L, 
    1056L, 1056L, 1056L, 1057L, 1057L, 1058L, 1058L, 1058L, 1058L, 
    1058L, 1060L, 1060L, 1060L, 1061L, 1061L, 1061L, 1061L, 1061L, 
    1062L, 1062L, 1062L, 1062L, 1062L, 1063L, 1063L, 1063L, 1064L, 
    1064L, 1064L, 1064L, 1065L, 1065L, 1066L, 1066L, 1066L, 1066L, 
    1066L, 1066L, 1067L, 1067L, 1067L, 1068L, 1068L, 1068L, 1068L, 
    1068L, 1068L, 1068L, 1069L, 1070L, 1070L, 1070L, 1071L, 1071L, 
    1071L, 1072L, 1072L, 1072L, 1072L, 1072L, 1073L, 1073L, 1073L, 
    1073L, 1074L, 1074L, 1074L, 1075L, 1075L, 1075L, 1075L, 1075L, 
    1075L, 1076L, 1076L, 1076L, 1077L, 1077L, 1077L, 1077L, 1077L, 
    1077L, 1078L, 1078L, 1078L, 1078L, 1078L, 1078L, 1078L, 1080L, 
    1080L, 1080L, 1080L, 1081L, 1081L, 1082L, 1082L, 1082L, 1083L, 
    1083L, 1084L, 1085L, 1085L, 1085L, 1085L, 1085L, 1085L, 1086L, 
    1086L, 1086L, 1087L, 1087L, 1087L, 1087L, 1087L, 1087L, 1087L, 
    1087L, 1088L, 1088L, 1088L, 1088L, 1089L, 1089L, 1089L, 1089L, 
    1089L, 1090L, 1090L, 1091L, 1091L, 1091L, 1091L, 1091L, 1092L, 
    1092L, 1092L, 1092L, 1092L, 1093L, 1093L, 1093L, 1093L, 1094L, 
    1094L, 1094L, 1094L, 1094L, 1095L, 1095L, 1095L, 1095L, 1096L, 
    1097L, 1097L, 1098L, 1098L, 1098L, 1098L, 1098L, 1099L, 1099L, 
    1099L, 1099L, 1099L, 1099L, 1099L, 1099L, 1100L, 1100L, 1100L, 
    1101L, 1101L, 1101L, 1101L, 1103L, 1103L, 1103L, 1103L, 1103L, 
    1103L, 1103L, 1104L, 1104L, 1104L, 1104L, 1105L, 1105L, 1105L, 
    1106L, 1106L, 1106L, 1106L, 1106L, 1106L, 1106L, 1106L, 1106L, 
    1107L, 1108L, 1110L, 1111L, 1112L, 1117L, 1123L), y = c(1934.047646, 
    1075.598345, 1956.214821, 2000.38538, 2000.38538, 732.315937, 
    3119.86, 624.951231, 791.2764892, 1884.530826, 624.951231, 1047.57, 
    1047.57, 791.2764892, 1238.306103, 1555.042976, 2547.870529, 
    2547.870529, 2467.385, 1181.635212, 1181.635212, 565.306282, 
    2016.027874, 2016.027874, 712.6134567, 635.2537841, 2167.362267, 
    2575.574188, 2167.362267, 2480.028259, 2575.574188, 2875.363243, 
    1180.139938, 2828.037147, 3017.119362, 2722.940933, 2167.92, 
    2409.652458, 2245.442558, 724.1520328, 635.6034756, 1649.08326, 
    966.8182507, 865.2717723, 1570.23, 916.1300105, 1180.999973, 
    2351.32885, 2418.851707, 2290.038887, 2224.060562, 2509.52, 1174.589081, 
    1540.219376, 2692.26, 1300.899734, 1100.650177, 1786.628242, 
    1705.842979, 543.8596134, 1786.628242, 2115.374241, 2331.46, 
    875.949604, 2241.945103, 2319.666939, 2316.220234, 719.7139549, 
    2042.803307, 719.7139549, 1132.977503, 875.949604, 2316.220234, 
    1737.18, 1351.629826, 1291.44593, 1291.44593, 1108.26586, 1028.979719, 
    1291.44593, 2068.934227, 2440.784416, 1036.72, 894.6663704, 2449.184731, 
    1109.9, 672.9310664, 2072.320354, 2114.215416, 2114.215416, 1805.422001, 
    2461.18, 2101.374248, 2105.879, 1600.086481, 2866.84, 1600.086481, 
    2807.311, 3055.569931, 1600.086481, 2602.287521, 2690.007614, 
    620.5975037, 2608.4, 2722.3, 2713.66185, 2608.4, 1590.002, 2198.211, 
    2488.097725, 2198.211, 2322.616348, 2627.1, 2418.328346, 2601.661034, 
    531.7369251, 811.9494571, 884.31, 768.0526981, 652.1271248, 768.0526981, 
    2767.479, 1047.144354, 1047.144354, 1995.119, 1995.119, 707.6093158, 
    707.6093158, 1120.650104, 3036.591904, 3036.591904, 3081.86, 
    1193.583691, 2056.569244, 1823.155, 1238.948124, 2124.685, 887.20438, 
    1823.155, 2056.569244, 2056.569244, 2560.155342, 3095.923164, 
    3095.923164, 3003.729011, 2861.12, 2560.155342, 2735.26, 822.8209591, 
    1648.951, 1648.951, 1648.951, 822.8209591, 906.7692623, 582.787096, 
    1286.45, 797.2365359, 2566.770554, 2666.41, 2666.41, 2045.320816, 
    2401.21, 2401.21, 2583.2, 2581.32, 2622.357, 2581.32, 2588.462498, 
    442.433671, 1251.627064, 406.2565479, 2108.787437, 983.1101169, 
    2102.085403, 1155.713411, 1909.797131, 2871.55, 2711.07, 2883.22245, 
    2883.22245, 2711.07, 3027.103172, 3108.21537, 3007.87294, 3208.963631, 
    3108.21537, 2617.91, 2457.464466, 2890.51, 2698.48214, 2700.723, 
    2700.723, 2817.668579, 2700.723, 1349.90691, 1476.19994, 1552.95, 
    1349.90691, 925.8325004, 1258.28, 840.1875095, 2405.175911, 840.1875095, 
    1056.678543, 1571.936, 1210.89, 1210.89, 673.7005405, 687.7842464, 
    1016.86, 1217.866, 1493.791817, 2246.726913, 1054.821, 1054.821, 
    563.6580887, 1054.821, 1540.429863, 2209.006493, 1437.835186, 
    2191.308, 1412.128944, 2724.164597, 2791.705185, 2727.774208, 
    2070.451198, 866.7974147, 1661.082638, 2108.271309, 2411.515434, 
    2342.026085, 2071.06, 2258.321014, 1537.06, 760.6319065, 867.7596569, 
    1907.60466, 1770.658, 760.6319065, 912.8781966, 912.8781966, 
    912.8781966, 1257.222706, 2586.922356, 1608.28, 962.5674305, 
    1085.451181, 2539.218132, 2535.526085, 2561.60054, 1600.198, 
    2100.048149, 758.3851737, 758.3851737, 2643.373329, 367.7795143, 
    866.0683727, 718.5049658, 866.0683727, 1906.694649, 2291.48, 
    2190.560314, 744.1710777, 1498.981777, 2460.912292, 590.1345787, 
    2487.559135, 1855.601353, 660.9104843, 1116.08, 792.929533, 708.8373737, 
    2272.232933, 1801.729801, 2299.800095, 2272.232933, 2299.800095, 
    1895.828438, 1757.75, 1050.279345, 1757.75, 1326.09478, 1326.09478, 
    1633.119305, 1558, 1167.971405, 1828.16, 1788.571758, 2175.469, 
    1071.039494, 941.6030864, 2053.067215, 1461.02132, 1597.646778, 
    1885.321567, 2195.704372, 2195.704372, 1675.768558, 3157.550789, 
    1565.173126, 2195.704372, 3157.550789, 2404.836883, 2541.045593, 
    585.7223682, 2465.177761, 2678.462074, 500.3733997, 2465.177761, 
    781.342, 898.3551559, 2465.177761, 2465.177761, 1807.02, 1418.888027, 
    1797.36, 1807.02, 2200.06, 2218.369926, 2200.06, 1986.642735, 
    2088.292, 2069.139, 1507.901432, 2061.395798, 2075.164864, 2081.913219, 
    2081.913219, 483.8579493, 1857.88, 2578.772636, 1857.88, 1857.88, 
    1039.632153, 2288.28, 2288.28, 1831.349922, 2349.23, 933.1002788, 
    2626.298935, 1521.744, 933.1002788, 2626.298935, 1984.760715, 
    2450.333, 1732.339031, 1984.760715, 2731.9, 869.2320918, 1785.72, 
    1922.798, 3081.28, 1508.8, 2421.288597, 1922.798, 1268.074959, 
    1569.05, 1808.115, 1569.05, 1268.074959, 2165.724808, 2165.724808, 
    1808.115, 2084.149837, 2693.027184, 2464.489, 2607.653496, 1012.837271, 
    1012.837271, 2673.190872, 2635.290516, 2773.42, 2635.290516, 
    2654.772674, 2377.905655, 2679.014969, 2654.772674, 1226.40016, 
    1470.69, 1273.789799, 2294.926086, 1226.40016, 1470.69, 1273.789799, 
    1873.817, 2274.930534, 2317.429165, 959.1709613, 1328.159428, 
    1328.159428, 1328.159428, 959.1709613, 1630.28, 1610.54982, 2507.05302, 
    750.467966, 750.467966, 821.2255058, 802.8240452, 2829.47879), 
        age = c(31.54004107, 11.95071869, 27.88501027, 27.88501027, 
        25.07871321, 10.90759754, 25.70020534, 9.560574949, 11.17864476, 
        15.8384668, 9.560574949, 11.23613963, 14.01232033, 10.54620123, 
        12.89527721, 14.52977413, 24.96919918, 24.72005476, 23.95893224, 
        13.31690623, 11.52087611, 9.927446954, 22.10814511, 16.44353183, 
        10.90759754, 7.991786448, 17.26488706, 23.95893224, 15.66872005, 
        17.63723477, 24.72005476, 30.97330595, 11.52087611, 17.5633128, 
        30.11088296, 23.31279945, 17.26488706, 20.58590007, 28.27926078, 
        11.66324435, 9.927446954, 13.92744695, 11.20328542, 12.70362765, 
        13.52498289, 12.21355236, 13.80150582, 22.81724846, 39.3045859, 
        16.62696783, 22.63107461, 29.86447639, 12.54483231, 14.42299795, 
        34.27789185, 12.91170431, 12.25462012, 21.81245722, 21.81245722, 
        10.05065024, 23.6659822, 16.22450376, 28.74743326, 12.70362765, 
        35.43052704, 21.21013005, 19.28542094, 12.77207392, 16.59411362, 
        12.12867899, 11.29637235, 11.81930185, 19.04449008, 19.93429158, 
        16.14236824, 12.85420945, 13.21560575, 11.61396304, 11.85763176, 
        13.3798768, 17.42915811, 24.41341547, 13.08418891, 11.6659822, 
        24.41341547, 12.06297057, 10.22861054, 26.15468857, 21.71937029, 
        20.1889117, 12.60232717, 25.39904175, 30.72689938, 19.22245038, 
        14.45037645, 24.77207392, 13.47570157, 17.87816564, 27.52635181, 
        15.16221766, 19.68514716, 21.67282683, 9.062286105, 20.43805613, 
        21.67282683, 21.24024641, 20.70362765, 13.5687885, 17.13347023, 
        28.11498973, 24.16974675, 18.19575633, 27.73442847, 15.52361396, 
        20.70362765, 11.76728268, 10.98699521, 11.51540041, 9.902806297, 
        13.05407255, 8.703627652, 25.60164271, 10.59000684, 10.59000684, 
        14.45859001, 14.05886379, 10.88295688, 10.75427789, 10.59000684, 
        26.50513347, 18.83093771, 22.86379192, 11.8384668, 15.04449008, 
        15.42505133, 14.14099932, 28.06844627, 11.51540041, 14.66119097, 
        13.79055441, 15.37850787, 22.58179329, 22.86379192, 30.0752909, 
        21.85900068, 25.60164271, 15.29089665, 26.79534565, 11.68514716, 
        15.42505133, 15.58384668, 15.08555784, 14.11909651, 11.6659822, 
        10.21765914, 12.1670089, 10.50239562, 23.3045859, 15.92607803, 
        22.58179329, 16.65982204, 20.58590007, 39.3045859, 32.56947296, 
        16.90349076, 25.12799452, 17.88364134, 19.46338125, 8.736481862, 
        14.14099932, 8.736481862, 17.68104038, 14.54893908, 19.22245038, 
        12.98562628, 22.45311431, 18.83093771, 38.68856947, 26.50513347, 
        25.44010951, 28.70910335, 19.21697467, 30.0752909, 26.50513347, 
        29.45106092, 33.31690623, 16.68172485, 15.816564, 24.89801506, 
        15.816564, 18.7761807, 18.4366872, 19.45790554, 19.78370979, 
        14.98973306, 15.89869952, 29.06502396, 16.14236824, 10.74880219, 
        13.47843943, 10.5982204, 24.61875428, 10.74880219, 12.47364819, 
        16.95277207, 12.41889117, 13.44832307, 9.984941821, 9.451060917, 
        12.59137577, 13.38261465, 15.14852841, 21.65913758, 12.57494867, 
        12.40520192, 10.75701574, 15.16495551, 15.67419576, 22.52703628, 
        13.31143053, 16.71457906, 12.98288843, 32.16974675, 25.3798768, 
        30.57084189, 22.14647502, 11.43874059, 13.25119781, 18.48049281, 
        25.81519507, 24.78028747, 17.85626283, 27.70704997, 13.28952772, 
        8.703627652, 11.61396304, 35.04996578, 15.61943874, 8.703627652, 
        13.33333333, 10.56810404, 11.34017796, 13.5797399, 28.79671458, 
        12.56673511, 13.33333333, 12.55578371, 30.80082136, 23.63039014, 
        29.66461328, 13.25119781, 17.46748802, 8.703627652, 8.703627652, 
        21.21013005, 9.768651608, 13.46748802, 10.75427789, 13.24298426, 
        26.87474333, 27.43326489, 20.6899384, 10.0752909, 13.37713895, 
        28.38056126, 8.911704312, 24.62149213, 14.32443532, 10.24229979, 
        13.87268994, 10.54620123, 11.44421629, 21.68377823, 15.61943874, 
        27.97809719, 28.90075291, 28.90075291, 24.64339493, 14.32443532, 
        10.61190965, 15.8110883, 14.25051335, 14.25051335, 13.64818617, 
        26.05338809, 13.69746749, 23.98083504, 16.68172485, 20.42162902, 
        12.68172485, 11.51813826, 16.65982204, 14.32443532, 15.49897331, 
        35.04996578, 18.70225873, 17.47570157, 14.66666667, 26.83915127, 
        13.29226557, 18.14647502, 25.70020534, 14.67761807, 16.61601643, 
        9.812457221, 15.96714579, 24.41341547, 8.911704312, 17.61806982, 
        11.87953457, 11.80561259, 19.15400411, 17.61806982, 15.70704997, 
        12.35318275, 18.12457221, 16.8733744, 32.02464066, 32.02464066, 
        25.30047912, 16.13415469, 19.37850787, 26.50513347, 15.89869952, 
        13.79055441, 25.42368241, 16.05201916, 15.43874059, 9.158110883, 
        14.39014374, 22.12183436, 15.70704997, 15.35934292, 11.44421629, 
        28.45995893, 17.06502396, 14.39014374, 26.32991102, 12.38056126, 
        16.42436687, 13.37713895, 11.70978782, 17.62628337, 16.13415469, 
        17.61806982, 15.11019849, 14.09993155, 21.89185489, 13.80150582, 
        16.8733744, 17.73305955, 25.55509925, 14.75975359, 24.03559206, 
        14.36002738, 12.73100616, 16.09034908, 18.12457221, 15.11019849, 
        13.69472964, 23.03901437, 16.94182067, 15.70704997, 13.99315537, 
        21.89185489, 15.65776865, 19.25530459, 10.43394935, 12.72826831, 
        24.41341547, 24.25735797, 37.41820671, 37.41820671, 25.25393566, 
        24.78028747, 25.25393566, 37.41820671, 12.11772758, 14.19575633, 
        14.091718, 15.10746064, 13.16906229, 12.09856263, 13.3798768, 
        14.39014374, 36.3504449, 22.68035592, 11.21149897, 12.73100616, 
        13.34702259, 14.5982204, 11.31827515, 15.14579055, 15.44969199, 
        15.65776865, 12.12867899, 12.43531828, 12.72005476, 14.11909651, 
        24.25735797)), row.names = c(7L, 303L, 323L, 372L, 391L, 
    240L, 311L, 38L, 46L, 94L, 149L, 154L, 185L, 362L, 40L, 70L, 
    98L, 262L, 305L, 73L, 74L, 77L, 306L, 374L, 104L, 397L, 14L, 
    43L, 188L, 248L, 370L, 50L, 101L, 143L, 25L, 155L, 251L, 37L, 
    173L, 208L, 263L, 49L, 383L, 389L, 30L, 237L, 353L, 156L, 283L, 
    288L, 302L, 325L, 33L, 158L, 159L, 35L, 360L, 57L, 128L, 204L, 
    387L, 300L, 365L, 16L, 51L, 82L, 85L, 93L, 148L, 150L, 232L, 
    242L, 287L, 32L, 62L, 200L, 285L, 290L, 193L, 352L, 398L, 54L, 
    175L, 203L, 324L, 69L, 195L, 92L, 106L, 141L, 189L, 218L, 347L, 
    394L, 23L, 24L, 120L, 166L, 257L, 349L, 6L, 118L, 235L, 266L, 
    269L, 275L, 282L, 390L, 122L, 153L, 330L, 378L, 53L, 88L, 229L, 
    241L, 314L, 135L, 278L, 332L, 384L, 64L, 168L, 207L, 212L, 359L, 
    329L, 338L, 130L, 67L, 108L, 286L, 316L, 182L, 254L, 113L, 215L, 
    247L, 273L, 322L, 336L, 27L, 102L, 162L, 171L, 270L, 326L, 19L, 
    205L, 210L, 307L, 333L, 358L, 375L, 41L, 111L, 179L, 226L, 2L, 
    277L, 367L, 68L, 83L, 147L, 180L, 260L, 354L, 144L, 81L, 342L, 
    103L, 217L, 321L, 376L, 131L, 280L, 39L, 267L, 291L, 301L, 400L, 
    11L, 36L, 152L, 177L, 377L, 21L, 201L, 236L, 281L, 312L, 331L, 
    355L, 369L, 8L, 176L, 202L, 385L, 45L, 327L, 12L, 138L, 151L, 
    157L, 233L, 95L, 258L, 279L, 224L, 239L, 243L, 310L, 328L, 63L, 
    191L, 214L, 227L, 356L, 80L, 110L, 366L, 97L, 107L, 293L, 373L, 
    117L, 335L, 22L, 160L, 209L, 221L, 230L, 268L, 55L, 163L, 284L, 
    5L, 10L, 76L, 132L, 222L, 256L, 399L, 228L, 127L, 343L, 357L, 
    133L, 259L, 334L, 261L, 341L, 382L, 393L, 395L, 213L, 219L, 249L, 
    289L, 44L, 126L, 368L, 42L, 72L, 196L, 297L, 308L, 320L, 84L, 
    137L, 172L, 60L, 129L, 142L, 186L, 197L, 319L, 15L, 109L, 115L, 
    116L, 125L, 199L, 223L, 190L, 245L, 346L, 396L, 146L, 364L, 1L, 
    29L, 192L, 112L, 170L, 315L, 164L, 225L, 231L, 255L, 274L, 345L, 
    65L, 96L, 264L, 4L, 28L, 31L, 59L, 87L, 250L, 271L, 295L, 161L, 
    198L, 265L, 339L, 18L, 26L, 114L, 124L, 174L, 145L, 304L, 105L, 
    119L, 140L, 238L, 381L, 48L, 52L, 71L, 351L, 371L, 244L, 253L, 
    294L, 340L, 20L, 75L, 86L, 165L, 167L, 47L, 89L, 298L, 318L, 
    211L, 350L, 380L, 66L, 79L, 90L, 234L, 309L, 61L, 99L, 139L, 
    276L, 299L, 344L, 348L, 361L, 313L, 337L, 379L, 9L, 58L, 181L, 
    187L, 17L, 100L, 121L, 123L, 184L, 206L, 220L, 178L, 292L, 386L, 
    392L, 194L, 252L, 272L, 3L, 56L, 134L, 136L, 183L, 216L, 246L, 
    296L, 363L, 169L, 388L, 78L, 34L, 13L, 91L, 317L), class = "data.frame")

Это код для подгонки модели и построения средней прогнозируемой траектории.

    library(nlme)
    library(splines)
    library(tidyverse)

# LINEAR MIXED MODEL WITH NATURAL SPLINE FUNCTION FOR AGE
    nspline_model <- lme(y ~ ns(age, df = 4), data = dat, random = ~ age | id)
    
# PLOT MEAN FITTED NATURAL SPLINE TRAJECTORY
    dat$ns_pred <- fitted(nspline_model, level = 0)
    ggplot(data = dat, aes(x = age, y = ns_pred)) + geom_line()

Это возможное решение с использованием функции getPeakTrough и smooth.spline, но я не уверен, правильно ли оно - или как получить конкретные пиковые скорости объекта, а не только средние значения.

    # POSSIBLE SOLUTION FOR THE NATURAL SPLINE MODEL???
    x <- dat$age
    y <- fitted(nspline_model, level = 0)
    y <- predict(smooth.spline(x, y, df = 4), x, deriv = 1)$y
    sitar::getPeakTrough(x, y, peak = T)
    # mean age at peak velocity = 11.82447 years; 
    # mean peak velocity 187.13404 grams/year

r spline mixed-models cubic-spline bspline
person aelhak    schedule 15.12.2020    source источник

Ответы (2)


arrow_upward
3
arrow_downward

Примечание. Мой ответ здесь также служит учебным пособием / виньеткой для SplinesUtils.

Пакет SplinesUtils (https://github.com/ZheyuanLi/SplinesUtils) теперь расширен для поддержки линейных смешанных моделей, приспособленных nlme::lme. (Предыдущая версия мотивирована этим потоком SO: Экспорт подогнанных сплайнов регрессии (построенных с помощью 'bs' или 'ns') как кусочных полиномов < / а>).

library(SplinesUtils)

Описание модели

nspline_model <- lme(y ~ ns(age, df = 4), data = dat, random = ~ age | id)

Это линейная модель смешанного эффекта с:

  • фиксированный эффект ns(age, df = 4), который является естественным кубическим сплайном age:

  • случайный эффект ~ age | id, который является линейной функцией age для каждого id.

Фиксированный эффект отражает среднее соотношение между y и age по населению; в то время как случайный эффект - это специфическое для выборки отклонение от этой общей зависимости.

Используйте SplinesUtils для обработки общей взаимосвязи

spl_population <- RegSplineAsPiecePoly(nspline_model, "ns(age)")
#Error: 
#  Required SplineTerm not found! Available terms are:
#    * ns(age, df = 4)

Выше я намеренно написал строку неправильного кода. Эта ошибка сбила с толку некоторых пользователей (Определение коэффициента сплайн-регрессии в множественной регрессии в R). Проблема в том, что пользователи должны правильно указать имя сплайна, которое внутренне известно функции подгонки модели. Начиная с версии 0.2, сообщение об ошибке лучше отформатировано, перечисляя такие внутренние имена по одному в каждой строке. По сути, пользователям просто нужно скопировать правильный, а затем вставить его в RegSplineAsPiecePoly:

spl_population <- RegSplineAsPiecePoly(nspline_model, "ns(age, df = 4)")

На данный момент репараметризация завершена. Существуют различные способы проверки экспортированных кусочно-полиномов.

Самый принципиальный из них - получить коэффициенты в матричной форме:

spl_population$PiecePoly$coef
#             [,1]       [,2]        [,3]          [,4]
#[1,] 5.684342e-13 645.727356 1383.877112  1.871766e+03
#[2,] 9.419308e+01 220.369323  209.951537 -3.223616e-01
#[3,] 1.782352e-12  26.623843  -30.064255 -4.796068e-01
#[4,] 1.872590e+00  -6.240304    1.432464  9.595289e-03

Чтобы лучше понять значение этих чисел, мы можем распечатать отформатированные строки полиномиальных уравнений:

spl_population  ## or print(spl_polulation)
#4 piecewise polynomials of degree 3 are constructed!
#Use 'summary' to export all of them.
#The first 4 are printed below.
#5.68e-13 + 94.2 * (x - 7.99) + 1.78e-12 * (x - 7.99) ^ 2 + 1.87 * (x - 7.99) ^ 3
#646 + 220 * (x - 12.7) + 26.6 * (x - 12.7) ^ 2 - 6.24 * (x - 12.7) ^ 3
#1.38e+03 + 210 * (x - 15.8) - 30.1 * (x - 15.8) ^ 2 + 1.43 * (x - 15.8) ^ 3
#1.87e+03 - 0.322 * (x - 22.6) - 0.48 * (x - 22.6) ^ 2 + 0.0096 * (x - 22.6) ^ 3

Форматированные строки всегда используют x как имя переменной. В этом контексте мы знаем, что x представляет age.

Обратите внимание, что метод печати, написанный для объекта PiecePoly, предназначен для печати форматированных полиномиальных уравнений. Чтобы узнать, что такое spl_polulation, используйте:

str(spl_population)
#List of 2
# $ PiecePoly:List of 2
#  ..$ coef : num [1:4, 1:4] 5.68e-13 9.42e+01 1.78e-12 1.87 6.46e+02 ...
#  ..$ shift: logi TRUE
# $ knots    : num [1:5] 7.99 12.73 15.76 22.64 39.3
# - attr(*, "class")= chr [1:2] "PiecePoly" "RegSpline"

Фактически есть промежуточный шаг от этой необработанной формы до print, то есть summary. Функция summary создает форматированные строки с использованием необработанных коэффициентов. Проверим это:

summary(spl_population)
#4 piecewise polynomials of degree 3 are constructed!

Вы не видели уравнения, потому что по умолчанию они скрыты (потому что я не хочу распечатывать все части, чтобы очистить экран). Ты можешь попробовать:

str(summary(spl_population))
#4 piecewise polynomials of degree 3 are constructed!
#List of 4
# $ : chr "5.68e-13 + 94.2 * (x - 7.99) + 1.78e-12 * (x - 7.99) ^ 2 + 1.87 * (x - 7.99) ^ 3"
# $ : chr "646 + 220 * (x - 12.7) + 26.6 * (x - 12.7) ^ 2 - 6.24 * (x - 12.7) ^ 3"
# $ : chr "1.38e+03 + 210 * (x - 15.8) - 30.1 * (x - 15.8) ^ 2 + 1.43 * (x - 15.8) ^ 3"
# $ : chr "1.87e+03 - 0.322 * (x - 22.6) - 0.48 * (x - 22.6) ^ 2 + 0.0096 * (x - 22.6) ^ 3"

Если пользователи задаются вопросом, почему эти уравнения являются полиномами от (x - some_number), а не от x, обратитесь к ?PiecePoly, чтобы узнать о разнице между сдвинутой и несмещенной формами.

Построить этот сплайн довольно просто:

plot(spl_population)

Рисунок 1

Если вы обнаружите, что график визуально недостаточно гладкий, увеличьте параметр spread (значение по умолчанию - 3; увеличивайте его медленно!):

plot(spl_population, spread = 5)

Функция построения графика также может рисовать 1-ю (или любую) производную сплайна, используя аргумент deriv. Пытаться

plot(spl_population, spread = 5, deriv = 1)

Рисунок 2

Чтобы найти все локальные экстремумы сплайна, используйте solve + predict:

( xe <- solve(spl_population, b = 0, deriv = 1) )
#[1] 22.45925
( ye <- predict(spl_population, xe) )
#[1] 1871.8
plot(spl_population)
points(xe, ye)

Рисунок 3

Чтобы найти все локальные экстремумы 1-й производной сплайна (это то, что вы ищете), нам нужно взять дополнительную производную:

( xa <- solve(spl_population, b = 0, deriv = 2) )
#[1] 14.15315
( ya <- predict(spl_population, xa, deriv = 1) )
#[1] 258.2323
plot(spl_population, deriv = 1)
points(xa, ya)

Рисунок 4

Результат, который вы получили от sitar::getPeakTrough в своем вопросе, близок, но неверен. Потому что, используя smooth.spline, вы делаете некоторое приближение, тогда как результат, полученный с помощью SplinesUtils, является точным.

Используйте SplinesUtils для обработки отношения, зависящего от образца

Добавляя случайную линейную функцию age для конкретной выборки к сплайну генеральной совокупности, мы получаем сплайн для конкретной выборки.

Случайная линия может быть получена как подобранные коэффициенты случайного эффекта:

( random_line <- random.effects(nspline_model) )
#     (Intercept)         age
#1001  105.366791 -17.3226305
#1002    2.118663  -3.1252693
#1003 -124.348097  25.3333794
#1004   35.118293  -9.2338755
#...          ...         ...

На этом этапе удобнее экспортировать сплайн населенности в несмещенном виде:

spl_population_unshifted <- RegSplineAsPiecePoly(nspline_model, "ns(age, df = 4)", FALSE)
spl_population_unshifted
#4 piecewise polynomials of degree 3 are constructed!
#Use 'summary' to export all of them.
#The first 4 are printed below.
#-1.71e+03 + 453 * x - 44.9 * x ^ 2 + 1.87 * x ^ 3
#1.5e+04 - 3.49e+03 * x + 265 * x ^ 2 - 6.24 * x ^ 3
#-1.5e+04 + 2.22e+03 * x - 97.8 * x ^ 2 + 1.43 * x ^ 3
#1.52e+03 + 36.2 * x - 1.13 * x ^ 2 + 0.0096 * x ^ 3

Теперь, например, чтобы получить сплайн для конкретного образца для id = 1001, мы можем обновить коэффициенты полинома:

( coef_population <- spl_population_unshifted$PiecePoly$coef )
#            [,1]         [,2]          [,3]          [,4]
#[1,] -1708.58691 15031.740239 -14997.457442  1.521760e+03
#[2,]   452.99242 -3491.784683   2224.770784  3.615668e+01
#[3,]   -44.89601   264.959868    -97.787157 -1.131417e+00
#[4,]     1.87259    -6.240303      1.432464  9.595289e-03

coef_1001 <- coef_population  ## initialize
coef_1001[1, ] <- coef_1001[1, ] + random_line[1, 1]  ## update intercept
coef_1001[2, ] <- coef_1001[2, ] + random_line[1, 2]  ## update slope

Сама по себе матрица коэффициентов не является объектом PiecePoly, поэтому мы не можем ни plot, ни solve. Но мы можем взломать:

spl_1001_unshifted <- spl_population_unshifted
spl_1001_unshifted$PiecePoly$coef <- coef_1001

Теперь мы можем построить его, найти и наложить его экстремумы:

plot(spl_1001_unshifted)
( xe_1001 <- solve(spl_1001_unshifted, b = 0, deriv = 1) )
#[1] 20.72561
if (length(xe_1001)) {
  ( ye_1001 <- predict(spl_1001_unshifted, xe_1001) )
  #[1] 1606.861
  points(xe_1001, ye_1001)
}

(Примечание: обертка if обычно необходима, потому что xe_1001 может быть numeric(0), что означает, что сплайн монотонный и нет экстремума!)

Рисунок 5

Мы также можем построить его 1-ю производную, а затем наложить ее экстремумы (это то, что вы ищете):

plot(spl_1001_unshifted, deriv = 1)
( xa_1001 <- solve(spl_1001_unshifted, b = 0, deriv = 2) )
#[1] 7.991786 14.153151 39.304586
if (length(xa_1001)) {
  ( ya_1001 <- predict(spl_1001_unshifted, xa_1001, deriv = 1) )
  #[1] 76.87045 240.90965 -25.63581
  points(xa_1001, ya_1001)
}

Рисунок 6

Этот график показывает, что две границы также отмечены. Это связано с тем, что естественный кубический сплайн должен иметь нулевую 1-ю производную на двух концах. Так что неудивительно. Конечно, компьютерная программа из-за небольшой ошибки округления может включать или не включать эти два конца. Например, они не попадают в случай популяции (см. Последний рисунок в предыдущем разделе). Но если они все же появятся, просто исключите их и оставьте средние.

В заключение, как только вы научитесь делать это для одного образца, вы можете написать цикл для обработки всех образцов.

В ответ на ваш комментарий

Единственная проблема заключается в том, как получить стандартные ошибки / доверительные интервалы для среднего возраста при максимальной скорости (т. Е. Оценить значения для каждого человека).

Это еще один вопрос, который стоит изучить, ответ на который требует таких же (а может и большего) усилия, как и этот. Политика переполнения стека может ожидать, что вы зададите для этого новый вопрос, иначе наше обсуждение здесь никогда не закончится.

person Zheyuan Li    schedule 16.12.2020
comment
Большое спасибо, глядя на эскиз 1-го производного сплайна, созданного с помощью plot (spl_population, spread = 5, diver = 1) - это, похоже, кривая скорости, а ее пики - это возраст при максимальной скорости. Единственная проблема заключается в том, как получить стандартные ошибки / доверительные интервалы для среднего возраста при максимальной скорости (т.е. оценить значения для каждого человека) ... Я считаю, что в этой статье обсуждается их получение onlinelibrary.wiley.com/doi/full/10.1002/sim.7694 - person aelhak; 17.12.2020

arrow_upward
0
arrow_downward

Один из способов решить эту проблему - вычислить наклоны между каждой последовательной точкой для каждого человека. Затем вы можете взять возраст верхней точки для периода с максимальным наклоном. В приведенном ниже коде это достигается. Обратите внимание, что вы можете заменить ns_pred на ls_pred на этапе изменения для переключения моделей.

data2 <- data %>%
  group_by(id) %>%
  filter(n() > 2) %>%
  arrange(age) %>%
  mutate(slope = (ns_pred - lag(ns_pred)) / (age - lag(age))) %>%
  arrange(slope) %>%
  top_n(n = 1)

Похоже, у вас есть несколько редких записей, поэтому я также отфильтровал людей с двумя или меньшим количеством наблюдений. Кроме того, поскольку кажется, что есть люди с большими промежутками между наблюдениями, возможно, было бы более подходящим усреднить два возраста, которые соответствуют максимальному периоду наклона. Однако это скорее дискреционные вопросы исследователя.

person mikeHoncho    schedule 15.12.2020
comment
Спасибо за предложение, однако я хотел бы оценить их, решив производные сплайна - я отредактировал свой вопрос - person aelhak; 15.12.2020