La Régression linéaire et de groupe par dans la R

Je veux faire une régression linéaire dans R à l'aide de la lm() fonction. Mes données est un temps annuel de la série avec un champ pour l'année (22 ans) et un autre pour l'état (50 membres). Je veux ajustement d'une droite de régression pour chaque état, de sorte qu'à la fin j'ai un vecteur de réponses lm. Je peux imaginer faire de boucle pour chaque état ensuite de faire de la régression à l'intérieur de la boucle, et en additionnant les résultats de chaque régression à un vecteur. Cela ne semble pas très R-comme, cependant. Dans le SAS, je ferais une 'par' énoncé, et en SQL, je voudrais faire un 'group by'. Quel est le R moyen de faire cela?

InformationsquelleAutor JD Long | 2009-07-23
  1. 42

    Voici un moyen à l'aide de la lme4 paquet.

     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
    • Est-il un moyen de la liste de R2 pour ces deux modèles? par exemple, ajouter un R2 de la colonne après année. Également ajouter de la p-valeur pour chacun des coeff?
    • ici vous pouvez trouver une feaseble solution de (mieux vaut tard que jamais): Votre.df[,summary(lm(Y~X))$r.au carré,par=Votre.facteur] où: Y, X et Votre.facteur sont variables. Veuillez garder à l'esprit que Votre.df doit être un ensemble de données.le tableau de la classe
    InformationsquelleAutor ars
  2. 57

    Ici est une approche à l'aide de la plyr package:

    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)
    • Dire que vous avez ajouté un supplémentaire variable indépendante, qui n'est pas disponible dans tous les états (c'est à dire des miles.de.de l'océan.littoral), qui était représentée par NA dans vos données. Ne serait pas le lm appel échoue? Comment pourrait-elle être traitée?
    • L'intérieur de la fonction, vous aurez besoin de test pour ce cas, et d'utiliser une formule différente
    • Est-il possible d'ajouter le nom du sous-groupe à chaque appel dans le résumé (dernière étape)?
    • si vous exécutez layout(matrix(c(1,2,3,4),2,2)) # optional 4 graphs/page et puis l_ply(models, plot) vous obtenez chacun des résidus de parcelles trop. Est-il possible de l'étiquette de chacune des parcelles avec le groupe (par exemple, un "état" dans ce cas)?
    InformationsquelleAutor hadley
  3. 38

    Depuis 2009, dplyr a été publié, qui fournit en fait une façon très agréable de faire ce genre de regroupement, qui ressemblent beaucoup à ce que SAS n'.

    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

    Extraire les coefficients et Rsquared/p.valeur, on peut utiliser la broom paquet. Ce paquet fournit:

    trois S3 génériques: bien rangé, qui résume un modèle de
    résultats statistiques tels que les coefficients de régression;
    augmenter, ce qui ajoute des colonnes à l'origine des données telles que
    les prédictions, les résidus et le cluster des travaux; et le regard, qui
    fournit une ligne de résumé de modèle au niveau des statistiques.

    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)
    • J'ai eu à faire rowwise(fitted_models) %>% tidy(model) pour obtenir le balai paquet de travailler, mais sinon, super réponse.
    • 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)
    InformationsquelleAutor Paul Hiemstra
  4. 23

    À mon avis, est un modèle linéaire mixte une meilleure approche pour ce type de données. Le code ci-dessous dans le fixe, l'effet de la tendance générale. Les effets aléatoires indiquer comment la tendance pour chaque individu de l'état diffèrent de la tendance mondiale. La structure de corrélation prend l'auto-corrélation temporelle en compte. Jetez un oeil à Pinheiro & Bates (à Effets Mixtes Modèles S et S-Plus).

    library(nlme)
lme(response ~ year, random = ~year|state, correlation = corAR1(~year))
    • C'est une très bonne stats théorie de la réponse qui m'a fait penser à certaines choses que je n'avais pas envisagé. La demande qui m'a amené à poser la question ne serait pas applicable à cette solution, mais je suis heureux de l'avoir fait. Je vous remercie.
    • Ce n'est pas une bonne idée de commencer avec un modèle mixte - comment savez-vous que les hypothèses sont garantis?
    • Il faut vérifier l'hypothèse par la validation du modèle (et la connaissance des données). BTW, vous ne pouvez pas justifier l'hypothèse sur la personne de lm soit. Vous devez valider tous les modèles séparément.
    InformationsquelleAutor Thierry
  5. 12

    Une belle solution à l'aide de data.table a été posté ici dans CrossValidated par @Zach.
    Je voudrais juste ajouter qu'il est possible d'obtenir de manière itérative également 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
    InformationsquelleAutor FraNut
  6. 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
    InformationsquelleAutor Eduardo Leoni
  7. 5

    J'ai maintenant ma réponse vient un peu tard, mais je cherchais une fonctionnalité similaire. Il semble que la fonction intégrée " par " dans R peut aussi faire le regroupement facilement:

    ?par contient l'exemple suivant, qui correspond par groupe et extraits de 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)
    InformationsquelleAutor Matthijs Cox
  8. 5

    Je pense qu'il est intéressant d'ajouter la purrr::map approche à 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 @Paul Hiemstra de réponse pour d'autres idées sur l'utilisation de la broom paquet avec ces résultats.

    InformationsquelleAutor ngm
  9. 3

    La lm() la fonction ci-dessus est un exemple simple. En passant, j'imagine que votre base de données comporte les colonnes comme dans la forme suivante:

    année de l'état var1 var2 y...

    À 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)
    InformationsquelleAutor Zack Mendes
  10. 0

    La question semble être sur la façon d'appeler les fonctions de régression avec des formules qui sont modifiés à l'intérieur d'une boucle.

    Voici comment vous pouvez le faire à l'aide de diamants dataset):

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