utilisation des vecteurs d'indice

[herve chapuis]

Bonjour,

je pars d'un tableau contenant les informations suivantes :

Code : Tout sélectionner

> TOT[1:5]
   famille ind pere mere Y
1       19   1    7    7 1
2        8   2    3    2 1
3       16   3    6    4 1
4        9   4    3    3 1
5       17   5    6    5 1


autrement dit un identifiant, un numéro de famille, un père, une mère et une variable binaire Y.

Pour obtenir ce tableau j'ai au préalable établi mes règles de correspondance entre famille et parents :

Code : Tout sélectionner

> Nfamille
[1] 27
> dam_fam
 [1] 1 2 3 1 2 3 1 2 3 4 5 6 4 5 6 4 5 6 7 8 9 7 8 9 7 8 9
> sire_fam
 [1] 1 1 1 2 2 2 3 3 3 4 4 4 5 5 5 6 6 6 7 7 7 8 8 8 9 9 9

Dans l'exemple j'ai 27 familles, à partir de 9 pères et 9 mères.

Je génère mon tableau en tirant pour chacun de mes NIND0 individus une famille de façon aléatoire :

Code : Tout sélectionner

floor(Nfamille*runif(n=NIND0,min=0,max=1))+1


Je désire ensuite synthétiser l'information dans un tableau qui contienne une ligne par famille, avec le nombre de (Y=1) et le nombre de (Y=0).

Code : Tout sélectionner

table(TOT$Y,TOT$famille)
   
    1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 21 23 24 25 26 27
  0 0 1 1 1 2 0 3 1 0  1  2  1  0  2  0  0  1  0  0  1  1  1  0  0  2
  1 5 1 2 2 1 4 0 3 3  1  0  0  2  0  3  2  4  1  2  0  5  2  1  2  3

on voit qu'il peut y avoir des familles non représentées.
J'aimerais utiliser cette liste de familles comme vecteur d'indices.

Sauf que colnames (ou rownames si je prends la transposée) n'est pas numérique :

Code : Tout sélectionner

> colnames(table(TOT$Y,TOT$famille))
 [1] "1"  "2"  "3"  "4"  "5"  "6"  "7"  "8"  "9"  "10" "11" "12" "13" "14" "15"
[16] "16" "17" "18" "19" "21" "23" "24" "25" "26" "27"

du coup je ne peux l'utiliser comme vecteur d'indice.

Code : Tout sélectionner

sire_fam[colnames(table(TOT$Y,TOT$famille))]
 [1] NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA


Peut être que ce n'est pas la bonne façon de procéder ?

Cela étant, ensuite je voudrais utiliser les valeurs estimées pour les pères et les mères apres un modèle mixte

Code : Tout sélectionner

lm3<-glmer(cbind(vivant,mort)~1+(1|pere)+ (1|mere)+(1|famille),data=TOT,family=binomial)
EBV3<-ranef(lm3)
pere3<-as.vector(unlist(EBV3$pere))


Selon la taille du dispositif il n'est pas dit que tous les pères seront représentés et donc estimés. C'est pourquoi je désire pouvoir utiliser le nom des variables et non simplement l'ordre.
En effet, si j'ai 80 pères estimés sur 80 possibles, je sais dire que pere3[25] correspond au père n°25. Si je n'en n'ai que 75 représentés dans mon échantillon simulé, il me faut un moyen pour attribuer la bonne valeur au bon père.

Suis-je clair ?

Merci

[David Robelin]

Bonjour,

Pour ton premier problème, de nom de colonne numérique, tu peux les transformer en nombre en utilisant la fonction as.numeric(...)

David

[herve chapuis]

Merci

Je pensais avoir essayé et que cela ne marchait pas.
J'ai testé sans y croire (pour démontrer que ça n'allait pas le faire) et ... ça marche !!!
:wink:

Au moins pour la partie 1 de ma question.

[David Robelin]

Bonjour,

Je ne suis pas tout à fait sur d'avoir compris le deuxième point.

Il me semble que les identifiants sont données par la commande :

Code : Tout sélectionner

rownames(EBV3$pere)


et qu'il est possible d'accéder à la ligne désirée de la façon suivante :

Code : Tout sélectionner

EBV3$pere["25",]


Ne pas oublier les guillemets qui indiquent que tu veux le père numéro "25" et non la 25ème valeur.
En espérant que ça répond à la question.

[herve chapuis]

Je relance la question. En effet je n'ai pu m'y consacrer avec autant d'énergie que je le souhaitais.
J'ai un peu de temps libre donc me revoilà...

Je reprends :
Je génère une liste d'animaux soumis à challenge ==> pedi0, et une liste d'animaux candidats à la sélection ==> pedi1.

Les deux jeux d'animaux sont issus des mêmes parents (mêmes familles).

L'idée est d'utiliser l'information apportée par les animaux du challenge pour évaluer les pères et les mères. Ce qui permet d'attribuer une valeur génétique sur ascendance aux candidats.

Pour cela j'utilise lmer

Code : Tout sélectionner

lm0<-lmer(survie~1+(1|pere)+ (1|mere),data=TOT0,family=binomial)

L'idée est donc d'utiliser les valeurs estimées pour les pères et les mères

Code : Tout sélectionner

EBV0pere<-as.vector(EBV0$pere)
EBV0mere<-as.vector(EBV0$mere)

Mon souci est que sur les 100 pères possibles, je n'en n'ai que 99 de représentés. Grrr je n'arrive pas copier coller ma sélection.
Bon considérons que je passe du père 61 au père 63.
Donc il me manque le père 62.
Pour ceux qui se demanderaient pourquoi je n'ai pas davantage soigné mon échantillonnage, il s'agit d'un génotypage a posteriori. Les œufs sont fécondés selon un plan factoriel puis les larves éclosent et tout est mélangé. Selon le nombre d'animaux tirés parmi les candidats et les challengés, on n'est absolument pas certain de représenter tous les parents.


Or parmi mes candidats, ce père 62 est représenté. J'ai donc besoin de calculer une valeur génétique pour tous ses descendants et une solution consiste à considérer la valeur du père 62 comme nulle (i.e. non informative dans la sommation "père + mère").
Pour mes candidats je cherche à calculer la valeur génétique sur ascendance, sachant que le père est en colonne 4 dans pedi1 et la mère en colonne 5.

Code : Tout sélectionner

EBV1p<-EBV0pere[pedi1[,4]]
EBV1m<-EBV0mere[pedi1[,5]]
EBVcan<-EBV1p + EBV1m

Le problème c'est que ma formule est fausse dès qu'un parent au moins est non représenté, tout est décalé.
Ca passe du 61 au 63 sans sourciller mais les descendants du père 100 n'ont pas de valeur estimée (car ma liste n'a que 99 pères). C'est donc faux pour 38 pères sur 99.

Idéalement il faudrait créer une matrice (ou un tableau) de dimension npere * 2 (et nmere *2).
où npere et nmere représentent les nombres de parents théoriquement représentés. Puis initialiser cette matrice de 1 à npere en colonne 1 et avec des zero en colonne 2.
Enfin je "colle" (mais comment ????) mes solutions du modèle mixte de telle sorte que ma matrice contienne soit un 0 (si le parent n'est pas représenté) soit sa valeur estimée.
Est-ce que c'est clairement posé ?
Et surtout qui aurait une solution ?

Merci !

[David Robelin]

Bonjour,


Peut-être qu'il vaut mieux éviter de passer par les deux variables : EBV0pere et EBV0mere (j'imagine que tu utilise la fonction 'ranef')

Par exemple :

Code : Tout sélectionner

EBV0 <- ranef(lm0)
EBVcan <-  EBV0[pedi1[,4],"pere"] + EBV0[pedi1[,5],"mere"]


Pas facile de tester sans un exemple de données (qu'est-ce qui est dans pedi1 ?).

[herve chapuis]

bonjour et merci,

d'abord la première question : qu'y a t il dans pedi1 ?

Code : Tout sélectionner

1> pedi1[1:5,]
      famille ind sexe pere mere   true_res
 [1,]     513   1    2   52   53 -1.5822285
 [2,]     385   2    1   39   35 -4.5450651
 [3,]      27   3    2    3    7  0.1415672
 [4,]     571   4    1   58   51 -4.9827598
 [5,]      21   5    2    3    1 -1.4416595


et accessoirement dans pedi0

Code : Tout sélectionner

1> pedi0[1:5,]
      famille ind pere mere   true_res
 [1,]      36   1    4    6 -2.1803001
 [2,]     290   2   29   30 -0.7211150
 [3,]     645   3   65   65  0.8746883
 [4,]     946   4   95   96  1.4700812
 [5,]     622   5   63   62  4.1573899


d'où l'on déduit TOT0, le fichier soumis à l'analyse de lmer :

Code : Tout sélectionner

1> TOT0[1:5,]
   famille pere mere survie
1       36    4    6      0
2      290   29   30      0
3      645   65   65      1
4      946   95   96      0
5      622   63   62      1


survie étant le résultat de la transformation en variable binaire d'une variable continue obtenue en ajoutant un bruit de fond à la variable "true_res" (susceptible d'être transmise à la descendance).

Je ne peux utiliser la suggestion car j'obtiens le message suivant :

Code : Tout sélectionner

EBVcan<-EBV0[pedi1[,4],"pere"]
Erreur dans EBV0[pedi1[, 4], "pere"] : nombre de dimensions incorrect


Et oui j'obtiens bien EBV0 par la fonction ranef.
A ce propos je viens de remarquer un truc bizarre : Si j'ai 2000 animaux dans pedi0 j'obtiens :

Code : Tout sélectionner

1> lm0
Generalized linear mixed model fit by the Laplace approximation
Formula: survie ~ 1 + (1 | pere) + (1 | mere)
   Data: TOT0
  AIC  BIC logLik deviance
 2746 2763  -1370     2740
Random effects:
 Groups Name        Variance Std.Dev.
 pere   (Intercept) 0.10011  0.31640
 mere   (Intercept) 0.15816  0.39769
Number of obs: 2000, groups: pere, 100; mere, 100

Fixed effects:
            Estimate Std. Error z value Pr(>|z|) 
(Intercept) -0.11649    0.06871  -1.695     0.09 .
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1


Autrement dit j'ai bien une variance père et une variance mère.
Mais si je n'ai que 500 animaux dans le challenge :

Code : Tout sélectionner

1> lm0
Generalized linear mixed model fit by the Laplace approximation
Formula: survie ~ 1 + (1 | pere) + (1 | mere)
   Data: TOT0
   AIC BIC logLik deviance
 694.3 707 -344.2    688.3
Random effects:
 Groups Name        Variance Std.Dev.
 mere   (Intercept) 0.00000  0.00000
 pere   (Intercept) 0.23072  0.48033
Number of obs: 500, groups: mere, 100; pere, 99

Fixed effects:
            Estimate Std. Error z value Pr(>|z|)
(Intercept)  -0.1369     0.1039  -1.317    0.188


pfff! la variance mère a disparu...
En plus, vu qu'il s'agit d'un schéma factoriel, les pères et mères devraient avoir une contribution à peu près équivalente.

[David Robelin]

Si j'ai bien suivi, seule la variable TOT0 est utile pour l'exemple. Voilà une proposition de code, avec un mini-exemple de données (4 individus, 2 pères, 2 mères)

Code : Tout sélectionner

#Lecture des donnees exemples
toto = read.table(textConnection("
    famille ind sexe pere mere   survie
[1,]     513   1    2   6   8 0
[2,]     385   2    1   6   8 1
[3,]      27   3    2   7   9 0
[4,]     571   4    1   7   9 1"))

#Ajustement de la regression
lm0<-lmer(survie~1+(1|pere)+ (1|mere),data=toto,family=binomial)

#Extraction des effets aléatoires
EBV0<-ranef(lm0)

#Ajout de deux colonnes au data.frame 'toto' (attention petite subtilite explique plus bas)
toto$EBV1p = EBV0$pere[as.character(toto$pere),1]
toto$EBV1m = EBV0$mere[as.character(toto$mere),1]

# Affichage du debut du jeu donnees pour voir ou on en est
head(toto)

# On remplace les NA par des 0
is.na(toto$EBV1p)<-0
is.na(toto$EBV1m)<-0

# Ajout d'une nouvelle colonne contenant la somme
toto$EBVcan = toto$EBV1p+toto$EBV1m

# Affichage du resultat
toto


La petite subtilité est l'utilisation de 'as.character'. En effet, 'ranef' produit une liste de matrices dont les lignes sont nommés par le numéro du père (ou de la mère). Attention à ne pas confondre :

Code : Tout sélectionner

EBV0$pere[18,1]     # => renvoie la 18eme ligne
EBV0$pere["18",1]     # => renvoie la ligne nomme "18"


En espèrant que ca aide,

David

[herve chapuis]

Merci !
On tourne autour, je veux croire que l'on s'approche... mais on est encore à côté. Hélas.

Je n'ai pas dû être très clair.
Le TOT0 me sert à faire tourner le lmer pour obtenir une valeur pour les pères et les mères.
Ensuite je dois attribuer une valeur génétique aux candidats qui sont dans pedi1, que j'ai déclaré comme une matrice (s'il faut en faire data.frame et si ça change quelque chose en bien, je suis open).
Le TOT0 c'est un data.frame à partir des colonnes idoines de pedi0, qui est du même type que pedi1.

Code : Tout sélectionner

pedi1<-matrix(nrow=NIND1,ncol=6)
colnames(pedi1)<-c("famille","ind","sexe","pere","mere","true_res")

si je prends un petit exemple:

Code : Tout sélectionner

pedi1[1:5,]
      famille ind sexe pere mere   true_res
 [1,]     870   1    1   87   90  -1.149654
 [2,]     836   2    2   84   86   0.314609
 [3,]     951   3    1   96   91   4.042575
 [4,]     610   4    1   61   70 -16.119775
 [5,]     745   5    1   75   75  22.568000


et si je regarde les valeurs pères obtenues par lmer sur le fichier TOT0 j'ai :
* père 87 : -0.282
* père 84 : -0.621
* père 96 : -0.070
* père 61 : -0.369
* père 75 : 0.174

Normalement ça devrait marcher avec

Code : Tout sélectionner

EBV1p = EBV0$pere[as.character(pedi1[,3]),1]

sauf que non :

Code : Tout sélectionner

1> EBV1p[1:5]
[1] -0.1076915  0.5413462 -0.1076915 -0.1076915 -0.1076915


et en fait si je creuse un peu :

Code : Tout sélectionner

as.character(pedi1[1:10,3])
 [1] "1" "2" "1" "1" "1" "1" "2" "1" "2" "1"


au lieu d'obtenir :
"87" "84" "96" "61" "75"

Peut être que la clé est dans cette direction ?

[dicko ahmadou]

Salut,
Je ne crois pas pouvoir apporter une solution, je voulais juste te dire que la colonne pere c'est la "quatrième" et non la "troisième".
cordialement

[herve chapuis]

Houps !!
:oops:

Merci car je viens de vérifier en faisant bien attention et déjà la correspondance a l'air bien meilleure.

Je dois creuser