Extraire les coefficients d’une équation de courbe pour chaque individu [Résolu]

Postez ici vos questions, réponses, commentaires ou suggestions - Les sujets seront ultérieurement répartis dans les archives par les modérateurs

Modérateur : Groupe des modérateurs

Aurélie Quinard
Messages : 14
Enregistré le : 12 Mar 2013, 11:32

Extraire les coefficients d’une équation de courbe pour chaque individu [Résolu]

Messagepar Aurélie Quinard » 15 Aoû 2017, 06:38

Bonjour,

Dans le cadre d’une étude sur les courbes de croissance, il me faut faire ajuster l’équation d’une courbe de croissance (ici le modèle MMF de la forme Weight = (a*b+c*Age^d)/(b+Age^d) ) sur un très grand nombre d’individus, plus clairement chaque individu doit avoir sa propre courbe avec ses propres coefficients a, b, c et d. Mon souci est de pouvoir récupérer les valeurs de ces 4 coefficients pour chaque individu.
Peut-être les informations complémentaires ci-dessous pourront vous être utiles :
Pour une analyse parallèle à celle que je souhaite faire à présent, j’ai pu développer un script me permettant de déterminer les coefficients pour la courbe d’une sous-population (script qui plus développé devrait m’amener à terme à déterminer si certains facteurs ont une influence sur les coefficients). J’ai essayé de réfléchir à un moyen de résoudre mon problème en partant de cet ébauche mais sans succès.

Voici le script – petite précision, les coefficients utilisés dans le « start » sont ceux correspondant à la courbe de croissance de la population entière ajustée grâce au logiciel CurveExpert, le script ayant besoin d’une base de coefficient pour tourner:


Code : Tout sélectionner

dataFemales<-read.table("GrowthCurveFemales.txt",header=T)

library(nlme)

MMFfunc<-function(a,b,c,d,Age){
(a*b+c*Age^d)/(b+Age^d)}

startFemales=c(a=51.59,b=1528,c=1342,d=1.8)

model <- nls(Weight~(a*b+c*Age^d)/(b+Age^d),start=startFemales, data=dataFemales)
summary(model)


Mes données étant de la forme ci-dessous, mais pouvant être remaniée si nécessaire :

Code : Tout sélectionner

Year   ID         Sex Age  Weight
2010   M10N00006   F   1   31
2010   M10N00006   F   2   33.2
2010   M10N00006   F   3   35.1
2010   M10N00006   F   4   43.1
2010   M10N00006   F   5   50
2010   M10N00006   F   6   61.4
2010   M10N00006   F   9   95
2010   M10N00006   F   10   105.1
2010   M10N00006   F   11   117.2
2010   M10N00006   F   12   126.6
2010   M10N00006   F   14   148
2010   M10N00006   F   76   818
2010   M10N00006   F   173   1002
2010   M10N00006   F   187   1073
2010   M10N00006   F   201   1082
2010   M10N00006   F   228   1200
2010   M10N00006   F   252   1210
2010   M10N00032   F   2   54
2010   M10N00032   F   3   58.2
2010   M10N00032   F   4   65.9
2010   M10N00032   F   5   73.4
2010   M10N00032   F   6   82.6
2010   M10N00032   F   8   103.7
2010   M10N00032   F   9   115.7
2010   M10N00032   F   10   128
2010   M10N00032   F   11   139.5
2010   M10N00032   F   13   166



Je vous remercie par avance pour votre aide.
Bonne journée à vous,

Aurélie

Stéphane Adamowicz
Messages : 206
Enregistré le : 07 Mar 2012, 10:13
Contact :

Re: Extraire les coefficients d’une équation de courbe pour chaque individu

Messagepar Stéphane Adamowicz » 17 Aoû 2017, 10:10

Bonjour,

vous avez raison de penser que vos données peuvent être remaniées. Vous avez noté "F" les femelles. Or, F est une valeur logique que R interprète comme signifiant FALSE. Cela n'a pas d'importance pour votre problème, mais mérite d'être signalé. Ce qui m'amène à une deuxième remarque. Il est impératif de faire systématiquement le summary d'un dataframe. Cela permet de faire un examen rapide des données, et vous vous seriez alors aperçue de la bévue. J'ai remplacé F par Female dans votre extrait de données, qui devient donc :

Code : Tout sélectionner

dataFemales <- structure(list(Year = c(2010L, 2010L, 2010L, 2010L, 2010L, 2010L,
2010L, 2010L, 2010L, 2010L, 2010L, 2010L, 2010L, 2010L, 2010L,
2010L, 2010L, 2010L, 2010L, 2010L, 2010L, 2010L, 2010L, 2010L,
2010L, 2010L, 2010L), ID = structure(c(1L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L,
2L, 2L, 2L, 2L, 2L), .Label = c("M10N00006", "M10N00032"), class = "factor"),
    Sex = c("Female", "Female", "Female", "Female", "Female",
    "Female", "Female", "Female", "Female", "Female", "Female",
    "Female", "Female", "Female", "Female", "Female", "Female",
    "Female", "Female", "Female", "Female", "Female", "Female",
    "Female", "Female", "Female", "Female"), Age = c(1L, 2L,
    3L, 4L, 5L, 6L, 9L, 10L, 11L, 12L, 14L, 76L, 173L, 187L,
    201L, 228L, 252L, 2L, 3L, 4L, 5L, 6L, 8L, 9L, 10L, 11L, 13L
    ), Weight = c(31, 33.2, 35.1, 43.1, 50, 61.4, 95, 105.1,
    117.2, 126.6, 148, 818, 1002, 1073, 1082, 1200, 1210, 54,
    58.2, 65.9, 73.4, 82.6, 103.7, 115.7, 128, 139.5, 166)), .Names = c("Year",
"ID", "Sex", "Age", "Weight"), row.names = c(NA, -27L), class = "data.frame")

summary(dataFemales)
      Year              ID         Sex                 Age             Weight     
 Min.   :2010   M10N00006:17   Length:27          Min.   :  1.00   Min.   :  31.0 
 1st Qu.:2010   M10N00032:10   Class :character   1st Qu.:  4.50   1st Qu.:  59.8 
 Median :2010                  Mode  :character   Median :  9.00   Median : 105.1 
 Mean   :2010                                     Mean   : 46.85   Mean   : 304.4 
 3rd Qu.:2010                                     3rd Qu.: 13.50   3rd Qu.: 157.0 
 Max.   :2010                                     Max.   :252.00   Max.   :1210.0 
>


Pour sélectionner un individu dans un dataframe, utilisez l'argument subset de la fonction nls, et pour extraire les coefficients, utilisez coef :

Code : Tout sélectionner

> startFemales=c(a=51.59,b=1528,c=1342,d=1.8)
> model <- nls(Weight~(a*b+c*Age^d)/(b+Age^d),start=startFemales, data=dataFemales, subset= ID=="M10N00006")
> summary(model)
Formula: Weight ~ (a * b + c * Age^d)/(b + Age^d)

Parameters:
   Estimate Std. Error t value Pr(>|t|)   
a   16.2396    20.8564   0.779    0.450   
b  356.7994   226.2898   1.577    0.139   
c 1305.5352    71.9026  18.157 1.28e-10 ***
d    1.4399     0.1752   8.218 1.66e-06 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 34.31 on 13 degrees of freedom

Number of iterations to convergence: 9
Achieved convergence tolerance: 6.11e-06

>  coef(model)
          a           b           c           d
  16.239563  356.799399 1305.535207    1.439901


Enfin, si vous voulez remplir automatiquement un dataframe avec les coefficients de chaque individu, préparez d'abord ce dataframe avec une ligne pour chaque individu et des NA pour les coefficients :

Code : Tout sélectionner

> individus <- unique(dataFemales$ID) ; dummy <- rep(NA, length(individus))
> resultat <- data.frame(ID=individus, a=dummy, b=dummy, c=dummy, d=dummy)
> resultat
         ID  a  b  c  d
1 M10N00006 NA NA NA NA
2 M10N00032 NA NA NA NA


Puis calculez dans une boucle (ce n'est qu'une solution parmi d'autres) les coefficients de chaque individu :

Code : Tout sélectionner

> compteur <- 0
> for(i in individus)
+     {model <- nls(Weight~(a*b+c*Age^d)/(b+Age^d),start=startFemales, data=dataFemales, subset= ID==i)
+     compteur <- compteur + 1
+     resultat[compteur,] <- c(i,coef(model))}
> resultat
         ID                a                b                c                d
1 M10N00006 16.2395632915754 356.799398578187 1305.53520677308 1.43990147008223
2 M10N00032  49.030270721997 405.497708156275 581.118646640876 1.84670902621133


Tout ceci en priant pour que les valeurs de départ conviennent à tous les individus ...

cordialement,

Stéphane
Stéphane Adamowicz
INRA, UR 1115 Plantes et Systèmes de Culture Horticoles (PSH)
domaine St Paul, site agroparc
84914 Avignon, cedex 9

Aurélie Quinard
Messages : 14
Enregistré le : 12 Mar 2013, 11:32

Re: Extraire les coefficients d’une équation de courbe pour chaque individu

Messagepar Aurélie Quinard » 20 Aoû 2017, 10:23

Bonjour Stéphane,

je vous remercie pour votre réponse. J'ai pu tester votre proposition ce matin.

Le code pour la récupération des coefficients pour un seul individu marche très bien. Lorsque j'ai comparé ce que R me sortait avec ce que CurveExpert me calculait, les coefficients étaient légèrement différents mais en visualisant sur un graphique, les deux courbes se chevauchent parfaitement bien, donc le calcul est bon.

Petit hic, lorsque j'ai voulu testé le second code sur mon jeu de donné incluant tous mes individus (5718 indiv), un message d'erreur est apparu:

Code : Tout sélectionner

Error in nls(Weight ~ (a * b + c * Age^d)/(b + Age^d), start = startTestM,  :
  singular gradient


De ce que j'ai pu trouver sur internet, cela pourrait être dû à des estimateurs de départ (startTestM=c(a=51.59,b=1528,c=1342,d=1.8)) trop mal choisis.

En testant sur un sous-jeu d'une trentaine d'individus, le code marchait bien. J'ai ensuite testé sur des individus aux croissances extrêmes (plus petits poids ou plus grands poids à l'âge adulte) ou avec des points qui sortaient de leur courbe mais là, encore le code réussi à calculer des coefficients. Je ne vois donc pas pour le moment ce qui pourrait poser problème.

Merci également pour votre conseil sur ma façon de désigner les femelles, je ne savais pas du tout que l'utilisation de la lettre "F" pouvait poser des problèmes.

Bonne journée,

Aurélie

Stéphane Adamowicz
Messages : 206
Enregistré le : 07 Mar 2012, 10:13
Contact :

Re: Extraire les coefficients d’une équation de courbe pour chaque individu

Messagepar Stéphane Adamowicz » 21 Aoû 2017, 09:02

Bonjour,

je vous l'avais bien dit qu'il fallait prier (le dieu des algorithmes) ...

Pour ce type d'erreur, on peut effectivement tenter de changer les valeurs de départ, avec le risque que de nouvelles valeurs ne conviennent plus aux autres individus.

On peut aussi changer d'algorithme d'optimisation. Par défaut, l'algorithme utilisé est celui de Gauss-Newton (voir l'aide de la fonction nls). On peut essayer un autre algortihme proposé par cette fonction :

Code : Tout sélectionner

model <- nls(Weight~(a*b+c*Age^d)/(b+Age^d),start=startFemales, data=dataFemales, subset= ID==i, algorithm= "port")


Autre message d'erreur que vous pouvez rencontrer :

Code : Tout sélectionner

Convergence failure: iteration limit reached without convergence (50)


Il signifie que l'algorithme n'a pas trouvé de solution avant d'avoir atteint le nombre maximal d'essais (fixé à 50 par défaut). Vous pouvez augmenter cette dernière valeur. Par exemple :

Code : Tout sélectionner

model <- nls(Weight~(a*b+c*Age^d)/(b+Age^d),start=startFemales, data=dataFemales, subset= ID==i, algorithm= "port", control=list(maxiter=500))


Bonne chance !

Stéphane
Stéphane Adamowicz

INRA, UR 1115 Plantes et Systèmes de Culture Horticoles (PSH)

domaine St Paul, site agroparc

84914 Avignon, cedex 9

Aurélie Quinard
Messages : 14
Enregistré le : 12 Mar 2013, 11:32

Re: Extraire les coefficients d’une équation de courbe pour chaque individu

Messagepar Aurélie Quinard » 23 Aoû 2017, 04:53

Bonjour Stéphane,

merci pour le temps que vous avez pris à me répondre. Après avoir testé vos suggestions, des erreurs sont de nouveau apparu avec mon jeu de données principales. Néanmoins, en cherchant encore de mon côté, j'ai pu découvrir une fonction remplaçant nls() qui m'a permis de faire tourner votre code sans soucis. Il s'agit de la fonction nlsLM() du package {minpack.lm}. Ce qui me donne donc le code suivant:

Code : Tout sélectionner

library(minpack.lm)

individus <- unique(dataHealthy$ID) ; dummy <- rep(NA, length(individus))
resultat <- data.frame(ID=individus, a=dummy, b=dummy, c=dummy, d=dummy)
resultat

compteur <- 0
for(i in individus)
      {model <- nlsLM(Weight~(a*b+c*Age^d)/(b+Age^d),start=list(a=46,b=1288,c=1599,d=1.7), data=dataHealthy, subset= ID==i)
      compteur <- compteur + 1
     resultat[compteur,] <- c(i,coef(model))}
resultat


Encore merci pour votre aide.
En vous souhaitant une belle journée,

Aurélie


Retourner vers « Questions en cours »

Qui est en ligne

Utilisateurs parcourant ce forum : Google [Bot] et 1 invité