Regresja liniowa i grupowanie według w R.

97

Chcę wykonać regresję liniową w R za pomocą lm()funkcji. Moje dane to roczny szereg czasowy z jednym polem dla roku (22 lata) i drugim dla stanu (50 stanów). Chcę dopasować regresję dla każdego stanu, tak aby na końcu mieć wektor odpowiedzi lm. Mogę sobie wyobrazić wykonanie pętli for dla każdego stanu, a następnie wykonanie regresji wewnątrz pętli i dodanie wyników każdej regresji do wektora. Nie wydaje się to jednak zbyt podobne do R. W SAS zrobiłbym instrukcję „by”, aw SQL zrobiłbym „group by”. Jak to zrobić w R?

r  regression  linear-regression  lm 
— JD Long
źródło
1
Chcę tylko powiedzieć ludziom, że chociaż istnieje wiele funkcji grupowania w R, nie wszystkie z nich są odpowiednie do regresji grupowania. Na przykład aggregatenie jest właściwe ; nie jesttapply .
— 李哲源

Odpowiedzi:

50

Oto jeden sposób korzystania z lme4pakietu.

 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
źródło
2
Czy istnieje sposób na wyświetlenie R2 dla obu tych dwóch modeli? np. dodaj kolumnę R2 po roku. Dodaj także wartość p dla każdego ze współczynników?
— ToToRo
3
@ToToRo tutaj możesz znaleźć wykonalne rozwiązanie (lepiej późno niż wcale): Your.df [, summary (lm (Y ~ X)) $ r.squared, by = Your.factor] gdzie: Y, X i Your.factor są twoimi zmiennymi. Należy pamiętać, że Your.df musi być klasą
— data.table
60

Oto podejście wykorzystujące pakiet 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
źródło
Załóżmy, że dodałeś dodatkową niezależną zmienną, która nie była dostępna we wszystkich stanach (tj. Mile.of.ocean.shoreline), która była reprezentowana przez NA w Twoich danych. Czy połączenie lm nie zawiodło? Jak można sobie z tym poradzić?
— MikeTP
Wewnątrz funkcji musisz przetestować ten przypadek i użyć innej formuły
— hadley
Czy można dodać nazwę podgrupy do każdego zaproszenia w podsumowaniu (ostatni krok)?
— burak
jeśli uciekniesz, layout(matrix(c(1,2,3,4),2,2)) # optional 4 graphs/page a następnie l_ply(models, plot)otrzymasz również każdy z wykresów reszt. Czy można oznaczyć każdą z działek grupą (np. „Stan” w tym przypadku)?
— Brian D
51

Od 2009 roku dplyrzostał wydany, co w rzeczywistości zapewnia bardzo przyjemny sposób na tego rodzaju grupowanie, bardzo przypominający to, co robi 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

Aby pobrać współczynniki i Rsquared / p.value, można użyć broompakietu. Ten pakiet obejmuje:

trzy typy generyczne S3: schludny, który podsumowuje wyniki statystyczne modelu, takie jak współczynniki regresji; rozszerzenie, które dodaje kolumny do pierwotnych danych, takich jak prognozy, reszty i przypisania klastrów; i rzut oka, który zawiera jednowierszowe podsumowanie statystyk na poziomie modelu.

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
źródło
2
Musiałem zrobić, rowwise(fitted_models) %>% tidy(model)aby pakiet miotły działał, ale poza tym świetna odpowiedź.
— pedram
3
Działa świetnie ... można to wszystko zrobić bez wychodzenia z rury:d %>% group_by(state) %>% do(model = lm(response ~ year, data = .)) %>% rowwise() %>% tidy(model)
— holastello
1
@pedram i @holastello, to już nie działa, przynajmniej z 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

Moim zdaniem mieszany model liniowy jest lepszym podejściem do tego rodzaju danych. Poniższy kod podaje stały efekt ogólnej tendencji. Losowe efekty wskazują, jak trend dla każdego stanu różni się od trendu globalnego. Struktura korelacji uwzględnia czasową autokorelację. Spójrz na Pinheiro & Bates (modele z efektami mieszanymi w S i S-Plus).

library(nlme)
lme(response ~ year, random = ~year|state, correlation = corAR1(~year))
— Thierry
źródło
3
To jest naprawdę dobra ogólna teoria statystyk, która każe mi myśleć o niektórych rzeczach, których nie brałem pod uwagę. Aplikacja, która spowodowała, że ​​zadałem to pytanie nie miałaby zastosowania do tego rozwiązania, ale cieszę się, że ją poruszyłaś. Dziękuję Ci.
— JD Long
1
Nie jest dobrym pomysłem rozpoczynanie od modelu mieszanego - skąd wiesz, że którekolwiek z założeń jest uzasadnione?
— hadley
8
Założenie należy sprawdzić poprzez walidację modelu (i znajomość danych). A tak przy okazji, nie możesz również zagwarantować założenia na temat poszczególnych lm. Musisz osobno zweryfikować wszystkie modele.
— Thierry,
17

Przyjemne rozwiązanie wykorzystujące data.tablezostało opublikowane tutaj w CrossValidated przez @Zach. Dodam tylko, że można iteracyjnie otrzymać także współczynnik regresji 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

a także wszystkie inne dane wyjściowe z 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
źródło
8

Myślę, że warto dodać purrr::mappodejście do tego problemu.

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 = .)))

Zobacz odpowiedź @Paul Hiemstra, aby zapoznać się z dalszymi pomysłami na wykorzystanie broompakietu z tymi wynikami.

— ngm
źródło
Małe rozszerzenie na wypadek, gdybyś potrzebował kolumny dopasowanych wartości lub reszt: zawiń wywołanie lm () w wywołanie reszty (), a następnie potokuj wszystko w ostatniej linii do wywołania unnest (). Oczywiście chciałbyś zmienić nazwę zmiennej z „model” na bardziej odpowiednią.
— 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
źródło
7

Teraz moja odpowiedź przychodzi trochę za późno, ale szukałem podobnej funkcjonalności. Wydawałoby się, że wbudowana funkcja „by” w R może również łatwo wykonać grupowanie:

? by zawiera następujący przykład, który pasuje do grupy i wyodrębnia współczynniki za pomocą 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
źródło
4

Powyższa lm()funkcja jest prostym przykładem. Swoją drogą wyobrażam sobie, że Twoja baza danych ma kolumny jak w następującej formie:

rok stan var1 var2 y ...

Z mojego punktu widzenia możesz użyć następującego kodu:

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
źródło
0

Wydaje się, że pytanie dotyczy tego, jak wywołać funkcje regresji za pomocą formuł, które są modyfikowane w pętli.

Oto, jak możesz to zrobić w (używając zestawu danych diamentów):

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
źródło
Korzystając z naszej strony potwierdzasz, że przeczytałeś(-aś) i rozumiesz nasze zasady używania plików cookie i zasady ochrony prywatności.
Licensed under cc by-sa 3.0 with attribution required.