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

Как вы вытаскиваете 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-квадрата, но я хочу иметь возможность привязать их к другим переменным.

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

Модель p-value: если вы хотите получить 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) аналогично суммарному объекту выше.

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

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

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

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

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

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

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

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

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

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

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

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

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

Для gls() созданных 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)) 

Я рассматриваю этот вопрос, исследуя предлагаемые решения для аналогичной проблемы; Я полагаю, что для будущей справки может оказаться целесообразным обновить ansible список ответов с помощью решения, использующего пакет 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" 

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

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

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

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

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

Приведем пример. У вас здесь

Здесь воспроизводимый пример с разными переменными:

 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) 

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

 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 

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

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

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

Другой вариант – использовать функцию 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 
  • Получить «встроенные nul (s), найденные во вводе» при чтении csv, используя read.csv ()
  • конвертировать письма в номера
  • Есть ли способ сделать R звуковой сигнал / воспроизвести звук в конце скрипта?
  • Получение OVER QUERY LIMIT после одного запроса с геокодом
  • Как назначить значения динамическим именам переменных
  • Глобальные и локальные переменные в R
  • R создать идентификатор внутри группы
  • Использование R для отображения всех файлов с указанным расширением
  • Преобразовать значения в столбце в имена строк в существующем кадре данных в R
  • Каковы форматы «стандартной четкой даты»?
  • Открытие всех файлов в папке и применение функции
  • Давайте будем гением компьютера.