Я пытаюсь подогнать под очень простую нелинейную смешанную модель (кривая Гомперца) без какой-либо иерархической структуры (просто повторяющиеся измерения). Сначала я хочу попробовать только с фиксированными, а затем со случайными эффектами. Это данные (я взял подмножество из 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 R
library(drc)
sv <- drm(NDVIr ~ CumGDD, data = data, fct = G.3())
sv
Я пробовал с пакетом nlme R
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 ((value + t (value)) / 2): ведущий минор порядка 1 не является положительно определенным
traceback()
Я тоже пробовал с пакетом lme4 R
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 находится в кадре данных.
Любое понимание более чем приветствуется. Заранее спасибо!
