вытащить р-значения и г-квадрат из линейной регрессии

179

Как вы извлекаете значение p (для значимости коэффициента единственной объясняющей переменной, отличной от нуля) и значение R-квадрата из простой модели линейной регрессии? Например...

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

Я знаю, что summary(fit) отображает значение p и значение R в квадрате, но я хочу иметь возможность вставить их в другие переменные.

grautur
источник
Он отображает значения только в том случае, если вы не назначаете вывод для объекта (например, r <- summary(lm(rnorm(10)~runif(10)))ничего не показывает).
Джошуа Ульрих

Ответы:

157

r-квадрат : вы можете вернуть значение r-квадрат непосредственно из итогового объекта summary(fit)$r.squared. Смотрите names(summary(fit))список всех предметов, которые вы можете извлечь непосредственно.

Значение p модели. Если вы хотите получить значение p общей модели регрессии, в этом сообщении в блоге описывается функция, возвращающая значение 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

В случае простой регрессии с одним предиктором, p-значение модели и p-значение для коэффициента будут одинаковыми.

P-значения коэффициентов: если у вас есть более одного предиктора, то приведенное выше вернет p-значение модели, и p-значение для коэффициентов можно извлечь с помощью:

summary(fit)$coefficients[,4]  

В качестве альтернативы, вы можете получить p-значение коэффициентов из anova(fit)объекта аналогично итоговому объекту выше.

гнаться
источник
13
Это немного лучше использовать inherits, чем classнапрямую. А может ты хочешь unname(pf(f[1],f[2],f[3],lower.tail=F))?
Хэдли
150

Обратите внимание, что summary(fit)генерируется объект со всей необходимой вам информацией. В нем хранятся векторы бета, se, t и p. Получите p-значения, выбрав 4-й столбец матрицы коэффициентов (хранится в итоговом объекте):

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

Попробуйте str(summary(fit))просмотреть всю информацию, которая содержится в этом объекте.

Изменить: я неправильно прочитал ответ Чейза, который в основном говорит вам, как добраться до того, что я здесь даю.

Винсент
источник
11
Примечание: это единственный метод, который дает вам легкий доступ к p-значению перехвата, а также к другим предикторам. Безусловно лучшее из вышеперечисленного.
Даниэль Иган
2
Это правильный ответ. Ответ с самым высоким рейтингом НЕ работал для меня.
Крис
8
ЕСЛИ ВЫ ХОТИТЕ ЛЕГКИЙ ДОСТУП К P-VALUE, ИСПОЛЬЗУЙТЕ ЭТОТ ОТВЕТ. Зачем вам заниматься написанием многострочных функций или созданием новых объектов (например, выводов anova), когда вам просто нужно немного сложнее найти p-значение в самом итоговом выводе. Чтобы изолировать отдельное значение p непосредственно, вы должны добавить номер строки к ответу Винсента: например, summary(fit)$coefficients[1,4] для их инцепции
лесничий
2
Примечание: этот метод работает для моделей, созданных с использованием, lm()но не работает для gls()моделей.
лесовод
3
Ответ Чейза возвращает p-значение модели, этот ответ возвращает p-значение коэффициентов. В случае простой регрессии они одинаковы, но в случае модели с несколькими предикторами они не одинаковы. Таким образом, оба ответа полезны в зависимости от того, что вы хотите извлечь.
Джером Энглим
44

Вы можете увидеть структуру объекта, возвращенного summary()путем вызова str(summary(fit)). Каждая часть может быть доступна с помощью $. Значение p для статистики F легче получить из возвращаемого объекта anova.

Вкратце, вы можете сделать это:

rSquared <- summary(fit)$r.squared
pVal <- anova(fit)$'Pr(>F)'[1]
jberg
источник
10
это работает только для одномерных регрессий, где
реальная
23

Хотя оба приведенных выше ответа хороши, процедура извлечения частей объектов носит более общий характер.

Во многих случаях функции возвращают списки, и к отдельным компонентам можно получить доступ, используя str()которые будут печатать компоненты вместе с их именами. Вы можете получить доступ к ним с помощью $ оператора, то есть myobject$componentname.

В случае объектов lm есть несколько предопределенных методов, которые можно использовать, например coef(), resid()и summary()т. Д., Но вам не всегда повезет.

richiemorrisroe
источник
23

Я сталкивался с этим вопросом, изучая предложенные решения для подобной проблемы; Я предполагаю, что для дальнейшего использования может быть целесообразно обновить доступный список ответов решением, использующим broomпакет.

Образец кода

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

Полученные результаты

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

Примечания стороны

Я считаю, что glanceфункция полезна, так как она аккуратно суммирует ключевые значения. Результаты сохраняются в виде, data.frameчто облегчает дальнейшие манипуляции:

>> class(glance(fit))
[1] "data.frame"
Konrad
источник
Это отличный ответ!
Андрей Бреза
9

Расширение @Vincent «s ответ :

Для lm()сгенерированных моделей:

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

Для gls()сгенерированных моделей:

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

Чтобы изолировать отдельное p-значение, вы должны добавить номер кода в код:

Например, чтобы получить доступ к p-значению перехвата в обеих сводках модели:

summary(fit)$coefficients[1,4]
summary(fit)$tTable[1,4]  
  • Обратите внимание, что вы можете заменить номер столбца на имя столбца в каждом из приведенных выше случаев:

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

Если вы все еще не знаете, как получить доступ к значению из сводной таблицы, используйте ее str()для определения структуры сводной таблицы:

str(summary(fit))
theforestecologist
источник
7

Это самый простой способ получить значения p:

coef(summary(modelname))[, "Pr(>|t|)"]
RTrain3k
источник
1
Я попробовал этот метод, но он потерпит неудачу, если линейная модель содержит какие-либо термины NA
j_v_wow_d
5

Я использовал эту функцию LMP довольно много раз.

И в один прекрасный момент я решил добавить новые функции для улучшения анализа данных. Я не эксперт в R или статистике, но люди обычно смотрят на различную информацию о линейной регрессии:

  • р-значение
  • а и б
  • и, конечно, аспект распределения точек

Давайте иметь пример. У вас есть здесь

Вот воспроизводимый пример с разными переменными:

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)

Конечно, есть более быстрое решение, чем эта функция, но оно работает.

Дориан Грв
источник
2

Для конечного p-значения, отображаемого в конце summary(), функция использует pf()для вычисления summary(fit)$fstatisticзначений.

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

Источник: [1] , [2]

Saftever
источник
1
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
Jojo
источник
1
Не забудьте дать объяснение, даже если вкратце, почему этот код работает?
Aribeiro
как это улучшает существующие ответы (и в частности принятый ответ)?
Бен Болкер
0

Другой вариант - использовать функцию cor.test вместо 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
dalloliogm
источник
0

Использование:

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

где numчисло, которое обозначает строку матрицы коэффициентов. Это будет зависеть от того, сколько функций у вас есть в вашей модели и для какой вы хотите получить значение p. Например, если у вас есть только одна переменная, для пересечения будет одно p-значение, которое будет [1,4], и следующее для вашей фактической переменной, которое будет [2,4]. Таким образом, ваш numбудет 2.

Tirtha
источник