loop de regression et stockage des sorties

Version imprimable

22/06/2017, 20h04
Tsüki

1 pièce(s) jointe(s)

loop de regression et stockage des sorties

Bonjour à tous,
J’aimerais faire un loop de régression. En effet je veux expliquer la variable (duree_croissance) en première colonne (j’ai créé un petit fichier en pièce joint pour exemple) en fonction des autres variables (ici colonnes 2 à 6) (en vrai j’en ai plus de 200 c’est pourquoi je veux faire un loop).
Je veux donc faire des régressions et plot avec à chaque fois les données de la colonne 1 « duree_croissance » avec la colonne 2, puis avec la colonne 3, puis la 4, etc… J’utilise les fonctions « lm », « anova », « summary » et « sigma » afin d’avoir les résultats désirés. Sachant que j'ai des "NA" dans mon fichier donc j'en ai mis dans le fichier en pièce jointe, je pense que ça va poser probleme et qu'il faut éviter les NA pour les régressions.

Est-il possible de faire un loop pour avoir donc toutes ces reg (col 1 ~ col2, col 1 ~ col3, col1 ~ col4 , etc..) sans avoir à écrire toutes les régressions lm et que la boucle prenne elle-même les colonnes l’une après l’autre pour faire des régressions en fonction de la colonne 1 puis à chaque fois un plot avec la droite de régression, l’équation et le R² ?
J’aimerais aussi créer un nouveau data. frame avec les coefficients suivants : a, b, R² adjusted, R², N, pvalue, RMSE et n (nombre de ligne) de toutes les régressions. Je sais comment trouver ces paramètres et les imprimer pour une seule régression mais pas pour une boucle.

Merci par avance et bonne soirée.Pièce jointe 288086

Citation:

Est-il possible de faire un loop pour avoir donc toutes ces reg (col 1 ~ col2, col 1 ~ col3, col1 ~ col4 , etc..)

Citation:

J’aimerais aussi créer un nouveau data. frame avec les coefficients suivants : a, b, R² adjusted, R², N, pvalue, RMSE et n (nombre de ligne) de toutes les régressions. Je sais comment trouver ces paramètres et les imprimer pour une seule régression mais pas pour une boucle.

Une idée de solution
Code:

1 2 3 4 5 6 7 8 9 10 11 12 13 14 noms <- names( donnees) # donnees est le data.frame des donnees df.resultat <- data.frame( R2=rep( NA, length( noms) - 1), P=NA, ...) rownames( df.resultat) <- noms[-1] list.lm <- setNames( vector( "list", length( noms) - 1)), noms[-1]) # On suppose des regressions de la premiere colonne sur les autres for( nom in noms[-1]) { form <- formula( paste( noms[1], nom, sep="~") lm.prov <- lm( form, df.resultat) # na.omit par defaut list.lm[[nom]] <- lm.prov sum.prov <- summary( lm.prov) df.resultat[nom,"R2"] <- sum.prov$r.squared # meme chose pour les autres valeurs }
Par contre, je ne conseille pas le plot dans la boucle car la boucle va trop vite et dans le cas d'une nouvelle fenêtre par régression, il va y avoir trop de fenêtre d'où la liste 'list.lm' qui retourne les différentes analyses qui pourront être réutilisée si nécessaire pour les plot.

Bonjour,

je vous remercie pour votre réponse c'est très gentil à vous. J'ai testé votre proposition et les première commandes fonctionnent mais je suis bloquée pour la boucle.
Code:

1 2 3 4 5 6 7 8 9 setwd ("D:\\forums") fleurs<-read.csv2("fleurs.csv", na.strings = "NA",header = TRUE) noms <- names(fleurs) # donnees est le data.frame des donnees df.resultat <- data.frame( R2=rep( NA, length( noms) - 1), P=NA) rownames( df.resultat) <- noms[-1]
Jusque ici ça fonctionne

mon fichier de test ressemble à ça:

duree_croissance lieu annee plante temperature hum
1 105 lieu1 1990 iris 20 90
2 50 lieu2 1991 rose 10 20
3 108 lieu1 1990 iris 22 90
4 50 lieu2 1991 rose 8 20
5 20 <NA> 1992 tulipe NA 60
6 30 <NA> 1992 tulipe 2 NA

A partir de la boucle ça ne fonctionne pas, comme résultats j'obtient cela:

> df.resultat
R2 P
lieu NA NA
annee NA NA
plante NA NA
temperature NA NA
hum NA NA

Je dois sans doute faire quelque chose de faux, pouvez-vous me dire où je fais faux s'il vous plait? J'ai le message d'erreur suivant:
Code:

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 > for( noms in noms[-1]) { + form <- formula( paste( noms[1], nom, sep="~") + lm.prov <- lm( form, df.resultat) # na.omit par defaut Error: unexpected symbol in: " form <- formula( paste( noms[1], nom, sep="~") lm.prov" > list.lm[[nom]] <- lm.prov Error: object 'lm.prov' not found > sum.prov <- summary( lm.prov) Error in summary(lm.prov) : object 'lm.prov' not found > df.resultat[nom,"R2"] <- sum.prov$r.squared Error: object 'sum.prov' not found > # meme chose pour les autres valeurs > } Error: unexpected '}' in "}"
Merci par avance.

26/06/2017, 16h13
faubry

En tapant le code, j'ai oublié la seconde parenthèse fermante de la ligne :

Code:

form <- formula( paste( noms[1], nom, sep="~")

après l'argument 'sep', la parenthèse ferme les arguments du paste() et revoie une chaîne de caractères qui est l'argument de 'formula()'

Merci, je viens d'ajouter une parenthese mais j'ai maintenant un autre message d'erreur:

Citation:

Error in eval(expr, envir, enclos) : object 'duree_croissance' not found

le code est le suivant

Code:

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
setwd ("D:\\forums")
 
fleurs<-read.csv2("fleurs.csv", na.strings = "NA",header = TRUE)
 
 
noms <- names(fleurs) # donnees est le data.frame des donnees
 
df.resultat <- data.frame( R2=rep( NA, length( noms) - 1), P=NA)
 
rownames( df.resultat) <- noms[-1]
 
 
 
list.lm <- setNames(( vector( "list", length( noms) - 1)), noms[-1])
 
# On suppose des regressions de la premiere colonne sur les autres
 
 
                for( nom in noms[-1]) {
                  form <- formula( paste( noms[1], nom, sep="~"))
                                   lm.prov <- lm( form, df.resultat) # na.omit par defaut
                                   list.lm[[nom]] <- lm.prov
                                   sum.prov <- summary( lm.prov)
                                   df.resultat[nom,"R2"] <- sum.prov$r.squared
                                   # meme chose pour les autres valeurs
                }

26/06/2017, 18h49
faubry

Autre erreur de retranscription dans mon script : dans 'lm', il faut référencer les données initiales (ici fleurs) et non le résultat.

super ça fonctionne parfaitement, merci beaucoup!! :) :D

voici le resultat
Code:

1 2 3 4 5 6 7 > df.resultat A B R2ad R2 N P RMSE lieu -56.5000000 106.50000 0.9978885 0.9985923 4 0.0007040877 1.500000 annee -40.7500000 81193.75000 0.9314746 0.9451797 4 0.0011482294 9.813893 plante -56.5000000 106.50000 0.9870746 0.9922448 4 0.0009277088 4.262237 temperature 4.1765537 16.81073 0.9685926 0.9764444 4 0.0015453485 6.302706 hum 0.7565041 24.23577 0.3017793 0.4763344 4 0.1971184538 32.122120
juste une petite question sur cette boucle: j’extrais le nombre de lignes avec

Code:

df.resultat[nom,"N"] <- c(nrow(na.omit(fleurs)))

, je veux le nombre de lignes sans les na. Le soucis c'est qu'il me met 4 partout meme pour la colonne "annee" par exemple qui elle n'a pas de données manquantes et à donc un N= 6. Je pense qu'il prendre le nrow minimum sans les na de mon fichier fleurs, comment avoir le nombre de lignes sans les na pour chacune des colonnes?

Savez vous si dans la même boucle je peux ajouter l'effet de plusieurs variables. C à d que la boucle me calcule la corrélation de la colonne 1 avec les autres colonnes une par une comme le fait la boucle ci-dessous avec les régressions de la col 1 en fonction des autres, puis aussi d'ajouter par exemple une régression col 1 ~ col 2 + col 3 et ensuite col 1 ~ col 2 + col3 + col 4 et puis col 1 ~ col3+ col4 + col 2 afin de déterminer quels sont les facteurs qui influences le plus la croissance "col 1" ou les combinaisons de facteurs. comment faire un loop de régressions multiples automatiques et extraire ces données aussi? Ou peut être sous forme d'anova?

Merci pour tout et bonne soirée.
Tsuki.

26/06/2017, 22h03
faubry
Le nombre de lignes sans NA pour la colonne de nom 'nom' peut s'écritre

Code:

df.resultat[nom,"N"] <- sum( !is.na( fleurs[,nom]))

Si on fait

Code:

lm( y ~ x1 + x2, donnees)

l'analyse ne se fera que sur les lignes complètes, c'est-à-dire celle pour lesquelles il n'y a ni NA dans la colonne y, ni dans la colonne x1, ni dans la colonne x2. Donc, par exemple si x1 est la plante et x2, l'année, l'analyse ne prendra pas en compte les tulipes.

La recherche des prédicteurs les plus influents se fait en dehors de la boucle. Une solution est d'utiliser la fonction step en partant du modèle nul :
Code:

1 2 3 4 upper <- formula( paste0( " ~ ", paste( nom[-1], collapse=" + "))) nul <- formula( paste( nom[1], " ~ 1")) res <- step( lm( nul, fleurs), scope=list( upper=upper)) anova( res)
L'ordre de la formule donnera l'ordre des influences des prédicteurs et le résultat de anova, donnera ce qu'ajoute chacun des prédicteurs à l'ajustement par rapport aux précédents déjà introduits dans la formule. Cependant, étant donné que des lignes sont incomplètes, à chaque étape l'analyse risque de ne pas être fait sur le même jeu de données (cf. ma remarque précédente), ce qui rend cette estimation de l'influence des prédicteurs discutable.

Autres remarques :
- pourquoi utiliser un modèle sans interaction, ce qui revient à dire par exemple que le lieu ou l'humidité a la même influence sur toutes les plantes
- certains prédicteurs devraient être fortement corrélés entre eux, ce qu'il faut vérifier ; en effet, on peut penser que lieu et humidité (et température) sont très liés entre eux. L'utilisation dans une analyse de prédicteurs fortement corrélées fausse l'analyse.
27/06/2017, 08h48
Tsüki
merci, la commande pour le nombre de lignes sans NA fonctionne bien.

je reprends vos remarques une à une pour essayer de mieux comprendre:
- lm( y ~ x1 + x2, donnees) --> l'analyse ne se fera que sur les lignes complètes, c'est-à-dire celle pour lesquelles il n'y a ni NA dans la colonne y, ni dans la colonne x1, ni dans la colonne x2. Donc, par exemple si x1 est la plante et x2, l'année, l'analyse ne prendra pas en compte les tulipes.

donc pour une regression simple je peux correler la col 1 (duree de croissance = variable à expliquer / prédire) à la col 2 (humidité) même si je n'ai pas le même nombre de valeurs à cause des NA mais pour une regression multiple ce n'est pas possible c'est bien ça? N'y a t-il pas moyen de "corriger" le nombre de valeurs différentes de chaque colonne?

--> La recherche des prédicteurs les plus influents se fait en dehors de la boucle. Une solution est d'utiliser la fonction step en partant du modèle nul.

L'ordre de la formule donnera l'ordre des influences des prédicteurs et le résultat de anova, donnera ce qu'ajoute chacun des prédicteurs à l'ajustement par rapport aux précédents déjà introduits dans la formule. :

Si je veux trouver les prédicteurs les plus influents (parmis toutes les colonnes) sur ma croissance (col1) je dois appliquer la formule ci-dessous qui va me dire l'influence de chacune des collones sur la colonne 1 (grâce à la fonction step et lm?) dans l'ordre; et ensuite l'anova donne ce qu'ajoute chacun des prédicteurs (colonnes prises une à une) en tenant compte de l'ordre donné par la fonction step, c'est bien ça?
Code:

1 2 3 4 upper <- formula( paste0( " ~ ", paste( nom[-1], collapse=" + "))) nul <- formula( paste( nom[1], " ~ 1")) res <- step( lm( nul, fleurs), scope=list( upper=upper)) anova( res)
--> Cependant, étant donné que des lignes sont incomplètes, à chaque étape l'analyse risque de ne pas être fait sur le même jeu de données (cf. ma remarque précédente), ce qui rend cette estimation de l'influence des prédicteurs discutable. Ici il n'y a pas de solution pour "corriger" le fait que le N soit différent d'une colonne à l'autre?

--> pourquoi utiliser un modèle sans interaction, ce qui revient à dire par exemple que le lieu ou l'humidité a la même influence sur toutes les plantes
que voulez vous dire par un modèle sans interraction? J'ai du mal à orienter mon modèle. Pour mon exemple, le but est de déterminer l'influence des différents paramètres sur la durée de la croissance. Je sais déjà que certains paramètres sont très influents, le plus influents est le type de plante par exemple, puis ensuite je pense que ça sera la température. La première étape serait donc de détermine quels paramètres influencent significativement la durée de la croissance et ensuite de les classer non?

En fait dans une seconde étape j'aimerais (après avoir selectionné les paramètres influents) pouvoir prédire pour une année à venir la durée de la croissance pour un type de plante (tulipe) sur base des paramètres déterminés les plus influents (ex: température, humidité, vent). Pensez-vous que cela est possible? Sachant que dans mon vrai fichier j'ai plus de 200 colonnes (paramètres potentiellements influents avec un N qui varie) et j'ai plus de 3000 lignes (lieux, types de plantes, années d'essais différents). Par contre j'ai la durée de la croissance pour les 3000 lignes.

--> certains prédicteurs devraient être fortement corrélés entre eux, ce qu'il faut vérifier ; en effet, on peut penser que lieu et humidité (et température) sont très liés entre eux. L'utilisation dans une analyse de prédicteurs fortement corrélées fausse l'analyse.

je peux vérifier cela avec ma boucle de regressions simples?

Un grand merci pour tout, étant débutante en statistiques vraiments appliquées, vos connaissances sont d'une grande aide.
27/06/2017, 11h11
faubry

Citation:

donc pour une regression simple je peux correler la col 1 (duree de croissance = variable à expliquer / prédire) à la col 2 (humidité) même si je n'ai pas le même nombre de valeurs à cause des NA mais pour une regression multiple ce n'est pas possible c'est bien ça? N'y a t-il pas moyen de "corriger" le nombre de valeurs différentes de chaque colonne?

Non, sauf à utiliser une méthode d'imputation des valeurs manquantes. Mais cela est loin de résoudre tous les problèmes et même peut s'avérer plus catastrophique que de ne pas prendre ne compte les données manquantes. Cependant, le problème apparaît si les données sont vraiment très déséquilibrées, c'est-à-dire le nombre de sujets par cellule est trop dispersé. S'il est à peu près constant, on peut dans un premier temps négliger le biais.

Citation:

Si je veux trouver les prédicteurs les plus influents (parmis toutes les colonnes) sur ma croissance (col1) je dois appliquer la formule ci-dessous qui va me dire l'influence de chacune des collones sur la colonne 1 (grâce à la fonction step et lm?) dans l'ordre; et ensuite l'anova donne ce qu'ajoute chacun des prédicteurs (colonnes prises une à une) en tenant compte de l'ordre donné par la fonction step, c'est bien ça?

Oui, mais pour éviter les problèmes, je n'utiliserai que le colonnes qui, dans la première étude, donnent un p inférieur ou égal à 0.05 (ou 0.1 pour me donner une marge de sécurité).

Citation:

que voulez vous dire par un modèle sans interraction?

Comme je l'ai dit, le modèle fait l'hypothèse que le lieu à la même influence sur la durée de croissance quelque soit la plante, ainsi que la température... Or le nature du sol, l'hygrométrie, la température... devraient avoir des influences différentes selon les plantes. Si c'est le cas, les interactions tiennent compte de ces phénomènes.

Citation:

certains prédicteurs devraient être fortement corrélés entre eux, ce qu'il faut vérifier ; en effet, on peut penser que lieu et humidité (et température) sont très liés entre eux. L'utilisation dans une analyse de prédicteurs fortement corrélées fausse l'analyse.

je peux vérifier cela avec ma boucle de regressions simples?

Non. C'est plus compliqué que cela. Il y a plusieurs solutions. La technique que j'utilise personnellement et que je recommande aux étudiants est l'utilisation en première intention des statistiques exploratoires sur les variables indépendantes (en ajoutant comme variable supplémentaires, donc simplement illustrative, la variable dépendante). Ici comme il y a des facteurs et des variables quantitatives, l'analyse à facteurs multiples peut-être utilisé avec trois groupes, deux groupes actifs, l'un des données quantitatives, l'autre des facteurs, et un groupe supplémentaire, la variable dépendante. Voir par exemple, la fonction MFA du package FactoMineR.
Une manière aussi de savoir s'il y a des colinéarités entre prédicteurs est de regarder le VIF (variance inflation factor) de l'analyse linéaire (fonction vif de ). Des valeurs supérieures à 3 ou 5 (selon les auteurs) pointent l'existence de colinéarités. Cependant, le IF ne dit pas lesquelles.
27/06/2017, 19h29
Tsüki

merci pour ces informations, je vais chercher plus en détail selon vos propositions car cela dépasse mes capacités pour le moment, je viendrai reposer des questions quand j'aurai un peu mieux compris.

Je peux déjà terminé la première étape de boucle et régressions et sélectionner les paramètres les mieux corrélés comme vous le suggérez.
Alors effectivement dans ma boucle ou je corrèle la colonne 1 avec chacune des 200 colonnes prises séparément, j'obtiens ce message d'erreur:

Citation:

Error in lm.fit(x, y, offset = offset, singular.ok = singular.ok, ...) :
aucun cas ne contient autre chose que des valeurs manquantes (NA)

et quand je regarde les résultats j'ai beaucoup de NA, pour certains facteurs où j'ai normalement des valeurs, il me sort des "NA" pour le R² et même pour le nombre de valeurs N, est-ce que cela peut être dû au fait que j'ai trop de NA et que certaines régressions ne se font pas?
28/06/2017, 09h12
faubry
Pour se faire une idée du nombre de NA par colonne :

Code:

colMeans( is.na( fleurs))

Ensuite pour les autres facteurs, il peut n'exister qu'un niveau ou certains niveaux qui n'ont qu'un individu. Pour vérifier :
Code:

1 2 fac <- sapply( fleurs, is.factor) lapply( fleurs[,fac], table)

super cette fonction, merci!

y a t il un moyen d'écrire le resultats de lapply?
Code:

1 2 fac <- sapply( fleurs, is.factor) lapply( fleurs[,fac], table)
j'ai tenté de stocker la fonction et de l'écrire
Code:

1 2 3 4 5 stock<- lapply( fleurs[,fac], table) write.table(stock...) ou d'autres "write" mais le résultat est le suivant: Error in (function (..., row.names = NULL, check.rows = FALSE, check.names = TRUE, : arguments imply differing number of rows: 179, 5
je pense que cela est dû au fait que ça ne ressemble pas à une matrice