Quiero hacer una regresión lineal en R usando el lm() función. Mis datos es un evento anual de la serie de tiempo con un campo para el año (22 años) y otra para el estado (50 estados). Quiero ajuste de una regresión para cada estado, de modo que al final tengo un vector de respuestas de lm. Me puedo imaginar haciendo bucle for para cada estado, a continuación, haciendo la regresión dentro del bucle y la adición de los resultados de cada una regresión a un vector. Que no parece muy R-como, sin embargo. En el SAS me gustaría hacer un ‘por’ y de la instrucción en SQL me gustaría hacer un ‘grupo’. ¿Cuál es la R manera de hacerlo?

InformationsquelleAutor JD Long | 2009-07-23

10 Comentarios

  1. 42

    Aquí es una manera de usar el lme4 paquete.

     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
    • Es allí una manera a la lista de R2 para ambos de estos dos modelos? por ejemplo, para añadir un R2 de la columna tras año. Añadir también el valor de p para cada uno de los coef?
    • aquí usted puede encontrar una feaseble solución (mejor tarde que nunca): Su.df[,resumen(lm(Y~X))$r.cuadrado,por=Su.factor] donde: Y, X y Tu.factor son las variables. Por favor, tenga en cuenta que Su.el df debe ser una de datos.la tabla de la clase
  2. 57

    He aquí un método que utiliza la plyr paquete:

    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)
    • Decir que se agregó un adicional variable independiente que no estaba disponible en todos los estados (es decir millas.de.océano.de la costa), que fue representada por NA en sus datos. No la película error de llamada? ¿Cómo puede tratarse?
    • Dentro de la función que había necesidad de prueba para el caso y el uso de una fórmula diferente
    • Es posible añadir el nombre del subgrupo para cada llamada en el resumen (último paso)?
    • si ejecuta layout(matrix(c(1,2,3,4),2,2)) # optional 4 graphs/page y, a continuación, l_ply(models, plot) usted obtener cada uno de los grácos de residuos también. Es posible etiquetar cada una de las parcelas con el grupo (por ejemplo, «estado» en este caso)?
  3. 39

    Desde 2009, dplyr ha sido liberada, que en realidad proporciona una forma muy agradable de hacer este tipo de agrupación, muy parecidas a lo que SAS no.

    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  

    Recuperar los coeficientes y Rsquared/p.valor, se puede utilizar el broom paquete. Este paquete proporciona:

    tres S3 genéricos: ordenado, en el que se resume de un modelo
    resultados estadísticos tales como los coeficientes de una regresión;
    aumentar, lo que añade columnas para los datos originales, tales como
    las predicciones, los residuos y el clúster de las asignaturas; y la mirada, que
    proporciona una fila de resumen del modelo a nivel de estadísticas.

    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)
    • Tuve que hacer rowwise(fitted_models) %>% tidy(model) para obtener la escoba paquete para el trabajo, pero por lo demás, gran respuesta.
    • Funciona muy bien… puede hacer todo esto sin salir de la tubería: d %>% group_by(state) %>% do(model = lm(response ~ year, data = .)) %>% rowwise() %>% tidy(model)
  4. 23

    En mi opinión es un modelo lineal mixto un mejor enfoque para este tipo de datos. El código a continuación en el efecto fijo de la tendencia general. El de efectos aleatorios indican cómo la tendencia para cada uno de los estados difieren de la tendencia global. La correlación de la estructura de toma de la autocorrelación temporal en cuenta. Eche un vistazo a Pinheiro & Bates (Modelos de Efectos Mixtos en S y S-Plus).

    library(nlme)
    lme(response ~ year, random = ~year|state, correlation = corAR1(~year))
    • Esto es realmente un buen general de estadísticas de la teoría de la respuesta que me hace pensar en algunas cosas que no había considerado. La aplicación que me hizo la pregunta no sería aplicable para esta solución, pero me alegro de que te crió. Gracias.
    • No es una buena idea comenzar con un modelo mixto – ¿cómo sabe que alguno de los supuestos que están garantizados?
    • Debemos comprobar la hipótesis mediante la validación del modelo (y el conocimiento de los datos). Por CIERTO que no garantiza la asunción en el individuo, la película tampoco. Usted tendría que validar todos los modelos por separado.
  5. 12

    Una buena solución utilizando data.table fue publicado aquí en CrossValidated por @Zach.
    Solo quisiera agregar que es posible obtener de forma iterativa también el coeficiente de regresión 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

    así como todos los otros de salida 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
  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
  7. 6

    Yo ahora mi respuesta llega un poco tarde, pero yo estaba buscando una funcionalidad similar. Parecería que la función incorporada » por » en R también se puede hacer de la agrupación fácilmente:

    ?por contiene el siguiente ejemplo, que se adapta a cada grupo y extractos de los coeficientes con 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)
  8. 6

    Creo que merece la pena agregar el purrr::map enfoque a este problema.

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

    Ver @Pablo Hiemstra la respuesta para más ideas sobre el uso de la broom paquete con estos resultados.

    • Un poco de extensión en caso de que quieras una columna de valores ajustados o residuales: envoltura de la lm() la llamada en un resid() la llamada y, a continuación, tubo todo en la última línea en un unnest() la llamada. Por supuesto, usted desea cambiar el nombre de la variable de «modelo» a algo más relevante.
  9. 3

    La lm() función anterior es un ejemplo simple. Por cierto, me imagino que la base de datos de las columnas, como en la forma siguiente:

    año estado var1 var2 y…

    En mi punto de vista, puede utilizar el código siguiente:

    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)
  10. 0

    La pregunta parece ser acerca de cómo llamar a funciones de regresión con las fórmulas que se modifican dentro de un bucle.

    Aquí es cómo usted puede hacerlo en (el uso de los diamantes conjunto de datos):

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

Dejar respuesta

Please enter your comment!
Please enter your name here