modèle mixte non linéaire sans spécification de structure du groupe

Question

Je suis en train d'installer un modèle mixte non linéaire très simple (courbe Gompertz) sans aucune structure hiérarchique (seulement des mesures répétées). D'abord, je veux l'essayer seulement avec des effets fixes puis avec des effets aléatoires. Ce sont les données (I, a eu un sous-ensemble de 10 échantillons, mais je 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") 

Les valeurs initiales ont été extraites avec le paquet DRC R

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

I essayé avec le paquet de R 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)) 

Cependant, je suis arrivé cette erreur: erreur dans chol.default ((valeur + t (valeur))/2): le leader mineur de l'ordre 1 n'est pas positif défini

traceback() 

J'ai aussi essayé avec le paquet R 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)) 

Et dans ce cas, je reçois une erreur tout à fait différente: Erreur eval (forme [[2]]): objet 'NDVIr' non trouvé.

traceback() 

L'erreur du paquet nlme R semble être lié au calcul des paramètres, mais je ne sais pas si elle est une mauvaise spécification du modèle. Je m'attendais à obtenir des résultats similaires à la fonction nls (juste en tenant compte des effets fixes). De plus, je ne sais pas si la spécification des groupes (variable ID) est correcte car je ne veux pas de structure de regroupement (je suppose que tous les échantillons proviennent de la même population), cependant le package m'a demandé un groupement structure. D'un autre côté, je ne sais pas pourquoi l'erreur du paquet lme4 R est en train de se produire. La variable NDVIr est dans la trame de données.

Toute idée est plus que bienvenue. Merci d'avance!

Answer 1

Voici une solution avec une fonction d'auto-démarrage, SSgompertz, bien que sa spécification est différent de votre modèle de Gompertz. Dans tous les cas, le problème n'est probablement pas les valeurs de départ, donc vous pouvez probablement modifier cela pour qu'il fonctionne pour vous. Editer: cela fonctionne également avec votre version reparameterized, voir ci-dessous.

# 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) 
} 

Edit: avec votre version reparamétrés, nous pouvons également obtenir ceci pour adapter:

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"]))) 

Bien noter que l'utilisation de vos valeurs initiales arrondies directement, il ne converge pas (probablement pas assez précis). C'est pourquoi j'aime toujours utiliser nls d'abord pour obtenir des valeurs de départ, il semble être (généralement) une approche assez robuste.