WG: grouped anova from lm as data frame, tidyr broom

Previous Topic Next Topic
 
classic Classic list List threaded Threaded
1 message Options
Reply | Threaded
Open this post in threaded view
|

WG: grouped anova from lm as data frame, tidyr broom

Kluth, Christian
Dear R-project team,
I got a message that the e-mail from this morning was filtered away completely. So hopefully this is working.
I am trying to solve a probably easy problem with R using the tidyr package. Probably I did not get the concept very clearly.
I would like to get a nice tidy data frame of ANOVAs of several identical models (lm or lme (anovas then generated with the car::Anova() function) from grouped data.  As an example, you will find my beginnings of what I did in the attached R-file ("grouped anova from lm as data frame". The example data frame (auto.noise) I use for the example comes from the emmeans-package.
As an output, I just want something like this (Anova table grouped by size):
size Dependent HypothesisType Source DF SS MS FValue ProbF
L noise 1 side 1 833.3333333 833.3333333 50 0.000104954
L noise 1 type 1 75 75 4.5 0.066688
L noise 1 side*type 1 133.3333333 133.3333333 8 0.022203904
M noise 1 side 1 52.08333333 52.08333333 2.777777778 0.134140641
M noise 1 type 1 1752.083333 1752.083333 93.44444444 1.09266E-05
M noise 1 side*type 1 52.08333333 52.08333333 2.777777778 0.134140641
S noise 1 side 1 408.3333333 408.3333333 49 0.000112639
S noise 1 type 1 33.33333333 33.33333333 4 0.080516238
S noise 1 side*type 1 133.3333333 133.3333333 16 0.003949773

The output is generated with SAS (code below) using the 'by' statement within the glm procedure. (The model is probably not appropriate since there are repeated measures (factor side, however no subject is given in the data), but it works as a working example.  The example code comes from two sources:
1 https://dplyr.tidyverse.org/reference/do.html
2 https://cran.r-project.org/web/packages/broom/vignettes/broom_and_dplyr.html.
I tried to edit the code in order to get a data frame like the one above but without any success. At the end, it seems that I am too stupid to get the concept in my head. I hope that with your help I will better understand the concept and that I then can extract further information as a data frame like post hoc tukey tests and corresponding clds (compact letter display) and apply it to more complex models of lme-type.

I would be very happy if there is somebody who could help.
Thank you very much in advance
Regards
Christian
PS: here are the few lines that generates the above shown data set
proc glm data=auto_noise;
by size;
class side type;
model noise=side|type;
ods output ModelANOVA=ModelANOVA;
run;

And here comes the r-code as plain text:
# working example for grouped anova from lm as data frame
rm(list=ls())
library(emmeans)
library(broom)
library(dplyr)

str(auto.noise)
View(auto.noise)

# 1. I tried this corresponding to https://dplyr.tidyverse.org/reference/do.html

by_size <- auto.noise %>% group_by(size)
# do() with named arguments becomes nest_by() + mutate() & list()
models <- by_size %>% do(mod = lm(noise ~ type*side, data = .))
# ->
models %>% summarise(rsq = summary(mod)$r.squared)
# this works nicely but than I can't extract the anovatables =>

models %>% summary.aov(mod)
models %>% summarise(summary.aov(mod))



# 2. corresponding to  https://cran.r-project.org/web/packages/broom/vignettes/broom_and_dplyr.html
R.lmbysize<-auto.noise %>%
  nest(-size) %>%
  mutate(
    fit = map(data, ~ lm(noise ~ type*side, data = .x)),
    tidied = map(fit, tidy),
    glanced = map(fit, glance),
    augmented = map(fit, augment))

fit.bysize.df<-R.lmbysize%>%
  unnest(tidied)
View(fit.bysize.df)


stats.bysize.df<-R.lmbysize%>%
  unnest(glanced)
View(stats.bysize.df)


resid.pred.bysize.df<-R.lmbysize%>%
  unnest(augmented)
View(resid.pred.bysize.df)
View(R.lmbysize)

## how to get an extracted anova table by size? the following does not work
Anova.bysize.df<-R.lmbysize%>%
  unnest(summary.aov(fit))
View(fit.bysize.df)

summary.aov(fit.bysize.df)


# I tried to implement some thing like this: but it doesent work as well
Anova.tab=map(Anova.tab(fit),tidy)


******************************************
Dr. Christian Kluth
Georg-August-Universität Göttingen
Lehrkraft für besondere Aufgaben
Statistikberatung für Studierende
Department für Nutzpflanzenwissenschaften Carl-Sprengel-Weg 1
37075 Göttingen
Raum: 1.206
https://www.geodata.uni-goettingen.de/lageplan/?ident=4400_2_1.OG_1.206

Tel.: 0551/39-25582
E-Mail: mailto:[hidden email]


Termine nach Vereinbarung
*****************************************

______________________________________________
[hidden email] mailing list -- To UNSUBSCRIBE and more, see
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.