Reputation: 415
This thread continues with Get subject-specific peak velocity and age at peak velocity values from linear mixed spline models.
I am fitting a linear mixed effects model with a natural spline function for age. I would like to estimate age at peak velocity (apv - years) and peak velocity (pv - grams) for each person in the dataset by differentiating the spline terms. The model includes a random quadratic slope for age.
How can I estimate the person-specific apv and pv? I am using the SplinesUtils
package.
Example data:
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")
Code for fitting the linear mixed spline model:
library(nlme)
library(splines)
library(tidyverse)
nspline_model <- lme(y ~ ns(age, df = 3), data = dat, random = ~ poly(age, 2) | id)
This is code showing how I tried to calculate apv and pv for two individuals - but i get a negative value for pv, and an apv value of ~39 years
# LIST OF PERSON IDs
dat %>% distinct(id) %>% pull()
library(SplinesUtils)
(random_quadratic <- random.effects(nspline_model))
spl_population_unshifted <- RegSplineAsPiecePoly(nspline_model, "ns(age, df = 3)", FALSE)
# ID = 1001
(coef_1001 <- spl_population_unshifted$PiecePoly$coef)
coef_1001[1, ] <- coef_1001[1, ] + random_quadratic[1, 1]
coef_1001[2, ] <- coef_1001[2, ] + random_quadratic[1, 2] + random_quadratic[1, 3]
spl_1001_unshifted <- spl_population_unshifted
spl_1001_unshifted$PiecePoly$coef <- coef_1001
(apv_1001 <- solve(spl_1001_unshifted, b = 0, deriv = 2))
(pv_1001 <- predict(spl_1001_unshifted, apv_1001, deriv = 1))
(apv_pv_1001 <- as.data.frame(cbind(apv_1001, pv_1001)))
(apv_pv_1001 <- apv_pv_1001 %>% top_n(1, pv_1001))
# ID = 1002
(coef_1002 <- spl_population_unshifted$PiecePoly$coef)
coef_1002[1, ] <- coef_1002[1, ] + random_quadratic[2, 1]
coef_1002[2, ] <- coef_1002[2, ] + random_quadratic[2, 2] + random_quadratic[2, 3]
spl_1002_unshifted <- spl_population_unshifted
spl_1002_unshifted$PiecePoly$coef <- coef_1002
(apv_1002 <- solve(spl_1002_unshifted, b = 0, deriv = 2))
(pv_1002 <- predict(spl_1002_unshifted, apv_1002, deriv = 1))
(apv_pv_1002 <- as.data.frame(cbind(apv_1002, pv_1002)))
(apv_pv_1002 <- apv_pv_1002 %>% top_n(1, pv_1002))
Upvotes: 0
Views: 104
Reputation: 73385
Note: OP updated his question, replacing the random line (as used in the previous thread) with a random quadratic after being informed by my initial answer. My answer is then renewed to accommodate that.
New answer
When using a random quadratic, specify it with raw polynomial instead of (the default) orthogonal polynomial, otherwise the resulting polynomial coefficients have a different interpretation.
nspline_model <- lme(y ~ ns(age, df = 3), data = dat,
random = ~ poly(age, 2, raw = TRUE) | id)
random_quadratic <- random.effects(nspline_model)
Then when using SplinesUtils, make sure you add coefficients correctly:
library(SplinesUtils)
spl_population_unshifted <- RegSplineAsPiecePoly(nspline_model, "ns(age, df = 3)", FALSE)
# ID = 1001
(coef_1001 <- spl_population_unshifted$PiecePoly$coef)
coef_1001[1, ] <- coef_1001[1, ] + random_quadratic[1, 1] ## age ^ 0
coef_1001[2, ] <- coef_1001[2, ] + random_quadratic[1, 2] ## age ^ 1
coef_1001[3, ] <- coef_1001[3, ] + random_quadratic[1, 3] ## age ^ 2
Then you can proceed:
spl_1001_unshifted <- spl_population_unshifted
spl_1001_unshifted$PiecePoly$coef <- coef_1001
(apv_1001 <- solve(spl_1001_unshifted, b = 0, deriv = 2))
(pv_1001 <- predict(spl_1001_unshifted, apv_1001, deriv = 1))
(apv_pv_1001 <- as.data.frame(cbind(apv_1001, pv_1001)))
(apv_pv_1001 <- apv_pv_1001 %>% top_n(1, pv_1001))
Initial answer
This was indeed puzzling but I finally had my head turned around. There is nothing wrong with our code; what we see is mathematically guaranteed. Note that:
subject spline = population spline + linear line
After taking 2nd derivative, we have (since the 2nd derivative of a linear line is 0):
subject spline 2nd derivative = population spline 2nd derivative
Therefore, the age at peak velocity (where the 2nd derivative hits 0) is identical for all subjects! However, the peak velocity is different between subjects.
If we expect to derive a different age at peak velocity we need a random quadratic rather than a random line in the model. But building statistical model is not a mathematical game, so we need think twice about it.
Upvotes: 3