2017-08-06 37 views
0

我想沒有任何層次結構(只是重複測量),以適應一個非常簡單的非線性混合模型(姜氏曲線)。首先,我只想在固定和隨機效果下嘗試。 這是數據(我帶的10個樣本的子集,但我有396)非線性混合模型,而不羣結構規範

data <- structure(list(CumGDD = c(124.66, 124.66, 124.66, 124.66, 124.66, 
124.66, 124.66, 124.66, 124.66, 124.66, 198.127916666667, 198.127916666667, 
198.127916666667, 198.127916666667, 198.127916666667, 198.127916666667, 
198.127916666667, 198.127916666667, 198.127916666667, 198.127916666667, 
242.772916666667, 242.772916666667, 242.772916666667, 242.772916666667, 
242.772916666667, 242.772916666667, 242.772916666667, 242.772916666667, 
242.772916666667, 242.772916666667, 275.351666666667, 275.351666666667, 
275.351666666667, 275.351666666667, 275.351666666667, 275.351666666667, 
275.351666666667, 275.351666666667, 275.351666666667, 275.351666666667, 
304.095, 304.095, 304.095, 304.095, 304.095, 304.095, 304.095, 
304.095, 304.095, 304.095, 378.717916666667, 378.717916666667, 
378.717916666667, 378.717916666667, 378.717916666667, 378.717916666667, 
378.717916666667, 378.717916666667, 378.717916666667, 378.717916666667, 
443.024166666667, 443.024166666667, 443.024166666667, 443.024166666667, 
443.024166666667, 443.024166666667, 443.024166666667, 443.024166666667, 
443.024166666667, 443.024166666667, 483.205833333333, 483.205833333333, 
483.205833333333, 483.205833333333, 483.205833333333, 483.205833333333, 
483.205833333333, 483.205833333333, 483.205833333333, 483.205833333333, 
512.984583333333, 512.984583333333, 512.984583333333, 512.984583333333, 
512.984583333333, 512.984583333333, 512.984583333333, 512.984583333333, 
512.984583333333, 512.984583333333, 583.738333333333, 583.738333333333, 
583.738333333333, 583.738333333333, 583.738333333333, 583.738333333333, 
583.738333333333, 583.738333333333, 583.738333333333, 583.738333333333, 
658.390833333333, 658.390833333333, 658.390833333333, 658.390833333333, 
658.390833333333, 658.390833333333, 658.390833333333, 658.390833333333, 
658.390833333333, 658.390833333333, 685.78625, 685.78625, 685.78625, 
685.78625, 685.78625, 685.78625, 685.78625, 685.78625, 685.78625, 
685.78625, 713.718333333333, 713.718333333333, 713.718333333333, 
713.718333333333, 713.718333333333, 713.718333333333, 713.718333333333, 
713.718333333333, 713.718333333333, 713.718333333333, 835.992083333333, 
835.992083333333, 835.992083333333, 835.992083333333, 835.992083333333, 
835.992083333333, 835.992083333333, 835.992083333333, 870.053333333333, 
870.053333333333, 870.053333333333, 870.053333333333, 870.053333333333, 
870.053333333333, 870.053333333333, 870.053333333333, 870.053333333333, 
870.053333333333, 909.068333333333, 909.068333333333, 909.068333333333, 
909.068333333333, 909.068333333333, 909.068333333333, 909.068333333333, 
909.068333333333, 909.068333333333, 997.753333333333, 997.753333333333, 
997.753333333333, 997.753333333333, 997.753333333333, 997.753333333333, 
997.753333333333, 997.753333333333, 997.753333333333, 1061.60916666667, 
1061.60916666667, 1061.60916666667, 1061.60916666667, 1061.60916666667, 
1061.60916666667, 1061.60916666667, 1109.1775, 1109.1775, 1109.1775, 
1109.1775, 1109.1775, 1109.1775, 1109.1775, 1109.1775, 1164.07458333333, 
1164.07458333333, 1164.07458333333, 1164.07458333333, 1164.07458333333, 
1164.07458333333, 1164.07458333333, 1164.07458333333, 1223.085, 
1223.085, 1223.085, 1223.085, 1223.085, 1223.085, 1223.085, 1223.085, 
1223.085, 1223.085, 1268.74041666667, 1268.74041666667, 1268.74041666667, 
1268.74041666667, 1268.74041666667, 1268.74041666667, 1268.74041666667, 
1268.74041666667, 1268.74041666667, 1268.74041666667, 1308.92041666667, 
1308.92041666667, 1308.92041666667, 1308.92041666667, 1308.92041666667, 
1308.92041666667, 1308.92041666667, 1308.92041666667, 1353.40125, 
1353.40125, 1353.40125, 1353.40125, 1353.40125, 1353.40125, 1353.40125, 
1353.40125, 1353.40125, 1353.40125, 1386.49166666667, 1386.49166666667, 
1386.49166666667, 1386.49166666667, 1386.49166666667, 1386.49166666667, 
1386.49166666667, 1386.49166666667, 1386.49166666667, 1386.49166666667, 
1420.16833333333, 1420.16833333333, 1420.16833333333, 1420.16833333333, 
1420.16833333333, 1420.16833333333, 1420.16833333333, 1420.16833333333, 
1420.16833333333, 1420.16833333333, 1465.77083333333, 1465.77083333333, 
1465.77083333333, 1465.77083333333, 1465.77083333333, 1465.77083333333, 
1465.77083333333, 1509.21416666667, 1509.21416666667, 1509.21416666667, 
1509.21416666667, 1509.21416666667, 1509.21416666667, 1509.21416666667, 
1509.21416666667), NDVIr = c(0.326734537556686, 0.34635103511923, 
0.41413832498987, 0.351310749147721, 0.393142917023234, 0.346386159863839, 
0.470936633563708, 0.393683553738797, 0.477145162676157, 0.374332602869744, 
0.399656917554062, 0.405729803302398, 0.462392583145425, 0.385385634302403, 
0.284644802690205, 0.341708626865214, 0.373245280739098, 0.297082860395878, 
0.338992263970809, 0.328009446471202, 0.3864924745074, 0.338233219314218, 
0.414182270012836, 0.381044071285413, 0.27853092703881, 0.334093310912654, 
0.396283587652542, 0.333797210629104, 0.33288365893073, 0.338806141443146, 
0.396692170766243, 0.402010070712185, 0.477692555140439, 0.583054609863216, 
0.33248604169778, 0.388480831363342, 0.394041596269905, 0.396121088530724, 
0.445932737630167, 0.406255657534553, 0.480521164786982, 0.438116936717331, 
0.500943099169269, 0.574114274948364, 0.468465182355142, 0.438476710317985, 
0.508222915780569, 0.471488038890975, 0.494881991280694, 0.487205003301489, 
0.582759919048259, 0.580840031834031, 0.636409022987258, 0.668700748285102, 
0.569981692001408, 0.527429372767094, 0.585363727232429, 0.595424001202526, 
0.605682862491717, 0.576679575419451, 0.620146320429479, 0.614639916488238, 
0.694476703246357, 0.613780016162866, 0.638307007861515, 0.629767618439801, 
0.67601226305135, 0.647827001051488, 0.60251563176366, 0.652981338154494, 
0.545897498488816, 0.576573102801764, 0.651135693716886, 0.613480310032889, 
0.610301226779858, 0.656782465200423, 0.622287535397165, 0.632757970716384, 
0.606639447102299, 0.618769753052964, 0.587778836508016, 0.583678448342968, 
0.661228984828489, 0.592678976584954, 0.584443712907538, 0.640458831034116, 
0.634544825544102, 0.582358840543761, 0.62446711093034, 0.645849315193373, 
0.644237881152471, 0.670350391520271, 0.691637067201447, 0.695190444098271, 
0.691394387522623, 0.737855968215111, 0.697492990890306, 0.707429622258371, 
0.686456658647713, 0.738873457213227, 0.617541840377972, 0.706422161253019, 
0.71025496607701, 0.712495891471769, 0.709308696776658, 0.72408051102251, 
0.711917453325673, 0.673964854327079, 0.674452192699244, 0.692371436599349, 
0.533800749647061, 0.66819940121409, 0.659604328631773, 0.64363430526236, 
0.667660060188645, 0.670391061160267, 0.633111535203985, 0.639627330791422, 
0.658698545446175, 0.677339420961202, 0.648138134767384, 0.640074788026247, 
0.693005678880091, 0.657859308607955, 0.690383624907919, 0.71885911210521, 
0.708693883372722, 0.688129672503304, 0.694178199870893, 0.683109547482279, 
0.65096223415747, 0.645964783381146, 0.678481468832164, 0.717443844344107, 
0.741072132668907, 0.757662926321472, 0.742872179704767, 0.708541794658953, 
0.664657940009943, 0.700151429494567, 0.683095773693102, 0.675361184526528, 
0.71700081070095, 0.699138073359393, 0.72230869829755, 0.67662429008069, 
0.673751607473228, 0.725647610769991, 0.667739928913317, 0.708865651147583, 
0.67287378458334, 0.682184345860339, 0.740350385885337, 0.727258756451038, 
0.728589635252484, 0.689473748097639, 0.70551431523546, 0.657733014543568, 
0.673279525146335, 0.655096309990237, 0.732820064995832, 0.756668414106169, 
0.740753806231525, 0.725855875527774, 0.697966218530406, 0.726641271757911, 
0.703228834250003, 0.696255542746872, 0.758071641902502, 0.743409728013122, 
0.76129008710658, 0.768630092782366, 0.748669172624243, 0.674195927717801, 
0.752013797871578, 0.714829530594045, 0.747646135447338, 0.757461729402785, 
0.729581787835431, 0.728841943820421, 0.769516084221819, 0.703774076752448, 
0.702515263428803, 0.718605605350307, 0.766389826687477, 0.680061945544163, 
0.724215955387007, 0.732058965732438, 0.735653731478568, 0.652530832382112, 
0.72686370110552, 0.680375685387295, 0.736603377306099, 0.712724798670091, 
0.746652726085055, 0.725535181246623, 0.703603402186289, 0.723437714770724, 
0.721356258293984, 0.675886723363852, 0.688142017694115, 0.743067388281102, 
0.758275152601243, 0.730431173674391, 0.74650872966109, 0.7378875471765, 
0.737214160722912, 0.754093200845553, 0.74987379796935, 0.673571663557409, 
0.747632763176418, 0.701454462346379, 0.710761528481747, 0.701370454047384, 
0.695399954171373, 0.703877556668931, 0.714144530907548, 0.689335025711769, 
0.729373446048829, 0.71031954060566, 0.73159841576388, 0.717918969997454, 
0.73020737891142, 0.716490651754628, 0.673751702743956, 0.69102014860641, 
0.728911443402563, 0.700885207746844, 0.71642901303611, 0.736495017307433, 
0.717658913350811, 0.688042925533586, 0.720292021535713, 0.712222809826439, 
0.68173364585245, 0.705327637704497, 0.617025646437043, 0.741437556633414, 
0.725191637769468, 0.774702607053949, 0.762398649890891, 0.775791416700744, 
0.782116753721391, 0.786733589301968, 0.745576974597893, 0.751720639746142, 
0.743906463887347, 0.695837087044011, 0.731147780678493, 0.739828109231531, 
0.668630156814454, 0.68138726516326, 0.737585692878805, 0.719679687801702, 
0.666348512128941, 0.750737638141788, 0.708418042251955, 0.777239294856883, 
0.732865875474648, 0.724267563473574, 0.747340495998486, 0.735491104316784 
), ID = structure(c(1L, 2L, 3L, 4L, 5L, 6L, 7L, 8L, 9L, 10L, 
1L, 2L, 3L, 4L, 5L, 6L, 7L, 8L, 9L, 10L, 1L, 2L, 3L, 4L, 5L, 
6L, 7L, 8L, 9L, 10L, 1L, 2L, 3L, 4L, 5L, 6L, 7L, 8L, 9L, 10L, 
1L, 2L, 3L, 4L, 5L, 6L, 7L, 8L, 9L, 10L, 1L, 2L, 3L, 4L, 5L, 
6L, 7L, 8L, 9L, 10L, 1L, 2L, 3L, 4L, 5L, 6L, 7L, 8L, 9L, 10L, 
1L, 2L, 3L, 4L, 5L, 6L, 7L, 8L, 9L, 10L, 1L, 2L, 3L, 4L, 5L, 
6L, 7L, 8L, 9L, 10L, 1L, 2L, 3L, 4L, 5L, 6L, 7L, 8L, 9L, 10L, 
1L, 2L, 3L, 4L, 5L, 6L, 7L, 8L, 9L, 10L, 1L, 2L, 3L, 4L, 5L, 
6L, 7L, 8L, 9L, 10L, 1L, 2L, 3L, 4L, 5L, 6L, 7L, 8L, 9L, 10L, 
1L, 2L, 3L, 4L, 5L, 6L, 9L, 10L, 1L, 2L, 3L, 4L, 5L, 6L, 7L, 
8L, 9L, 10L, 1L, 2L, 3L, 4L, 5L, 6L, 7L, 8L, 10L, 1L, 2L, 3L, 
4L, 5L, 6L, 7L, 8L, 10L, 2L, 3L, 4L, 5L, 6L, 8L, 9L, 1L, 2L, 
3L, 5L, 6L, 7L, 8L, 10L, 1L, 2L, 3L, 4L, 6L, 8L, 9L, 10L, 1L, 
2L, 3L, 4L, 5L, 6L, 7L, 8L, 9L, 10L, 1L, 2L, 3L, 4L, 5L, 6L, 
7L, 8L, 9L, 10L, 1L, 2L, 3L, 5L, 6L, 7L, 8L, 9L, 1L, 2L, 3L, 
4L, 5L, 6L, 7L, 8L, 9L, 10L, 1L, 2L, 3L, 4L, 5L, 6L, 7L, 8L, 
9L, 10L, 1L, 2L, 3L, 4L, 5L, 6L, 7L, 8L, 9L, 10L, 1L, 2L, 3L, 
7L, 8L, 9L, 10L, 1L, 2L, 3L, 4L, 5L, 7L, 9L, 10L), .Label = c("1", 
"2", "3", "4", "5", "6", "7", "8", "9", "10", "11", "12", "13", 
"14", "16", "17", "19", "20", "21", "22", "23", "24", "25", "26", 
"27", "28", "29", "30", "31", "32", "33", "34", "35", "36", "37", 
"38", "39", "40", "41", "42", "43", "44", "45", "46", "48", "49", 
"50", "51", "53", "54", "55", "56", "57", "58", "59", "60", "61", 
"62", "63", "64", "65", "66", "67", "68", "69", "70", "71", "72", 
"73", "74", "75", "76", "77", "78", "79", "80", "81", "82", "83", 
"84", "85", "87", "88", "89", "90", "91", "92", "93", "94", "95", 
"96", "97", "98", "99", "100", "101", "102", "103", "104", "105", 
"106", "107", "108", "109", "110", "111", "112", "113", "115", 
"116", "117", "118", "119", "120", "121", "122", "123", "124", 
"125", "126", "127", "128", "129", "130", "131", "132", "133", 
"134", "135", "136", "137", "138", "139", "140", "141", "142", 
"143", "144", "145", "146", "148", "149", "150", "152", "153", 
"154", "155", "156", "157", "158", "159", "160", "161", "162", 
"163", "164", "165", "166", "167", "168", "169", "171", "172", 
"173", "174", "175", "176", "177", "178", "179", "180", "181", 
"182", "183", "184", "186", "187", "188", "189", "190", "191", 
"192", "193", "194", "195", "196", "197", "198", "199", "200", 
"201", "202", "203", "204", "205", "206", "207", "208", "209", 
"210", "211", "212", "213", "214", "215", "216", "218", "219", 
"220", "221", "222", "223", "224", "225", "226", "227", "228", 
"229", "230", "231", "232", "233", "234", "235", "236", "237", 
"238", "239", "240", "241", "242", "243", "246", "247", "248", 
"251", "252", "253", "254", "255", "256", "257", "258", "259", 
"260", "261", "262", "263", "264", "265", "266", "267", "268", 
"269", "270", "271", "272", "273", "274", "275", "276", "277", 
"278", "280", "281", "283", "284", "285", "286", "287", "288", 
"289", "290", "291", "292", "293", "294", "295", "296", "297", 
"298", "299", "300", "301", "302", "303", "304", "305", "306", 
"307", "308", "309", "310", "311", "312", "313", "314", "315", 
"316", "317", "318", "319", "320", "321", "322", "323", "324", 
"325", "326", "327", "328", "329", "330", "331", "332", "333", 
"334", "335", "336", "337", "338", "339", "340", "341", "342", 
"343", "344", "345", "346", "347", "348", "349", "350", "351", 
"352", "353", "354", "355", "356", "357", "358", "359", "360", 
"361", "362", "363", "364", "365", "366", "367", "368", "369", 
"370", "371", "372", "373", "374", "375", "376", "377", "379", 
"381", "382", "383", "384", "385", "386", "387", "388", "389", 
"390", "391", "392", "393", "394", "395", "396"), class = "factor")), .Names = c("CumGDD", 
"NDVIr", "ID"), row.names = c(NA, -262L), class = "data.frame") 

初始值與DRCř包提取

library(drc) 
sv <- drm(NDVIr ~ CumGDD, data = data, fct = G.3()) 
sv 

我試圖與NLMEř包

library(nlme) 
model.nlme <- nlme(NDVIr ~ d*exp(-exp(b*(CumGDD-e))), 
         data = data, 
         fixed = d + b + e ~ 1, 
         #random = d + b + e ~ 1|ID, 
         groups = ~ ID, 
         start = c(b=-0.004, d=0.728, e=91.236)) 

但是我得到這個錯誤:錯誤chol.default((值+ T(值))/ 2): 1階主導次要不爲正明確

traceback() 

我也與lme4ř包

library(lme4) 
gomp <- ~ d*exp(-exp(b*(CumGDD-e))) 
gomp.deriv <- deriv(gomp,namevec=c("d","b","e"), function.arg=c("CumGDD","d","b","e")) 
model.nlmer <- nlmer(NDVIr ~ gomp.deriv(CumGDD,d,b,e), 
         #(d|ID) + (b|ID) + (e|ID), 
         data= data, 
         start = c(b=-0.004, d=0.728, e=91.236)) 

在這種情況下我得到一個完全不同的錯誤嘗試:錯誤的eval(形式[[2]]):找不到對象'NDVIr'。

traceback() 

NLME [R包的錯誤似乎涉及到的參數的計算,但我不知道這是否是模型的錯誤規範。我期待得到與nls函數類似的結果(只考慮固定效應)。此外,我不知道組的規格(ID變量)是否正確,因爲我不想要任何分組結構(我假設所有樣本都來自同一羣體),但是該軟件包要求我進行分組結構體。另一方面,我不知道爲什麼lme4 R軟件包的錯誤正在發生。變量NDVIr在數據框中。

任何有識之士的歡迎。提前致謝!

回答

2

這是一個帶自啓動功能,SSgompertz一個解決方案,但它的規格是從姜氏模式不同。在任何情況下,問題可能不是與起始值,所以你可以修改這個爲你工作。 編輯:這也與您的重新參數化版本的作品,請看下文。

# Find starting values with nls to the entire dataset 
fit_nlme_0 <- nls(NDVIr ~ SSgompertz(CumGDD, Asym, b2, b3), data=data) 
p <- coef(fit_nlme_0) 

# Fit non-linear mixed-effects model 
library(nlme) 
fit_nlme <- nlme(NDVIr ~ SSgompertz(CumGDD, Asym, b2, b3) , 
       data = data, 
       fixed = list(Asym + b2 + b3 ~ 1), 
       random = Asym + b2 + b3 ~ 1 | ID, 
      start=list(fixed=c(Asym=p["Asym"], b2=p["b2"], b3=p["b3"]))) 

# Coefficients: 
coef(fit_nlme) 

     Asym  b2  b3 
1 0.7031033 1.342452 0.9959273 
2 0.7253884 1.532975 0.9956931 
3 0.7268661 1.133122 0.9959026 
4 0.7382484 1.319501 0.9957342 
5 0.7378579 1.774258 0.9954877 
6 0.7440437 1.739252 0.9954699 
7 0.7316711 1.312134 0.9957769 
8 0.7264944 1.536741 0.9956833 
9 0.7311114 1.391140 0.9957375 
10 0.7366816 1.578322 0.9956006 

# And finally a simple plot of the fits to judge 
# both the quality of the fit and the variance due to 
# the random effect 
p <- coef(fit_nlme) 
palette(rainbow(10)) 
with(data, plot(CumGDD, NDVIr, pch=16, col=as.factor(ID), ylim=c(0,1))) 
for(i in 1:10){ 
    curve(p[i,"Asym"]*exp(-p[i,"b2"]*p[i,"b3"]^x), col=palette()[i], add=T) 
} 

enter image description here

編輯:您重新參數化的版本,我們也可以得到這個適合:

fit_nlme_0 <- nls(NDVIr ~ d*exp(-exp(b*(CumGDD-e))), data=data, 
       start=list(b=-0.004, d=0.728, e=91.236)) 
p <- coef(fit_nlme_0) 

fit_nlme <- nlme(NDVIr ~ d*exp(-exp(b*(CumGDD-e))), 
       data = data, 
       fixed = list(d + b + e ~ 1), 
       random = d + b + e ~ 1 | ID, 
       start=list(fixed=c(d=p["d"], b=p["b"], e=p["e"]))) 

但要注意的是,直接使用圓潤的初始值,它不收斂(大概不夠準確)。這就是爲什麼我總是喜歡首先使用nls來獲得起始值,它似乎是(通常)一個相當可靠的方法。

+0

感謝您的解決方案。然而,我仍然想知道爲什麼我的Gompertz的reparametrization不起作用。 – Diego

+0

查看我的編輯;也適用於您的示例 - 但需要重新設置起始值。 –

+0

我很欣賞你的幫助。從nls直接提取係數似乎是最好的選擇,而不是四捨五入。 – Diego