entre - regresion lineal

extraer valores p y r-cuadrado de una regresión lineal

Question

¿Cómo extrae el valor p (para la significación del coeficiente de la variable explicativa única que no es cero) y el valor R-cuadrado de un modelo de regresión lineal simple? Por ejemplo...

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

Sé que el summary(fit) muestra el valor p y el valor R-cuadrado, pero quiero poder incluirlos en otras variables.

Answer 1

Abordo esta pregunta mientras exploro soluciones sugeridas para un problema similar; Supongo que para referencia futura puede valer la pena actualizar la lista de respuestas disponible con una solución que utiliza el paquete de broom .

Código de muestra

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

Resultados

>> 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

Notas laterales

Encuentro la función de glance útil ya que resume perfectamente los valores útiles. Como beneficio adicional, los resultados se almacenan como un data.frame que facilita la manipulación adicional:

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

Answer 2

Esta es la forma más fácil de extraer los valores p:

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

Answer 3

Extensión de la respuesta de @Vincent:

Para lm() modelos generados:

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

Para modelos generados por gls() :

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

Para aislar un valor p individual en sí mismo, debe agregar un número de fila al código:

Por ejemplo, para acceder al valor p de la intersección en ambos resúmenes de modelos:

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

• Tenga en cuenta que puede reemplazar el número de columna con el nombre de columna en cada una de las instancias anteriores:

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

Si todavía no está seguro de cómo acceder a un valor desde la tabla de resumen, use str() para descubrir la estructura de la tabla de resumen:

str(summary(fit))

Answer 4

Otra opción es usar la función cor.test, en lugar 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

Answer 5

Puede ver la estructura del objeto devuelto por summary() llamando a str(summary(fit)) . Se puede acceder a cada pieza usando $ . El valor p para la estadística F se obtuvo más fácilmente del objeto devuelto por anova .

De manera concisa, puedes hacer esto:

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

Answer 6

Si bien ambas respuestas anteriores son buenas, el procedimiento para extraer partes de objetos es más general.

En muchos casos, las funciones devuelven listas, y se puede acceder a los componentes individuales usando str() que imprimirá los componentes junto con sus nombres. A continuación, puede acceder a ellos utilizando el operador $, es decir, myobject$componentname .

En el caso de los objetos lm, hay una serie de métodos predefinidos que se pueden usar, como coef() , resid() , summary() , etc., pero no siempre será tan afortunado.

Answer 7

Tenga en cuenta que el summary(fit) genera un objeto con toda la información que necesita. Los vectores beta, se, t y p se almacenan en él. Obtenga los valores de p seleccionando la cuarta columna de la matriz de coeficientes (almacenada en el objeto de resumen):

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

Pruebe str(summary(fit)) para ver toda la información que contiene este objeto.

Editar: He leído mal la respuesta de Chase, que básicamente te dice cómo llegar a lo que doy aquí.

Answer 8

Usé esta función lmp bastantes veces.

Y en un punto decidí agregar nuevas características para mejorar el análisis de datos. No soy experto en R ni en estadísticas, pero las personas generalmente miran información diferente de una regresión lineal:

• valor p
• a y B
• r²
• y, por supuesto, el aspecto de la distribución de puntos

Vamos a tener un ejemplo. Usted tiene aquí

Aquí un ejemplo reproducible con diferentes 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)

Sin duda, hay una solución más rápida que esta función, pero funciona.

Answer 9

Utilizar:

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

donde num es un número que denota la fila de la matriz de coeficientes. Dependerá de la cantidad de características que tenga en su modelo y de la que desee extraer el valor p. Por ejemplo, si solo tiene una variable, habrá un valor p para la intersección que será [1,4] y el siguiente para su variable real que será [2,4]. Entonces tu num será 2.

Answer 10

r-cuadrado : puede devolver el valor r-cuadrado directamente desde el resumen del objeto summary(fit)$r.squared . Vea los names(summary(fit)) para obtener una lista de todos los elementos que puede extraer directamente.

Valor p del modelo: si desea obtener el valor p del modelo de regresión general, esta publicación de blog describe una función para devolver el valor 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

En el caso de una regresión simple con un predictor, el valor p del modelo y el valor p para el coeficiente serán los mismos.

Valores p del coeficiente: si tiene más de un predictor, lo anterior devolverá el valor p del modelo, y el valor p para los coeficientes se puede extraer usando:

summary(fit)$coefficients[,4]  

Alternativamente, puede tomar el valor p de los coeficientes del objeto anova(fit) de forma similar al objeto de resumen anterior.

Answer 11

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