Regresi Linier dan kelompok dalam R

97

Saya ingin melakukan regresi linier di R menggunakan lm()fungsi tersebut. Data saya adalah deret waktu tahunan dengan satu bidang untuk tahun (22 tahun) dan bidang lainnya untuk negara bagian (50 negara bagian). Saya ingin menyesuaikan regresi untuk setiap negara bagian sehingga pada akhirnya saya memiliki vektor tanggapan lm. Saya dapat membayangkan melakukan for loop untuk setiap negara bagian kemudian melakukan regresi di dalam loop dan menambahkan hasil setiap regresi ke vektor. Namun, itu tidak tampak seperti R. Dalam SAS saya akan melakukan pernyataan 'by' dan dalam SQL saya akan melakukan 'group by'. Apa cara R untuk melakukan ini?

r regression linear-regression lm JD Long
sumber
1
Hanya ingin memberi tahu orang-orang bahwa meskipun ada banyak fungsi kelompok-menurut di R, tidak semuanya tepat untuk regresi kelompok-demi. Misalnya, aggregatebukan orang yang benar ; tidak jugatapply .
李哲源

Jawaban:

51

Berikut salah satu cara menggunakan lme4paket tersebut.

 library(lme4)
 d <- data.frame(state=rep(c('NY', 'CA'), c(10, 10)),
                 year=rep(1:10, 2),
                 response=c(rnorm(10), rnorm(10)))

 xyplot(response ~ year, groups=state, data=d, type='l')

 fits <- lmList(response ~ year | state, data=d)
 fits
#------------
Call: lmList(formula = response ~ year | state, data = d)
Coefficients:
   (Intercept)        year
CA -1.34420990  0.17139963
NY  0.00196176 -0.01852429

Degrees of freedom: 20 total; 16 residual
Residual standard error: 0.8201316
ars
sumber
2
Apakah ada cara untuk membuat daftar R2 untuk kedua model ini? misalnya tambahkan kolom R2 setelah tahun. Tambahkan juga nilai p untuk masing-masing kopi?
ToToRo
3
@ToToRo di sini Anda dapat menemukan solusi feaseble (lebih baik terlambat daripada tidak sama sekali): Your.df [, ringkasan (lm (Y ~ X)) $ r.squared, by = Your.factor] di mana: Y, X dan Your.factor adalah variabel Anda. Harap diingat bahwa Your.df harus berupa kelas
data.table
60

Berikut pendekatan menggunakan paket plyr :

d <- data.frame(
  state = rep(c('NY', 'CA'), 10),
  year = rep(1:10, 2),
  response= rnorm(20)
)

library(plyr)
# Break up d by state, then fit the specified model to each piece and
# return a list
models <- dlply(d, "state", function(df) 
  lm(response ~ year, data = df))

# Apply coef to each model and return a data frame
ldply(models, coef)

# Print the summary of each model
l_ply(models, summary, .print = TRUE)
hadley
sumber
Misalnya Anda menambahkan variabel independen tambahan yang tidak tersedia di semua negara bagian (yaitu miles.of.ocean.shoreline) yang diwakili oleh NA dalam data Anda. Bukankah panggilan lm gagal? Bagaimana cara mengatasinya?
MikeTP
Di dalam fungsi, Anda perlu menguji kasus itu dan menggunakan rumus yang berbeda
hadley
Apakah mungkin untuk menambahkan nama subkelompok ke setiap panggilan dalam ringkasan (langkah terakhir)?
erc
jika Anda menjalankan layout(matrix(c(1,2,3,4),2,2)) # optional 4 graphs/page dan kemudian l_ply(models, plot)Anda mendapatkan setiap plot sisa juga. Apakah mungkin untuk memberi label pada setiap plot dengan grup (misalnya, "negara bagian" dalam kasus ini)?
Brian D
51

Sejak 2009, dplyrtelah dirilis yang sebenarnya memberikan cara yang sangat bagus untuk melakukan pengelompokan semacam ini, sangat mirip dengan apa yang dilakukan SAS.

library(dplyr)

d <- data.frame(state=rep(c('NY', 'CA'), c(10, 10)),
                year=rep(1:10, 2),
                response=c(rnorm(10), rnorm(10)))
fitted_models = d %>% group_by(state) %>% do(model = lm(response ~ year, data = .))
# Source: local data frame [2 x 2]
# Groups: <by row>
#
#    state   model
#   (fctr)   (chr)
# 1     CA <S3:lm>
# 2     NY <S3:lm>
fitted_models$model
# [[1]]
# 
# Call:
# lm(formula = response ~ year, data = .)
# 
# Coefficients:
# (Intercept)         year  
#    -0.06354      0.02677  
#
#
# [[2]]
# 
# Call:
# lm(formula = response ~ year, data = .)
# 
# Coefficients:
# (Intercept)         year  
#    -0.35136      0.09385

Untuk mengambil koefisien dan nilai Rsquared / p., seseorang dapat menggunakan broompaket tersebut. Paket ini menyediakan:

tiga obat generik S3: rapi, yang merangkum temuan statistik model seperti koefisien regresi; augment, yang menambahkan kolom ke data asli seperti prediksi, residual, dan tugas cluster; dan glance, yang memberikan ringkasan satu baris statistik tingkat model.

library(broom)
fitted_models %>% tidy(model)
# Source: local data frame [4 x 6]
# Groups: state [2]
# 
#    state        term    estimate  std.error  statistic   p.value
#   (fctr)       (chr)       (dbl)      (dbl)      (dbl)     (dbl)
# 1     CA (Intercept) -0.06354035 0.83863054 -0.0757668 0.9414651
# 2     CA        year  0.02677048 0.13515755  0.1980687 0.8479318
# 3     NY (Intercept) -0.35135766 0.60100314 -0.5846187 0.5749166
# 4     NY        year  0.09385309 0.09686043  0.9689519 0.3609470
fitted_models %>% glance(model)
# Source: local data frame [2 x 12]
# Groups: state [2]
# 
#    state   r.squared adj.r.squared     sigma statistic   p.value    df
#   (fctr)       (dbl)         (dbl)     (dbl)     (dbl)     (dbl) (int)
# 1     CA 0.004879969  -0.119510035 1.2276294 0.0392312 0.8479318     2
# 2     NY 0.105032068  -0.006838924 0.8797785 0.9388678 0.3609470     2
# Variables not shown: logLik (dbl), AIC (dbl), BIC (dbl), deviance (dbl),
#   df.residual (int)
fitted_models %>% augment(model)
# Source: local data frame [20 x 10]
# Groups: state [2]
# 
#     state   response  year      .fitted   .se.fit     .resid      .hat
#    (fctr)      (dbl) (int)        (dbl)     (dbl)      (dbl)     (dbl)
# 1      CA  0.4547765     1 -0.036769875 0.7215439  0.4915464 0.3454545
# 2      CA  0.1217003     2 -0.009999399 0.6119518  0.1316997 0.2484848
# 3      CA -0.6153836     3  0.016771076 0.5146646 -0.6321546 0.1757576
# 4      CA -0.9978060     4  0.043541551 0.4379605 -1.0413476 0.1272727
# 5      CA  2.1385614     5  0.070312027 0.3940486  2.0682494 0.1030303
# 6      CA -0.3924598     6  0.097082502 0.3940486 -0.4895423 0.1030303
# 7      CA -0.5918738     7  0.123852977 0.4379605 -0.7157268 0.1272727
# 8      CA  0.4671346     8  0.150623453 0.5146646  0.3165112 0.1757576
# 9      CA -1.4958726     9  0.177393928 0.6119518 -1.6732666 0.2484848
# 10     CA  1.7481956    10  0.204164404 0.7215439  1.5440312 0.3454545
# 11     NY -0.6285230     1 -0.257504572 0.5170932 -0.3710185 0.3454545
# 12     NY  1.0566099     2 -0.163651479 0.4385542  1.2202614 0.2484848
# 13     NY -0.5274693     3 -0.069798386 0.3688335 -0.4576709 0.1757576
# 14     NY  0.6097983     4  0.024054706 0.3138637  0.5857436 0.1272727
# 15     NY -1.5511940     5  0.117907799 0.2823942 -1.6691018 0.1030303
# 16     NY  0.7440243     6  0.211760892 0.2823942  0.5322634 0.1030303
# 17     NY  0.1054719     7  0.305613984 0.3138637 -0.2001421 0.1272727
# 18     NY  0.7513057     8  0.399467077 0.3688335  0.3518387 0.1757576
# 19     NY -0.1271655     9  0.493320170 0.4385542 -0.6204857 0.2484848
# 20     NY  1.2154852    10  0.587173262 0.5170932  0.6283119 0.3454545
# Variables not shown: .sigma (dbl), .cooksd (dbl), .std.resid (dbl)
Paul Hiemstra
sumber
2
Saya harus melakukannya rowwise(fitted_models) %>% tidy(model)agar paket sapu berfungsi, tetapi sebaliknya, jawaban yang bagus.
pedram
3
Bekerja dengan baik ... dapat melakukan ini semua tanpa meninggalkan pipa:d %>% group_by(state) %>% do(model = lm(response ~ year, data = .)) %>% rowwise() %>% tidy(model)
holastello
1
@pedram dan @holastello, ini tidak berfungsi lagi, setidaknya dengan R 3.6.1, broom_0.7.0, dplyr_0.8.3. d %>% group_by(state) %>% do(model=lm(response ~year, data = .)) %>% rowwise() %>% tidy(model) Error in var(if (is.vector(x) || is.factor(x)) x else as.double(x), na.rm = na.rm) : Calling var(x) on a factor x is defunct. Use something like 'all(duplicated(x)[-1L])' to test for a constant vector. In addition: Warning messages: 1: Data frame tidiers are deprecated and will be removed in an upcoming release of broom. ...
Chris Nolte
24

Menurut pendapat saya adalah model linier campuran pendekatan yang lebih baik untuk jenis data ini. Kode di bawah ini diberikan dalam efek tetap tren keseluruhan. Efek acak menunjukkan bagaimana tren untuk setiap negara bagian berbeda dari tren global. Struktur korelasi memperhitungkan autokorelasi temporal. Lihat Pinheiro & Bates (Model Efek Campuran di S dan S-Plus).

library(nlme)
lme(response ~ year, random = ~year|state, correlation = corAR1(~year))
Thierry
sumber
3
Ini adalah jawaban teori statistik umum yang sangat bagus yang membuat saya berpikir tentang beberapa hal yang belum saya pertimbangkan. Aplikasi yang membuat saya mengajukan pertanyaan tidak akan berlaku untuk solusi ini, tetapi saya senang Anda mengemukakannya. Terima kasih.
JD Long
1
Bukan ide yang baik untuk memulai dengan model campuran - bagaimana Anda tahu bahwa salah satu asumsi dijamin?
hadley
8
Seseorang harus memeriksa asumsi dengan validasi model (dan pengetahuan tentang data). BTW Anda juga tidak dapat menjamin asumsi pada individu lm. Anda harus memvalidasi semua model secara terpisah.
Thierry
17

Solusi yang bagus menggunakan data.tablediposting di sini di CrossValidated oleh @Zach. Saya hanya menambahkan bahwa dimungkinkan untuk mendapatkan secara iteratif juga koefisien regresi r ^ 2:

## make fake data
    library(data.table)
    set.seed(1)
    dat <- data.table(x=runif(100), y=runif(100), grp=rep(1:2,50))

##calculate the regression coefficient r^2
    dat[,summary(lm(y~x))$r.squared,by=grp]
       grp         V1
    1:   1 0.01465726
    2:   2 0.02256595

serta semua keluaran lainnya dari summary(lm):

dat[,list(r2=summary(lm(y~x))$r.squared , f=summary(lm(y~x))$fstatistic[1] ),by=grp]
   grp         r2        f
1:   1 0.01465726 0.714014
2:   2 0.02256595 1.108173
FraNut
sumber
8

Saya pikir ada baiknya menambahkan purrr::mappendekatan untuk masalah ini.

library(tidyverse)

d <- data.frame(state=rep(c('NY', 'CA'), c(10, 10)),
                                 year=rep(1:10, 2),
                                 response=c(rnorm(10), rnorm(10)))

d %>% 
  group_by(state) %>% 
  nest() %>% 
  mutate(model = map(data, ~lm(response ~ year, data = .)))

Lihat jawaban @Paul Hiemstra untuk ide lebih lanjut tentang menggunakan broompaket dengan hasil ini.

ngm
sumber
Sedikit ekstensi jika Anda menginginkan kolom berisi nilai atau residu: bungkus panggilan lm () dalam panggilan resid () dan kemudian pipa semua yang ada di baris terakhir menjadi panggilan tidak terestimasi (). Tentu saja, Anda ingin mengubah nama variabel dari "model" menjadi sesuatu yang lebih relevan.
randy
8 
## make fake data
 ngroups <- 2
 group <- 1:ngroups
 nobs <- 100
 dta <- data.frame(group=rep(group,each=nobs),y=rnorm(nobs*ngroups),x=runif(nobs*ngroups))
 head(dta)
#--------------------
  group          y         x
1     1  0.6482007 0.5429575
2     1 -0.4637118 0.7052843
3     1 -0.5129840 0.7312955
4     1 -0.6612649 0.9028034
5     1 -0.5197448 0.1661308
6     1  0.4240346 0.8944253
#------------ 
## function to extract the results of one model
 foo <- function(z) {
   ## coef and se in a data frame
   mr <- data.frame(coef(summary(lm(y~x,data=z))))
   ## put row names (predictors/indep variables)
   mr$predictor <- rownames(mr)
   mr
 }
 ## see that it works
 foo(subset(dta,group==1))
#=========
              Estimate Std..Error   t.value  Pr...t..   predictor
(Intercept)  0.2176477  0.1919140  1.134090 0.2595235 (Intercept)
x           -0.3669890  0.3321875 -1.104765 0.2719666           x
#----------
## one option: use command by
 res <- by(dta,dta$group,foo)
 res
#=========
dta$group: 1
              Estimate Std..Error   t.value  Pr...t..   predictor
(Intercept)  0.2176477  0.1919140  1.134090 0.2595235 (Intercept)
x           -0.3669890  0.3321875 -1.104765 0.2719666           x
------------------------------------------------------------ 
dta$group: 2
               Estimate Std..Error    t.value  Pr...t..   predictor
(Intercept) -0.04039422  0.1682335 -0.2401081 0.8107480 (Intercept)
x            0.06286456  0.3020321  0.2081387 0.8355526           x

## using package plyr is better
 library(plyr)
 res <- ddply(dta,"group",foo)
 res
#----------
  group    Estimate Std..Error    t.value  Pr...t..   predictor
1     1  0.21764767  0.1919140  1.1340897 0.2595235 (Intercept)
2     1 -0.36698898  0.3321875 -1.1047647 0.2719666           x
3     2 -0.04039422  0.1682335 -0.2401081 0.8107480 (Intercept)
4     2  0.06286456  0.3020321  0.2081387 0.8355526           x
Eduardo Leoni
sumber
7

Saya sekarang jawaban saya datang agak terlambat, tetapi saya sedang mencari fungsi serupa. Tampaknya fungsi bawaan 'oleh' di R juga dapat melakukan pengelompokan dengan mudah:

? by berisi contoh berikut, yang sesuai untuk setiap grup dan mengekstrak koefisien dengan sapply:

require(stats)
## now suppose we want to extract the coefficients by group 
tmp <- with(warpbreaks,
            by(warpbreaks, tension,
               function(x) lm(breaks ~ wool, data = x)))
sapply(tmp, coef)
Matthijs Cox
sumber
4

The lm()fungsi di atas adalah contoh sederhana. By the way, saya membayangkan bahwa database Anda memiliki kolom seperti pada formulir berikut:

tahun negara var1 var2 y ...

Dalam pandangan saya, Anda dapat menggunakan kode berikut:

require(base) 
library(base) 
attach(data) # data = your data base
             #state is your label for the states column
modell<-by(data, data$state, function(data) lm(y~I(1/var1)+I(1/var2)))
summary(modell)
Zack Mendes
sumber
0

Pertanyaannya sepertinya tentang bagaimana memanggil fungsi regresi dengan rumus yang dimodifikasi di dalam loop.

Berikut adalah bagaimana Anda dapat melakukannya di (menggunakan dataset berlian):

attach(ggplot2::diamonds)
strCols = names(ggplot2::diamonds)

formula <- list(); model <- list()
for (i in 1:1) {
  formula[[i]] = paste0(strCols[7], " ~ ", strCols[7+i])
  model[[i]] = glm(formula[[i]]) 

  #then you can plot the results or anything else ...
  png(filename = sprintf("diamonds_price=glm(%s).png", strCols[7+i]))
  par(mfrow = c(2, 2))      
  plot(model[[i]])
  dev.off()
  }
IVIM
sumber