Régression linéaire et grouper par dans R

Question 1

Je veux faire une régression linéaire dans R en utilisant la lm()fonction. Mes données sont une série chronologique annuelle avec un champ pour l'année (22 ans) et un autre pour l'état (50 états). Je veux ajuster une régression pour chaque état afin qu'à la fin j'aie un vecteur de réponses lm. Je peux imaginer faire une boucle for pour chaque état puis faire la régression à l'intérieur de la boucle et ajouter les résultats de chaque régression à un vecteur. Cela ne semble pas très semblable à R, cependant. En SAS, je ferais une instruction «par» et en SQL, je ferais un «groupe par». Quelle est la manière R de faire cela?

Question 2

Voici une façon d'utiliser le lme4package.

 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

Question 3

Voici une approche utilisant le package 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)

Question 4

Depuis 2009, dplyra été publié, ce qui offre en fait une très belle façon de faire ce type de regroupement, ressemblant étroitement à ce que fait 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

Pour récupérer les coefficients et Rsquared / p.value, on peut utiliser le broompackage. Ce package fournit:

trois génériques S3: tidy, qui résume les résultats statistiques d'un modèle tels que les coefficients d'une régression; augment, qui ajoute des colonnes aux données d'origine telles que les prédictions, les résidus et les affectations de cluster; et glance, qui fournit un résumé sur une ligne des statistiques au niveau du modèle.

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)

Question 5

À mon avis, un modèle linéaire mixte est une meilleure approche pour ce type de données. Le code ci-dessous donne à effet fixe la tendance générale. Les effets aléatoires indiquent comment la tendance pour chaque état individuel diffère de la tendance mondiale. La structure de corrélation prend en compte l'autocorrélation temporelle. Jetez un œil à Pinheiro & Bates (modèles d'effets mixtes en S et S-Plus).

library(nlme)
lme(response ~ year, random = ~year|state, correlation = corAR1(~year))

Question 6

Une belle solution utilisant a data.tableété publiée ici dans CrossValidated by @Zach. J'ajouterais simplement qu'il est possible d'obtenir itérativement aussi le coefficient de régression 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

ainsi que toutes les autres sorties de 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

Question 7

Je pense qu'il vaut la peine d'ajouter l' purrr::mapapproche à ce problème.

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

Voir la réponse de @Paul Hiemstra pour plus d'idées sur l'utilisation du broompackage avec ces résultats.

Question 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

Question 9

Maintenant, ma réponse arrive un peu tard, mais je cherchais une fonctionnalité similaire. Il semblerait que la fonction intégrée `` by '' dans R puisse également effectuer le regroupement facilement:

? by contient l'exemple suivant, qui tient par groupe et extrait les coefficients avec 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)

Question 10

La lm()fonction ci-dessus est un exemple simple. Au fait, j'imagine que votre base de données a les colonnes sous la forme suivante:

année état var1 var2 y ...

De mon point de vue, vous pouvez utiliser le code suivant:

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)

Question 11

La question semble être de savoir comment appeler des fonctions de régression avec des formules modifiées à l'intérieur d'une boucle.

Voici comment vous pouvez le faire (en utilisant un jeu de données de diamants):

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()
  }

Answer 1

Je veux faire une régression linéaire dans R en utilisant la lm()fonction. Mes données sont une série chronologique annuelle avec un champ pour l'année (22 ans) et un autre pour l'état (50 états). Je veux ajuster une régression pour chaque état afin qu'à la fin j'aie un vecteur de réponses lm. Je peux imaginer faire une boucle for pour chaque état puis faire la régression à l'intérieur de la boucle et ajouter les résultats de chaque régression à un vecteur. Cela ne semble pas très semblable à R, cependant. En SAS, je ferais une instruction «par» et en SQL, je ferais un «groupe par». Quelle est la manière R de faire cela?

Answer 2

1

Je veux juste dire aux gens que bien qu'il y ait beaucoup de fonctions group-by dans R, toutes ne sont pas les bonnes pour la régression group-by. Par exemple, aggregaten'est pas une bonne solution ; ni l'un ni l'autretapply .

李哲源

Answer 3

48

Voici une façon d'utiliser le lme4package.

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

2

Existe-t-il un moyen de lister R2 pour ces deux modèles? par exemple, ajoutez une colonne R2 après année. Ajoutez également la valeur p pour chacun des coeff?

ToToRo

3

@ToToRo ici vous pouvez trouver une solution flexible (mieux vaut tard que jamais): Your.df [, summary (lm (Y ~ X)) $ r.squared, by = Your.factor] où: Y, X et Your.factor sont vos variables. Veuillez garder à l'esprit que Your.df doit être une classe

data.table

Answer 4

2

Existe-t-il un moyen de lister R2 pour ces deux modèles? par exemple, ajoutez une colonne R2 après année. Ajoutez également la valeur p pour chacun des coeff?

ToToRo

Answer 5

3

@ToToRo ici vous pouvez trouver une solution flexible (mieux vaut tard que jamais): Your.df [, summary (lm (Y ~ X)) $ r.squared, by = Your.factor] où: Y, X et Your.factor sont vos variables. Veuillez garder à l'esprit que Your.df doit être une classe

data.table

Answer 6

60

Voici une approche utilisant le package 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
la source

Supposons que vous ayez ajouté une variable indépendante supplémentaire qui n'était pas disponible dans tous les états (c'est-à-dire miles.of.ocean.shoreline) qui était représentée par NA dans vos données. L'appel lm n'échouerait-il pas? Comment pourrait-il être traité?

MikeTP

À l'intérieur de la fonction, vous devez tester ce cas et utiliser une formule différente

hadley

Est-il possible d'ajouter le nom du sous-groupe à chaque appel dans le résumé (dernière étape)?

betterave

si vous exécutez layout(matrix(c(1,2,3,4),2,2)) # optional 4 graphs/page et que l_ply(models, plot)vous obtenez également chacun des tracés des résidus. Est-il possible d'étiqueter chacune des parcelles avec le groupe (par exemple, «état» dans ce cas)?

Brian D

Answer 7

Supposons que vous ayez ajouté une variable indépendante supplémentaire qui n'était pas disponible dans tous les états (c'est-à-dire miles.of.ocean.shoreline) qui était représentée par NA dans vos données. L'appel lm n'échouerait-il pas? Comment pourrait-il être traité?

MikeTP

Answer 8

À l'intérieur de la fonction, vous devez tester ce cas et utiliser une formule différente

hadley

Answer 9

Est-il possible d'ajouter le nom du sous-groupe à chaque appel dans le résumé (dernière étape)?

betterave

Answer 10

si vous exécutez layout(matrix(c(1,2,3,4),2,2)) # optional 4 graphs/page et que l_ply(models, plot)vous obtenez également chacun des tracés des résidus. Est-il possible d'étiqueter chacune des parcelles avec le groupe (par exemple, «état» dans ce cas)?

Brian D

Answer 11

Depuis 2009, dplyra été publié, ce qui offre en fait une très belle façon de faire ce type de regroupement, ressemblant étroitement à ce que fait 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

Pour récupérer les coefficients et Rsquared / p.value, on peut utiliser le broompackage. Ce package fournit:

trois génériques S3: tidy, qui résume les résultats statistiques d'un modèle tels que les coefficients d'une régression; augment, qui ajoute des colonnes aux données d'origine telles que les prédictions, les résidus et les affectations de cluster; et glance, qui fournit un résumé sur une ligne des statistiques au niveau du modèle.

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)

Answer 12

2

J'ai dû faire rowwise(fitted_models) %>% tidy(model)pour que le paquet de balai fonctionne, mais sinon, excellente réponse.

pedram

Answer 13

3

Fonctionne très bien ... peut faire tout cela sans quitter le tuyau:d %>% group_by(state) %>% do(model = lm(response ~ year, data = .)) %>% rowwise() %>% tidy(model)

holastello

Answer 14

@pedram et @holastello, cela ne fonctionne plus, du moins avec 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

Answer 15

23

À mon avis, un modèle linéaire mixte est une meilleure approche pour ce type de données. Le code ci-dessous donne à effet fixe la tendance générale. Les effets aléatoires indiquent comment la tendance pour chaque état individuel diffère de la tendance mondiale. La structure de corrélation prend en compte l'autocorrélation temporelle. Jetez un œil à Pinheiro & Bates (modèles d'effets mixtes en S et S-Plus).

library(nlme)
lme(response ~ year, random = ~year|state, correlation = corAR1(~year))

Thierry
la source

3

C'est une très bonne réponse à la théorie des statistiques générales qui me fait réfléchir à certaines choses que je n'avais pas envisagées. L'application qui m'a amené à poser la question ne serait pas applicable à cette solution, mais je suis heureux que vous en ayez parlé. Je vous remercie.

JD Long

1

Ce n'est pas une bonne idée de commencer avec un modèle mixte - comment savez-vous que l'une des hypothèses est justifiée?

hadley

7

Il convient de vérifier l'hypothèse par la validation du modèle (et la connaissance des données). BTW, vous ne pouvez pas non plus garantir l'hypothèse sur les lm individuels. Vous devrez valider tous les modèles séparément.

Thierry

Answer 16

3

C'est une très bonne réponse à la théorie des statistiques générales qui me fait réfléchir à certaines choses que je n'avais pas envisagées. L'application qui m'a amené à poser la question ne serait pas applicable à cette solution, mais je suis heureux que vous en ayez parlé. Je vous remercie.

JD Long

Answer 17

1

Ce n'est pas une bonne idée de commencer avec un modèle mixte - comment savez-vous que l'une des hypothèses est justifiée?

hadley

Answer 18

7

Il convient de vérifier l'hypothèse par la validation du modèle (et la connaissance des données). BTW, vous ne pouvez pas non plus garantir l'hypothèse sur les lm individuels. Vous devrez valider tous les modèles séparément.

Thierry

Answer 19

Une belle solution utilisant a data.tableété publiée ici dans CrossValidated by @Zach. J'ajouterais simplement qu'il est possible d'obtenir itérativement aussi le coefficient de régression 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

ainsi que toutes les autres sorties de 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

Answer 20

8

Je pense qu'il vaut la peine d'ajouter l' purrr::mapapproche à ce problème.

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

Voir la réponse de @Paul Hiemstra pour plus d'idées sur l'utilisation du broompackage avec ces résultats.

ngm
la source

Une petite extension au cas où vous voudriez une colonne de valeurs ajustées ou de résidus: encapsulez l'appel lm () dans un appel resid () puis redirigez tout ce qui se trouve dans la dernière ligne dans un appel unnest (). Bien sûr, vous voudrez changer le nom de la variable de "modèle" à quelque chose de plus pertinent.

randy

Answer 21

Une petite extension au cas où vous voudriez une colonne de valeurs ajustées ou de résidus: encapsulez l'appel lm () dans un appel resid () puis redirigez tout ce qui se trouve dans la dernière ligne dans un appel unnest (). Bien sûr, vous voudrez changer le nom de la variable de "modèle" à quelque chose de plus pertinent.

randy

Answer 22

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

Answer 23

Maintenant, ma réponse arrive un peu tard, mais je cherchais une fonctionnalité similaire. Il semblerait que la fonction intégrée `` by '' dans R puisse également effectuer le regroupement facilement:

? by contient l'exemple suivant, qui tient par groupe et extrait les coefficients avec 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)

Answer 24

La lm()fonction ci-dessus est un exemple simple. Au fait, j'imagine que votre base de données a les colonnes sous la forme suivante:

année état var1 var2 y ...

De mon point de vue, vous pouvez utiliser le code suivant:

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)

Answer 25

La question semble être de savoir comment appeler des fonctions de régression avec des formules modifiées à l'intérieur d'une boucle.

Voici comment vous pouvez le faire (en utilisant un jeu de données de diamants):

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()
  }

Régression linéaire et grouper par dans R

Réponses: