Discussion :

[robert duparvis]

Bonjour,

J'ai un petit problème avec la fonction nls de R pour le modèle suivant : a-b*exp(-c/t^d), les coefficients recherchés sont a, b, c et d.
J'obtiens en générale le message d'erreur suivant "gradient singulier" ou "valeur manquante ou infinie obtenue au cours du calcul du modèle, et ce pour différentes valeurs de points initiaux.
Quelqu'un peut m'aider ?

Merci en avance.

[A. D.]

Bonjour,

Je n'ai pas la réponse à votre question mais peut-être pourriez-vous poster un petit exemple de code reproductible sur lequel vous rencontrez l'erreur que vous mentionnez ? Cela aidera peut-être d'autres personnes à mieux vous répondre.

Bonne continuation


Cordialement,


A.D.

[robert duparvis]

Bonjour,

Voilà le code utilisé pour la résolution :

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

 nls( G1 ~ a1-b1*exp(-c1/Y1^d1), start=list(a1=1,b1=1,c1=1,d1=1))

Avec Y1 des dates et G1 des mesures :

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
Y1=[1]  1  2  2  3  4  4  4  5  5  6  6  6  6  6  7  7  7  8  8  8  8  9 10 10 11 11 12 12 12 12
G1=[1] 4.9995 4.9962 4.9962 4.9899 4.9812 4.9812 4.9812 4.9709 4.9709 4.9593 4.9593 4.9593 4.9593 4.9593 4.9469 4.9469 4.9469 4.9339 4.9339 4.9339 4.9339
[22] 4.9204 4.9067 4.9067 4.8927 4.8927 4.8787 4.8787 4.8787 4.8787

Quand je le lance j'obtiens les erreurs indiquées plus haut.

[juliatheric]

Bonjour,
L'erreur est due aux paramètres de départ. La stratégie habituelle est d'explorer analytiquement et graphiquement le modèle (trouver des points d'inflexion de son graphe, ses asymptotes, etc.) et se servir des constats faits pour initialiser les paramètres de départ. Si cela ne donne pas assez de succès, on peut procéder par la linéarisation du modèle et faire une régression linéaire sur le modèle linéarisé (voir ?lm). Dans le cas présent, ce serait quelque chose du genre :


(3) est de la forme de (4), donc la régression linéaire est applicable en (3). Après, on jongle avec les résidus pour remonter dans le modèle initial, mais ça reste de la gymnastique. C'est là où la bibliothèque drc s'avère puissante. Votre modèle est visiblement de la famille Weibull (voir ?drc::weibull1). Une fois ce constat fait, la méthode drc::drm peut faire la régression et ses résultats peuvent servir à trouver les paramètres de départ de nls. Le retour à nls n'est pas indispensable, il est simplement une des stratégies qu'on peut adopter lorsqu'on manipule nls. Alors, le code :

Code R :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
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
 
library("drc")
 
ourData <- list(
  x = c(1, 2, 2, 3, 4, 4, 4, 5, 5, 6, 6, 6, 6, 6, 7, 7, 7, 8, 8, 8, 8, 9, 10,
        10, 11, 11, 12, 12, 12, 12),
  y = c(4.9995, 4.9962, 4.9962, 4.9899, 4.9812, 4.9812, 4.9812, 4.9709, 4.9709,
        4.9593, 4.9593, 4.9593, 4.9593, 4.9593, 4.9469, 4.9469, 4.9469, 4.9339,
        4.9339, 4.9339, 4.9339, 4.9204, 4.9067, 4.9067, 4.8927, 4.8927, 4.8787,
        4.8787, 4.8787, 4.8787)
)
 
#drc::drm
drcObj <- drm(
  formula = y ~ x,
  data = ourData,
  fct = weibull1(names = c("b=-d", "c=a", "d=a-b", "e=exp(-log(c) / b)"))
)
summary(drcObj)
plot(drcObj)
 
#stats:nls
coeffStart <- list(
  a = drcObj$coefficients[2],
  b = drcObj$coefficients[2] - drcObj$coefficients[3],
  c = drcObj$coefficients[4] ^ (-drcObj$coefficients[1]),
  d = -drcObj$coefficients[1]
)
coeffStart <- sapply(coeffStart, function(xxx) drop(xxx))
names(coeffStart) <- letters[seq_len(4)]
print(coeffStart)
 
res <- nls(
  formula = y ~ a - b * exp(-c / x ^ d),
  data = ourData,
  start = coeffStart,
  control = nls.control(printEval = TRUE)
)
 
summary(res)
print(res$convInfo$finTol)
print(res$convInfo$finIter)
plot(res)

Avec ces paramètres de départ (la variable coeffStart), le modèle converge en seulement 5 itérations et le likelihood est assez faible.

Formula: y ~ a - b * exp(-c/x^d)

Parameters:
    Estimate Std. Error t value Pr(>|t|)    
a  4.744e+00  8.876e-03  534.51   <2e-16 ***
b -2.575e-01  9.045e-03  -28.47   <2e-16 ***
c  6.529e-03  9.761e-05   66.89   <2e-16 ***
d -1.852e+00  2.018e-02  -91.75   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.0002651 on 26 degrees of freedom

Number of iterations to convergence: 5 
Achieved convergence tolerance: 6.209e-06

[robert duparvis]

Merci beaucoup pour votre réponse détaillée.

[juliatheric]

De rien.

Addendum
nls possède un mécanisme de détermination automatique des paramètres de départ, il s'agit du concept de self start ; voir help("selfStart"). Dans la section "See Also" de cette page de doc, il y a actuellement 10 fonctions self-starting et on les utilise telles qu'elles sont définies : lors d'un appel il n'y a rien à changer à leurs arguments formels, ils deviennent les arguments effectifs. Dans le cas de votre modèle, 5 d'entre elles convergent et c'est à vous de voir si elles répondent à vos besoins. Si j'utilise des données pseudo-aléatoires avec les résultats qu'elles donnent, elles s'avèrent être dépourvues de capacités prédictives alors que drc s'en sort haut la main. À toutes fins utiles, voilà le code complet :

Code R :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
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
47
48
49
50
51
52
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
86
87
88
89
90
91
92
93
94
95
96
97
 
invisible(sapply(c("drc","foreach", "magrittr"),
                 function(z) library(z, character.only = TRUE)))
 
ourData <- list(
  x = c(1, 2, 2, 3, 4, 4, 4, 5, 5, 6, 6, 6, 6, 6, 7, 7, 7, 8, 8, 8, 8, 9, 10,
        10, 11, 11, 12, 12, 12, 12),
  y = c(4.9995, 4.9962, 4.9962, 4.9899, 4.9812, 4.9812, 4.9812, 4.9709, 4.9709,
        4.9593, 4.9593, 4.9593, 4.9593, 4.9593, 4.9469, 4.9469, 4.9469, 4.9339,
        4.9339, 4.9339, 4.9339, 4.9204, 4.9067, 4.9067, 4.8927, 4.8927, 4.8787,
        4.8787, 4.8787, 4.8787)
)
 
visFitResult <- function(fitObj = NULL) {
  oldPar <- par(las = 1)
  if (is.null(fitObj)) {
    plot(y ~ x, data = ourData, pch = 19L, col = "red")
  } else {
    plot(y ~ x, data = ourData, pch = 19L, col = "red")
    plot(fitObj, add = TRUE, broken = TRUE, lty = 2L, lwd = 1, col = "red")
  }
  grid(20L, col = "gray")
  par(oldPar)
  invisible(TRUE)
}
 
visFitResult()
 
#drc::drm
drcObj <- drm(
  formula = y ~ x,
  data = ourData,
  fct = weibull1(names = c("b=-d", "c=a", "d=a-b", "e=exp(-log(c) / b)"))
)
summary(drcObj)
visFitResult(drcObj)
 
 
#stats:nls: Self Start.
nlsObj <- list(
  resWei = nls(
    formula = y ~ SSweibull(x, Asym, Drop, lrc, pwr),
    data = ourData,
    control = nls.control(printEval = TRUE)
  ),
  resBiexp = nls(
    formula = y ~ SSbiexp(x, A1, lrc1, A2, lrc2),
    data = ourData,
    control = nls.control(printEval = TRUE)
  ),
  resFpl = nls(
    formula = y ~ SSfpl(x, A, B, xmid, scal),
    data = ourData,
    control = nls.control(printEval = TRUE)
  ),
  resLog = nls(
    formula = y ~ SSlogis(x, Asym, xmid, scal),
    data = ourData,
    control = nls.control(printEval = TRUE)
  ),
  resMiMe = nls(
    formula = y ~ SSmicmen(x, Vm, K),
    data = ourData,
    control = nls.control(printEval = TRUE)
  )
)
 
summary(nlsObj$resWei)
nlsObj$resWei$m$getPars()
visFitResult(nlsObj$resWei)
 
 
#Erreurs. Celle de drc vient de summary(drcObj) mais on peut aussi la calculer
#manuellement en se servant des résultats obtenus
sapply(nlsObj, function(x) x$convInfo$finTol) %>% c(drc = 3.010583e-4) %>% sort
 
 
#Comportements des modèles obtenus avec des résultats pseudo-aléatoires
randomize <- function(size = 30L) {
  oldPar <- par(las = 1L, mfcol = c(2L, 3L))
  arg <-rlnorm(size)
 
  plot(arg, drcObj$curve[[1]](arg),
       xlab = "drc", ylab = "", pch = 20L)
 
  foreach(i = nlsObj, 
          j = c("Weibull", "Biexponential", "Logistic-3par", "Logistic-4par",
                "Michaelis-Menten")) %do% {
                  plot(arg, i$m$predict(arg),
                       xlab = j, pch = 20L, ylab = "", new = TRUE)
                }
 
  par(oldPar)
  invisible(TRUE)
}
 
randomize()