Comment extraire la valeur p (pour la signification du coefficient de la variable explicative unique étant non nulle) et la valeur R au carré d'un modèle de régression linéaire simple? Par exemple...

x = cumsum(c(0, runif(100, -1, +1)))
y = cumsum(c(0, runif(100, -1, +1)))
fit = lm(y ~ x)
summary(fit)

Je sais que cela summary(fit) affiche la valeur p et la valeur R au carré, mais je veux pouvoir les coller dans d'autres variables.

r-squared : vous pouvez renvoyer la valeur r-squared directement à partir de l'objet de résumé summary(fit)$r.squared. Consultez names(summary(fit))la liste de tous les éléments que vous pouvez extraire directement.

Valeur p du modèle: Si vous souhaitez obtenir la valeur p du modèle de régression global, cet article de blog décrit une fonction pour renvoyer la valeur p:

lmp <- function (modelobject) {
    if (class(modelobject) != "lm") stop("Not an object of class 'lm' ")
    f <- summary(modelobject)$fstatistic
    p <- pf(f[1],f[2],f[3],lower.tail=F)
    attributes(p) <- NULL
    return(p)
}

> lmp(fit)
[1] 1.622665e-05

Dans le cas d'une régression simple avec un seul prédicteur, la p-value du modèle et la p-value du coefficient seront les mêmes.

Valeurs p du coefficient: Si vous avez plus d'un prédicteur, alors ce qui précède renverra la valeur p du modèle, et la valeur p des coefficients peut être extraite en utilisant:

summary(fit)$coefficients[,4]  

Vous pouvez également saisir la valeur p des coefficients de l' anova(fit)objet de la même manière que l'objet récapitulatif ci-dessus.

Remarquez que summary(fit)génère un objet avec toutes les informations dont vous avez besoin. Les vecteurs bêta, se, t et p y sont stockés. Obtenez les p-values ​​en sélectionnant la 4ème colonne de la matrice des coefficients (stockée dans l'objet de synthèse):

summary(fit)$coefficients[,4] 
summary(fit)$r.squared

Essayez str(summary(fit))de voir toutes les informations que cet objet contient.

Edit: j'avais mal lu la réponse de Chase qui vous dit essentiellement comment arriver à ce que je donne ici.

Vous pouvez voir la structure de l'objet renvoyé par summary()en appelant str(summary(fit)). Chaque pièce est accessible en utilisant $. La valeur p de la statistique F est plus facilement obtenue à partir de l'objet renvoyé par anova.

En résumé, vous pouvez le faire:

rSquared <- summary(fit)$r.squared
pVal <- anova(fit)$'Pr(>F)'[1]

Bien que les deux réponses ci-dessus soient bonnes, la procédure d'extraction de parties d'objets est plus générale.

Dans de nombreux cas, les fonctions renvoient des listes et les composants individuels sont accessibles en utilisant str()ce qui imprimera les composants avec leurs noms. Vous pouvez ensuite y accéder en utilisant l'opérateur $, c'est-à-dire myobject$componentname.

Dans le cas des objets lm, il y a un certain nombre de méthodes prédéfinies , on peut utiliser par exemple coef(), resid(), summary()etc, mais vous ne serez pas toujours aussi chanceux.

Je suis tombé sur cette question en explorant des solutions suggérées pour un problème similaire; Je suppose que pour référence future, il peut être utile de mettre à jour la liste de réponses disponibles avec une solution utilisant le broompackage.

Exemple de code

x = cumsum(c(0, runif(100, -1, +1)))
y = cumsum(c(0, runif(100, -1, +1)))
fit = lm(y ~ x)
require(broom)
glance(fit)

Résultats

>> glance(fit)
  r.squared adj.r.squared    sigma statistic    p.value df    logLik      AIC      BIC deviance df.residual
1 0.5442762     0.5396729 1.502943  118.2368 1.3719e-18  2 -183.4527 372.9055 380.7508 223.6251          99

Notes annexes

Je trouve que la glancefonction est utile car elle résume parfaitement les valeurs clés. Les résultats sont stockés sous forme de, data.framece qui facilite la manipulation ultérieure:

>> class(glance(fit))
[1] "data.frame"

Extension de l » @Vincent réponse :

Pour les lm()modèles générés:

summary(fit)$coefficients[,4]   ##P-values 
summary(fit)$r.squared          ##R squared values

Pour les gls()modèles générés:

summary(fit)$tTable[,4]         ##P-values
##R-squared values are not generated b/c gls uses max-likelihood not Sums of Squares

Pour isoler une valeur de p individuelle elle-même, ajoutez un numéro de ligne au code:

Par exemple, pour accéder à la valeur p de l'intersection dans les deux résumés du modèle:

summary(fit)$coefficients[1,4]
summary(fit)$tTable[1,4]  

• Notez que vous pouvez remplacer le numéro de colonne par le nom de la colonne dans chacune des instances ci-dessus:

  summary(fit)$coefficients[1,"Pr(>|t|)"]  ##lm 
  summary(fit)$tTable[1,"p-value"]         ##gls 

Si vous ne savez toujours pas comment accéder à un formulaire de valeur, le tableau récapitulatif utilise str()pour déterminer la structure du tableau récapitulatif:

str(summary(fit))

C'est le moyen le plus simple d'extraire les valeurs p:

coef(summary(modelname))[, "Pr(>|t|)"]

J'ai utilisé cette fonction lmp plusieurs fois.

Et à un moment donné, j'ai décidé d'ajouter de nouvelles fonctionnalités pour améliorer l'analyse des données. Je ne suis pas expert en R ou en statistiques mais les gens recherchent généralement des informations différentes d'une régression linéaire:

• valeur p
• a et b
• r²
• et bien sûr l'aspect de la distribution des points

Prenons un exemple. Vous avez ici

Voici un exemple reproductible avec différentes variables:

Ex<-structure(list(X1 = c(-36.8598, -37.1726, -36.4343, -36.8644, 
-37.0599, -34.8818, -31.9907, -37.8304, -34.3367, -31.2984, -33.5731
), X2 = c(64.26, 63.085, 66.36, 61.08, 61.57, 65.04, 72.69, 63.83, 
67.555, 76.06, 68.61), Y1 = c(493.81544, 493.81544, 494.54173, 
494.61364, 494.61381, 494.38717, 494.64122, 493.73265, 494.04246, 
494.92989, 494.98384), Y2 = c(489.704166, 489.704166, 490.710962, 
490.653212, 490.710612, 489.822928, 488.160904, 489.747776, 490.600579, 
488.946738, 490.398958), Y3 = c(-19L, -19L, -19L, -23L, -30L, 
-43L, -43L, -2L, -58L, -47L, -61L)), .Names = c("X1", "X2", "Y1", 
"Y2", "Y3"), row.names = c(NA, 11L), class = "data.frame")


library(reshape2)
library(ggplot2)
Ex2<-melt(Ex,id=c("X1","X2"))
colnames(Ex2)[3:4]<-c("Y","Yvalue")
Ex3<-melt(Ex2,id=c("Y","Yvalue"))
colnames(Ex3)[3:4]<-c("X","Xvalue")

ggplot(Ex3,aes(Xvalue,Yvalue))+
          geom_smooth(method="lm",alpha=0.2,size=1,color="grey")+
          geom_point(size=2)+
          facet_grid(Y~X,scales='free')


#Use the lmp function

lmp <- function (modelobject) {
  if (class(modelobject) != "lm") stop("Not an object of class 'lm' ")
  f <- summary(modelobject)$fstatistic
    p <- pf(f[1],f[2],f[3],lower.tail=F)
    attributes(p) <- NULL
    return(p)
    }

# create function to extract different informations from lm

lmtable<-function (var1,var2,data,signi=NULL){
  #var1= y data : colnames of data as.character, so "Y1" or c("Y1","Y2") for example
  #var2= x data : colnames of data as.character, so "X1" or c("X1","X2") for example
  #data= data in dataframe, variables in columns
  # if signi TRUE, round p-value with 2 digits and add *** if <0.001, ** if < 0.01, * if < 0.05.

  if (class(data) != "data.frame") stop("Not an object of class 'data.frame' ")
  Tabtemp<-data.frame(matrix(NA,ncol=6,nrow=length(var1)*length(var2)))
  for (i in 1:length(var2))
       {
  Tabtemp[((length(var1)*i)-(length(var1)-1)):(length(var1)*i),1]<-var1
  Tabtemp[((length(var1)*i)-(length(var1)-1)):(length(var1)*i),2]<-var2[i]
  colnames(Tabtemp)<-c("Var.y","Var.x","p-value","a","b","r^2")

  for (n in 1:length(var1))
  {
  Tabtemp[(((length(var1)*i)-(length(var1)-1))+n-1),3]<-lmp(lm(data[,var1[n]]~data[,var2[i]],data))

  Tabtemp[(((length(var1)*i)-(length(var1)-1))+n-1),4]<-coef(lm(data[,var1[n]]~data[,var2[i]],data))[1]

  Tabtemp[(((length(var1)*i)-(length(var1)-1))+n-1),5]<-coef(lm(data[,var1[n]]~data[,var2[i]],data))[2]

  Tabtemp[(((length(var1)*i)-(length(var1)-1))+n-1),6]<-summary(lm(data[,var1[n]]~data[,var2[i]],data))$r.squared
  }
  }

  signi2<-data.frame(matrix(NA,ncol=3,nrow=nrow(Tabtemp)))
  signi2[,1]<-ifelse(Tabtemp[,3]<0.001,paste0("***"),ifelse(Tabtemp[,3]<0.01,paste0("**"),ifelse(Tabtemp[,3]<0.05,paste0("*"),paste0(""))))
  signi2[,2]<-round(Tabtemp[,3],2)
  signi2[,3]<-paste0(format(signi2[,2],digits=2),signi2[,1])

  for (l in 1:nrow(Tabtemp))
    {
  Tabtemp$"p-value"[l]<-ifelse(is.null(signi),
         Tabtemp$"p-value"[l],
         ifelse(isTRUE(signi),
                paste0(signi2[,3][l]),
                Tabtemp$"p-value"[l]))
  }

   Tabtemp
}

# ------- EXAMPLES ------

lmtable("Y1","X1",Ex)
lmtable(c("Y1","Y2","Y3"),c("X1","X2"),Ex)
lmtable(c("Y1","Y2","Y3"),c("X1","X2"),Ex,signi=TRUE)

Il existe certes une solution plus rapide que cette fonction mais cela fonctionne.

Pour la valeur p finale affichée à la fin de summary(), la fonction utilise pf()pour calculer à partir des summary(fit)$fstatisticvaleurs.

fstat <- summary(fit)$fstatistic
pf(fstat[1], fstat[2], fstat[3], lower.tail=FALSE)

Source: [1] , [2]

x = cumsum(c(0, runif(100, -1, +1)))
y = cumsum(c(0, runif(100, -1, +1)))
fit = lm(y ~ x)
> names(summary(fit))
[1] "call"          "terms"        
 [3] "residuals"     "coefficients" 
 [5] "aliased"       "sigma"        
 [7] "df"            "r.squared"    
 [9] "adj.r.squared" "fstatistic"   
[11] "cov.unscaled" 
    summary(fit)$r.squared

Une autre option consiste à utiliser la fonction cor.test, au lieu de lm:

> x <- c(44.4, 45.9, 41.9, 53.3, 44.7, 44.1, 50.7, 45.2, 60.1)
> y <- c( 2.6,  3.1,  2.5,  5.0,  3.6,  4.0,  5.2,  2.8,  3.8)

> mycor = cor.test(x,y)
> mylm = lm(x~y)

# r and rsquared:
> cor.test(x,y)$estimate ** 2
      cor 
0.3262484 
> summary(lm(x~y))$r.squared
[1] 0.3262484

# P.value 

> lmp(lm(x~y))  # Using the lmp function defined in Chase's answer
[1] 0.1081731
> cor.test(x,y)$p.value
[1] 0.1081731

Utilisation:

(summary(fit))$coefficients[***num***,4]

où numest un nombre qui désigne la ligne de la matrice des coefficients. Cela dépendra du nombre de fonctionnalités que vous avez dans votre modèle et de celle pour laquelle vous souhaitez extraire la valeur p. Par exemple, si vous n'avez qu'une seule variable, il y aura une valeur p pour l'interception qui sera [1,4] et la suivante pour votre variable réelle qui sera [2,4]. Donc, vous numserez 2.