#
options(digits=5, show.signif.stars=FALSE)
# install.packages("pacman")
pacman::p_load(alr4, tidyverse, magrittr)
# make sure data set is active
data(Stevens, package="alr4")
#
ot <- theme_set(theme_bw())
# 1 regression fits all
ggplot(data=Stevens, aes(loudness, y, group = subject, alpha = subject))+
geom_point() +
geom_line() +
stat_smooth(aes(group = 1), method = "lm", size = rel(.8)) +
coord_cartesian(ylim = c(-1, 5)) +
labs(x = "Loudness (db)", y = "Average log(Length)") +
theme(legend.position = "NONE")
## Warning: Using alpha for a discrete variable is not advised.
## `geom_smooth()` using formula 'y ~ x'
#m0 output
summary(m0 <- lm(y~loudness, data = Stevens))
##
## Call:
## lm(formula = y ~ loudness, data = Stevens)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.2669 -0.3179 -0.0093 0.3707 0.9219
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -3.20377 0.34429 -9.31 2.5e-12
## loudness 0.07901 0.00482 16.39 < 2e-16
##
## Residual standard error: 0.482 on 48 degrees of freedom
## Multiple R-squared: 0.848, Adjusted R-squared: 0.845
## F-statistic: 269 on 1 and 48 DF, p-value: <2e-16
# 1 regression model per individual
ggplot(data=Stevens, aes(loudness, y, group = subject, color = subject))+
geom_point() +
stat_smooth(method = "lm", se = F, size = rel(.5)) +
scale_color_grey()+
coord_cartesian(ylim = c(-1, 5)) +
labs(x = "Loudness (db)", y = "Average log(Length)") +
theme(legend.position = "NONE")
## `geom_smooth()` using formula 'y ~ x'
# grouped by subject before doing regression
Stevens %>%
group_by(subject) %>%
do(broom::tidy(lm(y ~ loudness, data = Stevens)))
## # A tibble: 20 x 6
## # Groups: subject [10]
## subject term estimate std.error statistic p.value
## <fct> <chr> <dbl> <dbl> <dbl> <dbl>
## 1 A (Intercept) -3.20 0.344 -9.31 2.53e-12
## 2 A loudness 0.0790 0.00482 16.4 2.69e-21
## 3 B (Intercept) -3.20 0.344 -9.31 2.53e-12
## 4 B loudness 0.0790 0.00482 16.4 2.69e-21
## 5 C (Intercept) -3.20 0.344 -9.31 2.53e-12
## 6 C loudness 0.0790 0.00482 16.4 2.69e-21
## 7 D (Intercept) -3.20 0.344 -9.31 2.53e-12
## 8 D loudness 0.0790 0.00482 16.4 2.69e-21
## 9 E (Intercept) -3.20 0.344 -9.31 2.53e-12
## 10 E loudness 0.0790 0.00482 16.4 2.69e-21
## 11 F (Intercept) -3.20 0.344 -9.31 2.53e-12
## 12 F loudness 0.0790 0.00482 16.4 2.69e-21
## 13 G (Intercept) -3.20 0.344 -9.31 2.53e-12
## 14 G loudness 0.0790 0.00482 16.4 2.69e-21
## 15 H (Intercept) -3.20 0.344 -9.31 2.53e-12
## 16 H loudness 0.0790 0.00482 16.4 2.69e-21
## 17 I (Intercept) -3.20 0.344 -9.31 2.53e-12
## 18 I loudness 0.0790 0.00482 16.4 2.69e-21
## 19 J (Intercept) -3.20 0.344 -9.31 2.53e-12
## 20 J loudness 0.0790 0.00482 16.4 2.69e-21
# fitting a model treating subject as fixed effects
summary(m1 <- lm(y~ subject/loudness - 1, data = Stevens))
##
## Call:
## lm(formula = y ~ subject/loudness - 1, data = Stevens)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.3318 -0.1073 0.0023 0.1485 0.3242
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## subjectA -4.49370 0.48664 -9.23 2.8e-10
## subjectB -4.58530 0.48664 -9.42 1.8e-10
## subjectC -2.29360 0.48664 -4.71 5.2e-05
## subjectD -2.91940 0.48664 -6.00 1.4e-06
## subjectE -5.38510 0.48664 -11.07 4.1e-12
## subjectF -2.42120 0.48664 -4.98 2.5e-05
## subjectG -2.81090 0.48664 -5.78 2.6e-06
## subjectH -3.42070 0.48664 -7.03 8.2e-08
## subjectI -1.94870 0.48664 -4.00 0.00038
## subjectJ -1.75910 0.48664 -3.61 0.00109
## subjectA:loudness 0.10265 0.00681 15.06 1.6e-15
## subjectB:loudness 0.09355 0.00681 13.73 1.8e-14
## subjectC:loudness 0.06950 0.00681 10.20 2.9e-11
## subjectD:loudness 0.07506 0.00681 11.02 4.6e-12
## subjectE:loudness 0.10391 0.00681 15.25 1.1e-15
## subjectF:loudness 0.07740 0.00681 11.36 2.2e-12
## subjectG:loudness 0.06497 0.00681 9.53 1.4e-10
## subjectH:loudness 0.08223 0.00681 12.07 4.9e-13
## subjectI:loudness 0.06561 0.00681 9.63 1.1e-10
## subjectJ:loudness 0.05525 0.00681 8.11 4.7e-09
##
## Residual standard error: 0.215 on 30 degrees of freedom
## Multiple R-squared: 0.996, Adjusted R-squared: 0.993
## F-statistic: 369 on 20 and 30 DF, p-value: <2e-16
knitr::spin("C:\\Users\\ASUS\\Desktop\\loudness.R", knit=FALSE)
## [1] "C:\\Users\\ASUS\\Desktop\\loudness.Rmd"