nonlinear mixed model without group structure specification

0
votes

I am trying to fit a very simple nonlinear mixed model (Gompertz curve) without any hierarchical structure (just repeated measures). First, I want to try it only with fixed and then with random effects. This is the data (I took a subset of 10 samples but I have 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")

The initial values were extracted with drc R package

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

I tried with the nlme R package

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

However I got this error: Error in chol.default((value + t(value))/2) : the leading minor of order 1 is not positive definite

traceback()

I was also trying with the lme4 R package

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

And in this case I am getting a totally different error: Error in eval(form[[2]]) : object 'NDVIr' not found.

traceback()

The error of nlme R package seems to be related to the calculation of parameters, however I do not know if it is a wrong specification of the model. I was expecting to get similar results as nls function (just considering the fixed effects). In addition, I do not know if the specification of the groups (ID variable) is correct since I do not want any grouping structure (I am assuming all the samples are coming from the same population), however the package asked me for a grouping structure. On the other hand, I have no idea why the error of lme4 R package is happening. The variable NDVIr is in the dataframe.

Any insight is more than welcome. Thanks in advance!

1
rlme4non-linear-regressionnlme

1 Answers

3
votes

Here is a solution with a self-starting function, SSgompertz, though its specification is different from your Gompertz model. In any case the issue is probably not with the starting values so you can probably modify this to work for you. Edit: this also works with your reparameterized version, see below.

# 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

Edit: with your reparameterized version, we can also get this to fit:

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

Though note that using your rounded starting values directly, it does not converge (presumably not quite accurate enough). This is why I always like to use nls first to get starting values, it seems to be (usually) a fairly robust approach.