如何解决来自线性混合自然样条模型的峰值速度下的个体年龄
此线程以 Get subject-specific peak velocity and age at peak velocity values from linear mixed spline models 继续。
我正在拟合具有年龄的自然样条函数的线性混合效应模型。我想通过区分样条项来估计数据集中每个人的峰值速度（apv - 年）和峰值速度（pv - 克）的年龄。该模型包括年龄的随机二次斜率。
如何估计特定于个人的 apv 和 pv？我正在使用 SplinesUtils 包。
示例数据：
dat <- structure(list(id = c(1001L,1001L,1002L,1003L,1004L,1005L,1006L,1007L,1008L,1009L,1010L,1011L,1012L,1013L,1014L,1015L,1016L,1017L,1018L,1020L,1021L,1022L,1023L,1024L,1025L,1026L,1027L,1028L,1029L,1030L,1031L,1032L,1033L,1034L,1035L,1036L,1037L,1039L,1040L,1041L,1042L,1043L,1044L,1045L,1047L,1048L,1049L,1051L,1052L,1053L,1054L,1056L,1057L,1058L,1060L,1061L,1062L,1063L,1064L,1065L,1066L,1067L,1068L,1069L,1070L,1071L,1072L,1073L,1074L,1075L,1076L,1077L,1078L,1080L,1081L,1082L,1083L,1084L,1085L,1086L,1087L,1088L,1089L,1090L,1091L,1092L,1093L,1094L,1095L,1096L,1097L,1098L,1099L,1100L,1101L,1103L,1104L,1105L,1106L,1107L,1108L,1110L,1111L,1112L,1117L,1123L),y = c(1934.047646,1075.598345,1956.214821,2000.38538,732.315937,3119.86,624.951231,791.2764892,1884.530826,1047.57,1238.306103,1555.042976,2547.870529,2467.385,1181.635212,565.306282,2016.027874,712.6134567,635.2537841,2167.362267,2575.574188,2480.028259,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,2115.374241,2331.46,875.949604,2241.945103,2319.666939,2316.220234,719.7139549,2042.803307,1132.977503,1737.18,1351.629826,1291.44593,1108.26586,1028.979719,2068.934227,2440.784416,1036.72,894.6663704,2449.184731,1109.9,672.9310664,2072.320354,2114.215416,1805.422001,2461.18,2101.374248,2105.879,1600.086481,2866.84,2807.311,3055.569931,2602.287521,2690.007614,620.5975037,2608.4,2722.3,2713.66185,1590.002,2198.211,2488.097725,2322.616348,2627.1,2418.328346,2601.661034,531.7369251,811.9494571,884.31,768.0526981,652.1271248,2767.479,1047.144354,1995.119,707.6093158,1120.650104,3036.591904,3081.86,1193.583691,2056.569244,1823.155,1238.948124,2124.685,887.20438,2560.155342,3095.923164,3003.729011,2861.12,2735.26,822.8209591,1648.951,906.7692623,582.787096,1286.45,797.2365359,2566.770554,2666.41,2045.320816,2401.21,2583.2,2581.32,2622.357,2588.462498,442.433671,1251.627064,406.2565479,2108.787437,983.1101169,2102.085403,1155.713411,1909.797131,2871.55,2711.07,2883.22245,3027.103172,3108.21537,3007.87294,3208.963631,2617.91,2457.464466,2890.51,2698.48214,2700.723,2817.668579,1349.90691,1476.19994,1552.95,925.8325004,1258.28,840.1875095,2405.175911,1056.678543,1571.936,1210.89,673.7005405,687.7842464,1016.86,1217.866,1493.791817,2246.726913,1054.821,563.6580887,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,912.8781966,1257.222706,2586.922356,1608.28,962.5674305,1085.451181,2539.218132,2535.526085,2561.60054,1600.198,2100.048149,758.3851737,2643.373329,367.7795143,866.0683727,718.5049658,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,1895.828438,1757.75,1050.279345,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,1675.768558,3157.550789,1565.173126,2404.836883,2541.045593,585.7223682,2465.177761,2678.462074,500.3733997,781.342,898.3551559,1807.02,1418.888027,1797.36,2200.06,2218.369926,1986.642735,2088.292,2069.139,1507.901432,2061.395798,2075.164864,2081.913219,483.8579493,1857.88,2578.772636,1039.632153,2288.28,1831.349922,2349.23,933.1002788,2626.298935,1521.744,1984.760715,2450.333,1732.339031,2731.9,869.2320918,1785.72,1922.798,3081.28,1508.8,2421.288597,1268.074959,1569.05,1808.115,2165.724808,2084.149837,2693.027184,2464.489,2607.653496,1012.837271,2673.190872,2635.290516,2773.42,2654.772674,2377.905655,2679.014969,1226.40016,1470.69,1273.789799,2294.926086,1873.817,2274.930534,2317.429165,959.1709613,1328.159428,1630.28,1610.54982,2507.05302,750.467966,821.2255058,802.8240452,2829.47879),age = c(31.54004107,11.95071869,27.88501027,25.07871321,10.90759754,25.70020534,9.560574949,11.17864476,15.8384668,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,7.991786448,17.26488706,15.66872005,17.63723477,30.97330595,17.5633128,30.11088296,23.31279945,20.58590007,28.27926078,11.66324435,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,10.05065024,23.6659822,16.22450376,28.74743326,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,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.24024641,20.70362765,13.5687885,17.13347023,28.11498973,24.16974675,18.19575633,27.73442847,15.52361396,11.76728268,10.98699521,11.51540041,9.902806297,13.05407255,8.703627652,25.60164271,10.59000684,14.45859001,14.05886379,10.88295688,10.75427789,26.50513347,18.83093771,22.86379192,11.8384668,15.04449008,15.42505133,14.14099932,28.06844627,14.66119097,13.79055441,15.37850787,22.58179329,30.0752909,21.85900068,15.29089665,26.79534565,11.68514716,15.58384668,15.08555784,14.11909651,10.21765914,12.1670089,10.50239562,23.3045859,15.92607803,16.65982204,32.56947296,16.90349076,25.12799452,17.88364134,19.46338125,8.736481862,17.68104038,14.54893908,12.98562628,22.45311431,38.68856947,25.44010951,28.70910335,19.21697467,29.45106092,33.31690623,16.68172485,15.816564,24.89801506,18.7761807,18.4366872,19.45790554,19.78370979,14.98973306,15.89869952,29.06502396,10.74880219,13.47843943,10.5982204,24.61875428,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,35.04996578,15.61943874,13.33333333,10.56810404,11.34017796,13.5797399,28.79671458,12.56673511,12.55578371,30.80082136,23.63039014,29.66461328,17.46748802,9.768651608,13.46748802,13.24298426,26.87474333,27.43326489,20.6899384,10.0752909,13.37713895,28.38056126,8.911704312,24.62149213,14.32443532,10.24229979,13.87268994,11.44421629,21.68377823,27.97809719,28.90075291,24.64339493,10.61190965,15.8110883,14.25051335,13.64818617,26.05338809,13.69746749,23.98083504,20.42162902,12.68172485,11.51813826,15.49897331,18.70225873,17.47570157,14.66666667,26.83915127,13.29226557,18.14647502,14.67761807,16.61601643,9.812457221,15.96714579,17.61806982,11.87953457,11.80561259,19.15400411,15.70704997,12.35318275,18.12457221,16.8733744,32.02464066,25.30047912,16.13415469,19.37850787,25.42368241,16.05201916,15.43874059,9.158110883,14.39014374,22.12183436,15.35934292,28.45995893,17.06502396,26.32991102,12.38056126,16.42436687,11.70978782,17.62628337,15.11019849,14.09993155,21.89185489,17.73305955,25.55509925,14.75975359,24.03559206,14.36002738,12.73100616,16.09034908,13.69472964,23.03901437,16.94182067,13.99315537,15.65776865,19.25530459,10.43394935,12.72826831,24.25735797,37.41820671,25.25393566,12.11772758,14.19575633,14.091718,15.10746064,13.16906229,12.09856263,36.3504449,22.68035592,11.21149897,13.34702259,14.5982204,11.31827515,15.14579055,15.44969199,12.43531828,12.72005476,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)
nspline_model <- lme(y ~ ns(age,df = 3),data = dat,random = ~ poly(age,2) | id)
这是显示我如何尝试计算两个人的 apv 和 pv 的代码 - 但我得到了 pv 的负值，以及 ~39 年的 apv 值
# 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,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,2] + random_quadratic[2,3]
spl_1002_unshifted <- spl_population_unshifted
spl_1002_unshifted$Piecepoly$coef <- coef_1002
(apv_1002 <- solve(spl_1002_unshifted,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))
解决方法
注意：OP 更新了他的问题，在得知我的初始答案后，将随机线（如 the previous thread 中使用的）替换为随机二次方。然后更新我的答案以适应这一点。
新答案
当使用随机二次式时，用原始多项式而不是（默认）正交多项式指定它，否则所得多项式系数有不同的解释。
nspline_model <- lme(y ~ ns(age,df = 3),data = dat,random = ~ poly(age,2,raw = TRUE) | id)
random_quadratic <- random.effects(nspline_model)
然后在使用 SplinesUtils 时，请确保正确添加系数：
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,2] ## age ^ 1
coef_1001[3,] <- coef_1001[3,3] ## age ^ 2
然后您可以继续：
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))
初步回答
这确实令人费解，但我终于转过头来。我们的代码没有任何问题；我们所看到的在数学上是有保证的。请注意：
subject spline = population spline + linear line
取二阶导数后，我们有（因为直线的二阶导数为0）：
subject spline 2nd derivative = population spline 2nd derivative
因此，所有受试者的峰值速度（二阶导数为 0）的年龄都是相同的！然而，受试者之间的峰值速度是不同的。
如果我们希望在峰值速度下得出不同的年龄，我们需要在模型中使用随机二次曲线而不是随机线。但是建立统计模型不是数学游戏，所以我们需要三思而后行。
版权声明：本文内容由互联网用户自发贡献，该文观点与技术仅代表作者本人。本站仅提供信息存储空间服务，不拥有所有权，不承担相关法律责任。如发现本站有涉嫌侵权/违法违规的内容， 请发送邮件至 dio@foxmail.com 举报，一经查实，本站将立刻删除。