Discussion :

[Aurélie Q.]

Bonjour à tous,

je reviens vers vous pour un nouveau problème concernant mon analyse sur les courbes de croissances, un peu plus compliqué cette fois.
Je possède un jeu de données avec de nombreux individus. Pour chacun de ces individus, je possède le poids pour différents âges suivant l’extrait de tableau ci-dessous :

Code :

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

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
ID         Sex   Age   Weight
M15N00010   Mal   0   41.2
M15N00010   Mal   1   39.3
M15N00010   Mal   2   42.3
M15N00010   Mal   3   46.3
M15N00010   Mal   4   53.4
M15N00010   Mal   5   71.5
…   …   ..   …
M15N00012   Mal   0   37.3
M15N00012   Mal   1   34.6
M15N00012   Mal   2   36.3
M15N00012   Mal   3   39.7
M15N00012   Mal   4   47.5
M15N00012   Mal   5   47.3
…   …   …   …

Sur ce jeu de données, j’ai réussi à obtenir un script me permettant de faire fitter la courbe d’un modèle de croissance sur chaque individu et de récupérer les paramètres du modèle.
Le script :

Code :

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
dataTestSubM<-read.table("TestSubsetM.txt",header=T) 
 
MMFfunc<-function(a,b,c,d,Age){(a*b+c*Age^d)/(b+Age^d)}
 
startTestM=c(a=60,b=2376,c=1904,d=1.8)
 
individus <- unique(dataTestSubM$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 <- nls(Weight~(a*b+c*Age^d)/(b+Age^d),start=startTestM, data=dataTestSubM, subset= ID==i,trace=TRUE)
      compteur <- compteur + 1
     resultat[compteur,] <- c(i,coef(model))}
resultat

Ce qui me donne alors la sortie:

Code :

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

1
2
3
ID                a                b                c                d
1 M15N00010 16.2395632915754 356.799398578187 1305.53520677308 1.43990147008223
2 M15N00012  49.030270721997 405.497708156275 581.118646640876 1.84670902621133

J’aimerai à présent obtenir les coordonnées du point d’inflexion, encore une fois pour chaque individu. Pour d’autres modèles que celui-ci que j’utilise, ce point d’inflexion est directement caractérisé par un des paramètres. Malheureusement pour moi, le modèle qui correspond le mieux à mes données est différent. Je pensais donc passer par la méthode de la dérivée seconde pour retrouver ce point d’inflexion.

Une idée que j’ai eu (peut-être en aurait vous une meilleur ou saurez-vous m’aider à développer celle-ci) était de coder un nouveau dataframe avec en colonne : le nom de l’individu, l’âge de l’individu de 0 à 365 et le poids « théorique » de l’individu pour chaque âge, calculé grâce à l’équation de mon modèle et les coefficients spécifiques récupérés avec le code précédent.
Je pensais ensuite récupérer les coordonnées du point lorsque la dérivée seconde s’annule.
J’ai tenté de me débrouiller mais je ne parviens pas d’entrée de jeu à créer le nouveau dataframe et en essayant de réfléchir à la manière de coder le reste, je n’ai pas eu plus de succès.

Merci pour l’aide que vous pourrez m’apporter,

Aurélie

[tototode]

Bonjour,

est-ce qu'un point d'inflexion autour de 18 pour le premier ID te semble cohérent ?

cordialement

[Aurélie Q.]

Bonjour,

est-ce qu'un point d'inflexion autour de 18 pour le premier ID te semble cohérent ?

cordialement

Bonjour,

j'ai tracé la courbe du premier ID avec le point pour x = 18 pour tenter de voir si visuellement cela faisait sens (n'ayant aucun moyen plus précis de déterminer le point d'inflexion). Cela m'a l'air d'être bon. La courbe paraît bien change de concavité vers ce point-ci.
Je te joins mon graphe ci-dessous:



Cordialement.

[tototode]

re,

alors voilà comment j'ai fait pour la première ligne de ID :

Code :

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
#fonction qui permet de déterminer la dérivée énième d'une fonction
DD <- function(expr, name, order = 1) {
  if(order < 1) stop("'order' must be >= 1")
  if(order == 1) D(expr, name)
  else DD(D(expr, name), name, order - 1)
}
 
# pour remplacer les paramètres a, b, c et d par les valeurs de la ième ligne du tableau ID (tu peux boucler dessus)
i <- 1
exp_mod <- substitute((a*b+c*x^d)/(b+x^d), ID[i,])
 
f2 <- DD(exp_mod, "x", 2) # dérivée seconde
f3 <- DD(exp_mod, "x", 3)# dérivée d'ordre 3 (si verif > 0)
 
d2 <- function(val) {
  res <- eval(f2, list(x = val))
  res
}
 
d3 <- function(val) {
  res <- eval(f3, list(x = val))
  res
}
 
# estimation du point d'inflexion
r1 <- uniroot(d2, c(1, 365))$root

Cordialement

[Aurélie Q.]

Bonjour tototode,

j'ai essayé de faire tourner le script que tu me proposes. Néanmoins, il coince.

A la ligne:

Code :

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

exp_mod <- substitute((a*b+c*x^d)/(b+x^d), ID[,i])

Il me renvoie ce message d'erreur:

Code :

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

Error in `[.default`(ID, , i) : incorrect number of dimensions

J'ai tenté de remplacer dans la ligne précédente "i <- 1", le 1 par le nombre d'individus, puis par le nombre de variables... mais non.

Bonne journée.

[mgdondon]

Bonjour,

Regardez bien le programme de tototode. Vous avez mis la variable i en indice de colonne et non de ligne.

Cordialement,

[Aurélie Q.]

Bonjour,

Regardez bien le programme de tototode. Vous avez mis la variable i en indice de colonne et non de ligne.

Cordialement,

Il est vrai que dans mon précédent message, à force de tester différentes choses, la virgule était passée devant le "i" au lieu de derrière. J'ai tenté de refaire tourner le script pour vérifier mais j'obtiens toujours la même erreur.

[tototode]

Bonjour,

Pour répondre efficacement a ta question il faudrait connaître exactement le message d'erreur que tu obtiens.
Pour que le code fonctionne il faut que ID soit un data.frame et que les colonnes a, b, c et d soient bien des numeric, sinon ça ne marchera pas.

Code :

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
str(ID)
'data.frame':   2 obs. of  5 variables:
 $ ID: Factor w/ 2 levels "M15N00010","M15N00012": 1 2
 $ a : num  16.2 49
 $ b : num  357 405
 $ c : num  1306 581
 $ d : num  1.44 1.85
 
r1 <- numeric(nrow(ID))
for (i in 1:nrow(ID)) {
  exp_mod <- substitute((a*b+c*x^d)/(b+x^d), ID[i,])
  f2 <- DD(exp_mod, "x", 2)
  d2 <- function(val) {
    res <- eval(f2, list(x = val))
    res
  }
  pt <- try(uniroot(d2, c(1, 365))$root)
  if (inherits(pt, "try-error")) {
    pt[i] <- NA
    next
    }
  r1[i] <- pt
  }
r1
[1] 18.02706 13.39907

Cordialement

[Aurélie Q.]

Bonjour totode,

effectivement, les colonnes a,b,c et d n'étaient pas en numeric. A partir de là, après les avoir transformé en numeric, j'obtenais toujours le message d'erreur précédent. Néanmoins en remplaçant "ID" par le data.frame "resultat" dans ton code, le script tourne et je retrouve les valeurs que tu avais trouvé.

Du coup, pour le moment je peux obtenir l'abscisse du point d'inflexion. Saurais-tu s'il est possible de récupérer l'ordonnée également?

Bien cordialement,

Aurélie

[tototode]

Bonjour,

pour récupérer les 'y' et pour la correction d'une coquille :

Code :

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
r1 <- numeric(nrow(ID))
y <- numeric(length(r1))
 
for (i in 1:nrow(ID)) {
  exp_mod <- substitute((a*b+c*x^d)/(b+x^d), resultat[i,])
  f2 <- DD(exp_mod, "x", 2)
  d2 <- function(val) {
    res <- eval(f2, list(x = val))
    res
  }
  pt <- try(uniroot(d2, c(1, 365))$root)
  if (inherits(pt, "try-error")) {
    r1[i] <- NA
    y[i] <- NA
    next
    }
  r1[i] <- pt
  y[i] <- eval(exp_mod, list(x = pt))
  }

Cordialement

[Aurélie Q.]

Bonjour,

je viens de tester la suite du script. Tout marche parfaitement bien.

Je te remercie encore pour ton aide précieuse.

Aurélie