Bagaimana Anda mengeluarkan nilai-p (untuk signifikansi koefisien variabel penjelas tunggal yang tidak nol) dan nilai R-kuadrat dari model regresi linier sederhana? Sebagai contoh...

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

Saya tahu bahwa summary(fit) menampilkan nilai p-nilai dan nilai R-kuadrat, tapi saya ingin dapat memasukkan ini ke variabel lain.

r-squared : Anda dapat mengembalikan nilai r-squared langsung dari objek ringkasan summary(fit)$r.squared. Lihat names(summary(fit))daftar semua barang yang bisa Anda ekstrak langsung.

Nilai p Model: Jika Anda ingin mendapatkan nilai p dari model regresi keseluruhan, posting blog ini menguraikan fungsi untuk mengembalikan nilai 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

Dalam kasus regresi sederhana dengan satu prediktor, model p-value dan p-value untuk koefisien akan sama.

Nilai p koefisien: Jika Anda memiliki lebih dari satu prediktor, maka nilai di atas akan mengembalikan nilai p model, dan nilai p untuk koefisien dapat diekstraksi menggunakan:

summary(fit)$coefficients[,4]  

Atau, Anda bisa mengambil nilai p dari koefisien dari anova(fit)objek dengan cara yang mirip dengan objek ringkasan di atas.

Perhatikan bahwa summary(fit)menghasilkan objek dengan semua informasi yang Anda butuhkan. Vektor beta, se, t dan p disimpan di dalamnya. Dapatkan nilai-p dengan memilih kolom ke-4 dari matriks koefisien (disimpan dalam objek ringkasan):

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

Cobalah str(summary(fit))untuk melihat semua info yang berisi objek ini.

Sunting: Saya salah membaca jawaban Chase yang pada dasarnya memberi tahu Anda cara mendapatkan apa yang saya berikan di sini.

Anda dapat melihat struktur objek yang dikembalikan summary()dengan menelepon str(summary(fit)). Setiap bagian dapat diakses menggunakan $. Nilai p untuk statistik F lebih mudah didapat dari objek yang dikembalikan oleh anova.

Secara ringkas, Anda dapat melakukan ini:

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

Meskipun kedua jawaban di atas baik, prosedur untuk mengekstraksi bagian-bagian objek lebih umum.

Dalam banyak kasus, daftar fungsi kembali, dan masing-masing komponen dapat diakses menggunakan str()yang akan mencetak komponen beserta namanya. Anda kemudian dapat mengaksesnya menggunakan $ operator, yaitu myobject$componentname.

Dalam kasus obyek lm, ada sejumlah metode yang telah ditetapkan yang bisa digunakan seperti coef(), resid(), summary()dll, tetapi Anda tidak akan selalu begitu beruntung.

Saya menemukan pertanyaan ini sambil mengeksplorasi solusi yang disarankan untuk masalah yang sama; Saya berasumsi bahwa untuk referensi di masa mendatang mungkin ada gunanya memperbarui daftar jawaban yang tersedia dengan solusi memanfaatkan broompaket.

Kode sampel

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

Hasil

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

Catatan samping

Saya menemukan glancefungsi ini berguna karena meringkas nilai-nilai kunci dengan rapi. Hasil disimpan sebagai data.frameyang membuat manipulasi lebih mudah:

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

Perpanjangan jawaban @Vincent :

Untuk lm()model yang dihasilkan:

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

Untuk gls()model yang dihasilkan:

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

Untuk mengisolasi nilai p individu itu sendiri, Anda akan menambahkan nomor baris ke kode:

Misalnya untuk mengakses nilai p intersep di kedua ringkasan model:

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

• Catatan, Anda dapat mengganti nomor kolom dengan nama kolom di masing-masing contoh di atas:

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

Jika Anda masih tidak yakin tentang cara mengakses nilai dari tabel ringkasan gunakan str()untuk mengetahui struktur tabel ringkasan:

str(summary(fit))

Ini adalah cara termudah untuk menarik nilai-p:

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

Saya menggunakan fungsi lmp ini cukup banyak.

Dan pada satu titik saya memutuskan untuk menambahkan fitur baru untuk meningkatkan analisis data. Saya tidak ahli dalam R atau statistik tetapi orang biasanya melihat informasi yang berbeda dari regresi linier:

• nilai p
• a dan b
• r²
• dan tentu saja aspek dari distribusi titik

Mari kita punya contoh. Anda punya di sini

Di sini contoh yang direproduksi dengan variabel yang berbeda:

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)

Tentu ada solusi yang lebih cepat daripada fungsi ini tetapi berfungsi.

Untuk nilai-p akhir yang ditampilkan di akhir summary(), fungsi digunakan pf()untuk menghitung dari summary(fit)$fstatisticnilai - nilai.

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

Sumber: [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

Pilihan lain adalah menggunakan fungsi cor.test, alih-alih 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

Menggunakan:

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

di mana numadalah angka yang menunjukkan deretan matriks koefisien. Ini akan tergantung pada berapa banyak fitur yang Anda miliki dalam model Anda dan yang mana Anda ingin menarik nilai p. Misalnya, jika Anda hanya memiliki satu variabel, akan ada satu nilai p untuk intersep yang akan menjadi [1,4] dan yang berikutnya untuk variabel aktual Anda yang akan menjadi [2,4]. Jadi Anda numakan 2.