Preliminaries

#library(car)
library(openintro)
library(broom)
library(mosaic)
library(Stat2Data)
data("Diamonds2")

Introduction

We’ll use the Diamonds2 data set as an example. It is a subset of 307 cases with the most frequent colors from the Diamonds data. We want to know if the color has a significant impact on the weight of the diamond.

You should use Carat as response and Color as predictor.

str(Diamonds2)
'data.frame':   307 obs. of  6 variables:
 $ Carat     : num  1.08 0.31 0.32 0.33 0.33 0.35 0.35 0.37 0.38 0.38 ...
 $ Color     : Factor w/ 4 levels "D","E","F","G": 2 3 3 1 4 3 3 3 1 2 ...
 $ Clarity   : Factor w/ 8 levels "IF","SI1","SI2",..: 5 7 7 1 7 5 5 7 1 8 ...
 $ Depth     : num  68.6 61.9 60.8 60.8 61.5 62.5 62.3 61.4 60 61.5 ...
 $ PricePerCt: num  6693 3159 3159 4759 2896 ...
 $ TotalPrice: num  7229 979 1011 1570 956 ...

lm analysis

Write the complete linear model:

Estimate and interpret all coefficients:

diamond_lm <- lm(Carat ~ Color, data = Diamonds2)
tidy(diamond_lm)

Contrasts

We can see a tabular representation of the choice made in coding the dummy regressors using the contrasts function:

contrasts(Diamonds2$Color)
  E F G
D 0 0 0
E 1 0 0
F 0 1 0
G 0 0 1

Check conditions

Calculate the ratio between the maximum sd and the minimum sd. What is your conclusion?

diamond_summrize <- Diamonds2 %>%
    group_by(Color) %>%
    summarise(mean = mean(Carat), sd = sd(Carat)) %>%
    arrange(sd)
`summarise()` ungrouping output (override with `.groups` argument)
diamond_summrize
diamond_summrize$sd[4]/diamond_summrize$sd[1]
[1] 2.073307

Draw a QQ plot. What is your conclusion?

diamond_aug <- augment(diamond_lm)
diamond_aug

Draw a side-by-side box plot. What is your conclusion?

ggplot(Diamonds2, aes(y=Carat, x=Color)) + geom_boxplot()

ggplot(diamond_aug, aes(sample = .resid)) +
  geom_qq()

Your conclusion:

The ratio between the maximum sd and the minimum sd is larger than 2 and the QQ plot shows that the distribution of the residuals is normal. So we need to transform the data to correct the skewness of the data.

Transformation

Draw an appropriate plot of Carat and then decide a transformation for Carat.

ggplot(Diamonds2, aes(x=log(Carat))) + geom_histogram()

Run lm over the transformed data

diamond_log_lm <- lm(log(Carat) ~ Color, data = Diamonds2)
tidy(diamond_log_lm)

Check conditions of transformed data.

Draw a residual plot.

diamond_log_aug <- augment(diamond_log_lm)
diamond_log_aug

Draw a side-by-side box plot.

ggplot(Diamonds2, aes(y=log(Carat), x=Color)) + geom_boxplot()

Draw a QQ plot.

ggplot(diamond_log_aug, aes(sample = .std.resid)) +
  geom_qq()

Calculate the ratio between the maximum sd and the minimum sd. What is your conclusion?

diamond_log_summrize <- Diamonds2 %>%
    group_by(Color) %>%
    summarise(mean = mean(log(Carat)), sd = sd(log(Carat))) %>%
    arrange(sd)
`summarise()` ungrouping output (override with `.groups` argument)
diamond_log_summrize
diamond_log_summrize$sd[4]/diamond_log_summrize$sd[1]
[1] 1.559727

Run ANOVA over the transformed data

diamond_log_aov <- aov(log(Carat) ~ Color, data = Diamonds2)
tidy(diamond_log_aov)

Your conclusion:

Given this very small p-value, we reject the null hyphothesis, which means that there is at least one color is significantly different from other colors.

Run pairwise comparision using TukeyHSD.

diamond_aov <- aov(log(Carat) ~ Color, data = Diamonds2)
TukeyHSD(diamond_aov)
  Tukey multiple comparisons of means
    95% family-wise confidence level

Fit: aov(formula = log(Carat) ~ Color, data = Diamonds2)

$Color
           diff          lwr       upr     p adj
E-D -0.02298611 -0.227423474 0.1814513 0.9914530
F-D  0.20781627  0.005671473 0.4099611 0.0412729
G-D  0.35757624  0.154992248 0.5601602 0.0000439
F-E  0.23080238  0.053304460 0.4083003 0.0048645
G-E  0.38056235  0.202564412 0.5585603 0.0000004
G-F  0.14975996 -0.025600091 0.3251200 0.1238450
LS0tCnRpdGxlOiAiREFUQSAzNTAtIE9uZS13YXkgQU5PVkEgd2l0aCBUcmFuc2Zvcm1hdGlvbiAtICBQcmFjdGljZSBBbnN3ZXIiCm91dHB1dDoKICAgIGh0bWxfbm90ZWJvb2s6CiAgICAgICAgdG9jOiB5ZXMKICAgICAgICB0b2NfZmxvYXQ6IHllcwotLS0KCiMjIFByZWxpbWluYXJpZXMKCmBgYHtyLCBtZXNzYWdlID0gRkFMU0V9CiNsaWJyYXJ5KGNhcikKbGlicmFyeShvcGVuaW50cm8pCmxpYnJhcnkoYnJvb20pCmxpYnJhcnkobW9zYWljKQpsaWJyYXJ5KFN0YXQyRGF0YSkKZGF0YSgiRGlhbW9uZHMyIikKYGBgCgoKIyMgSW50cm9kdWN0aW9uCgpXZSdsbCB1c2UgdGhlIGBEaWFtb25kczJgIGRhdGEgc2V0IGFzIGFuIGV4YW1wbGUuIEl0IGlzIGEgc3Vic2V0IG9mIDMwNyBjYXNlcyB3aXRoIHRoZSBtb3N0IGZyZXF1ZW50IGNvbG9ycyBmcm9tIHRoZSBEaWFtb25kcyBkYXRhLiBXZSB3YW50IHRvIGtub3cgaWYgdGhlIGNvbG9yIGhhcyBhIHNpZ25pZmljYW50IGltcGFjdCBvbiB0aGUgd2VpZ2h0IG9mIHRoZSBkaWFtb25kLiAKCllvdSBzaG91bGQgdXNlIENhcmF0IGFzIHJlc3BvbnNlIGFuZCBDb2xvciBhcyBwcmVkaWN0b3IuIAoKYGBge3J9CnN0cihEaWFtb25kczIpCmBgYAoKIyMgbG0gYW5hbHlzaXMKCldyaXRlIHRoZSBjb21wbGV0ZSBsaW5lYXIgbW9kZWw6CgoKRXN0aW1hdGUgYW5kIGludGVycHJldCBhbGwgY29lZmZpY2llbnRzOgoKYGBge3J9CmRpYW1vbmRfbG0gPC0gbG0oQ2FyYXQgfiBDb2xvciwgZGF0YSA9IERpYW1vbmRzMikKdGlkeShkaWFtb25kX2xtKQpgYGAKCgojIyBDb250cmFzdHMKCldlIGNhbiBzZWUgYSB0YWJ1bGFyIHJlcHJlc2VudGF0aW9uIG9mIHRoZSBjaG9pY2UgbWFkZSBpbiBjb2RpbmcgdGhlIGR1bW15IHJlZ3Jlc3NvcnMgdXNpbmcgdGhlIGBjb250cmFzdHNgIGZ1bmN0aW9uOgoKYGBge3J9CmNvbnRyYXN0cyhEaWFtb25kczIkQ29sb3IpCmBgYAoKIyMgQ2hlY2sgY29uZGl0aW9ucwoKQ2FsY3VsYXRlIHRoZSByYXRpbyBiZXR3ZWVuIHRoZSBtYXhpbXVtIHNkIGFuZCB0aGUgbWluaW11bSBzZC4gV2hhdCBpcyB5b3VyIGNvbmNsdXNpb24/CgpgYGB7cn0KZGlhbW9uZF9zdW1tcml6ZSA8LSBEaWFtb25kczIgJT4lCiAgICBncm91cF9ieShDb2xvcikgJT4lCiAgICBzdW1tYXJpc2UobWVhbiA9IG1lYW4oQ2FyYXQpLCBzZCA9IHNkKENhcmF0KSkgJT4lCiAgICBhcnJhbmdlKHNkKQpkaWFtb25kX3N1bW1yaXplCmBgYAoKYGBge3J9CmRpYW1vbmRfc3VtbXJpemUkc2RbNF0vZGlhbW9uZF9zdW1tcml6ZSRzZFsxXQpgYGAKCkRyYXcgYSBRUSBwbG90LiBXaGF0IGlzIHlvdXIgY29uY2x1c2lvbj8KCmBgYHtyfQpkaWFtb25kX2F1ZyA8LSBhdWdtZW50KGRpYW1vbmRfbG0pCmRpYW1vbmRfYXVnCmBgYAoKRHJhdyBhIHNpZGUtYnktc2lkZSBib3ggcGxvdC4gV2hhdCBpcyB5b3VyIGNvbmNsdXNpb24/IAoKYGBge3J9CmdncGxvdChEaWFtb25kczIsIGFlcyh5PUNhcmF0LCB4PUNvbG9yKSkgKyBnZW9tX2JveHBsb3QoKQpgYGAKCmBgYHtyfQpnZ3Bsb3QoZGlhbW9uZF9hdWcsIGFlcyhzYW1wbGUgPSAucmVzaWQpKSArCiAgZ2VvbV9xcSgpCmBgYAoKWW91ciBjb25jbHVzaW9uOgoKCjxkaXYgY2xhc3MgPSAiYW5zd2VyIj4KClRoZSByYXRpbyBiZXR3ZWVuIHRoZSBtYXhpbXVtIHNkIGFuZCB0aGUgbWluaW11bSBzZCBpcyBsYXJnZXIgdGhhbiAyIGFuZCB0aGUgUVEgcGxvdCBzaG93cyB0aGF0IHRoZSBkaXN0cmlidXRpb24gb2YgdGhlIHJlc2lkdWFscyBpcyBub3JtYWwuIFNvIHdlIG5lZWQgdG8gdHJhbnNmb3JtIHRoZSBkYXRhIHRvIGNvcnJlY3QgdGhlIHNrZXduZXNzIG9mIHRoZSBkYXRhLiAKCjwvZGl2PgoKIyMgVHJhbnNmb3JtYXRpb24KCkRyYXcgYW4gYXBwcm9wcmlhdGUgcGxvdCBvZiBDYXJhdCBhbmQgdGhlbiBkZWNpZGUgYSB0cmFuc2Zvcm1hdGlvbiBmb3IgQ2FyYXQuIAoKYGBge3J9CmdncGxvdChEaWFtb25kczIsIGFlcyh4PWxvZyhDYXJhdCkpKSArIGdlb21faGlzdG9ncmFtKCkKYGBgCgojIyBSdW4gbG0gb3ZlciB0aGUgdHJhbnNmb3JtZWQgZGF0YQoKYGBge3J9CmRpYW1vbmRfbG9nX2xtIDwtIGxtKGxvZyhDYXJhdCkgfiBDb2xvciwgZGF0YSA9IERpYW1vbmRzMikKdGlkeShkaWFtb25kX2xvZ19sbSkKYGBgCgojIyBDaGVjayBjb25kaXRpb25zIG9mIHRyYW5zZm9ybWVkIGRhdGEuCgpEcmF3IGEgcmVzaWR1YWwgcGxvdC4KCmBgYHtyfQpkaWFtb25kX2xvZ19hdWcgPC0gYXVnbWVudChkaWFtb25kX2xvZ19sbSkKZGlhbW9uZF9sb2dfYXVnCmBgYAoKCkRyYXcgYSBzaWRlLWJ5LXNpZGUgYm94IHBsb3QuCgpgYGB7cn0KZ2dwbG90KERpYW1vbmRzMiwgYWVzKHk9bG9nKENhcmF0KSwgeD1Db2xvcikpICsgZ2VvbV9ib3hwbG90KCkKYGBgCgpEcmF3IGEgUVEgcGxvdC4KCmBgYHtyfQpnZ3Bsb3QoZGlhbW9uZF9sb2dfYXVnLCBhZXMoc2FtcGxlID0gLnN0ZC5yZXNpZCkpICsKICBnZW9tX3FxKCkKYGBgCgoKQ2FsY3VsYXRlIHRoZSByYXRpbyBiZXR3ZWVuIHRoZSBtYXhpbXVtIHNkIGFuZCB0aGUgbWluaW11bSBzZC4gV2hhdCBpcyB5b3VyIGNvbmNsdXNpb24/CgpgYGB7cn0KZGlhbW9uZF9sb2dfc3VtbXJpemUgPC0gRGlhbW9uZHMyICU+JQogICAgZ3JvdXBfYnkoQ29sb3IpICU+JQogICAgc3VtbWFyaXNlKG1lYW4gPSBtZWFuKGxvZyhDYXJhdCkpLCBzZCA9IHNkKGxvZyhDYXJhdCkpKSAlPiUKICAgIGFycmFuZ2Uoc2QpCmRpYW1vbmRfbG9nX3N1bW1yaXplCmBgYAoKYGBge3J9CmRpYW1vbmRfbG9nX3N1bW1yaXplJHNkWzRdL2RpYW1vbmRfbG9nX3N1bW1yaXplJHNkWzFdCmBgYAoKIyMgUnVuIEFOT1ZBIG92ZXIgdGhlIHRyYW5zZm9ybWVkIGRhdGEKCmBgYHtyfQpkaWFtb25kX2xvZ19hb3YgPC0gYW92KGxvZyhDYXJhdCkgfiBDb2xvciwgZGF0YSA9IERpYW1vbmRzMikKdGlkeShkaWFtb25kX2xvZ19hb3YpCmBgYAoKWW91ciBjb25jbHVzaW9uOgoKCjxkaXYgY2xhc3MgPSAiYW5zd2VyIj4KCkdpdmVuIHRoaXMgdmVyeSBzbWFsbCBwLXZhbHVlLCB3ZSByZWplY3QgdGhlIG51bGwgaHlwaG90aGVzaXMsIHdoaWNoIG1lYW5zIHRoYXQgdGhlcmUgaXMgYXQgbGVhc3Qgb25lIGNvbG9yIGlzIHNpZ25pZmljYW50bHkgZGlmZmVyZW50IGZyb20gb3RoZXIgY29sb3JzLiAKCjwvZGl2PgoKIyMgUnVuIHBhaXJ3aXNlIGNvbXBhcmlzaW9uIHVzaW5nIFR1a2V5SFNELgoKYGBge3J9CmRpYW1vbmRfYW92IDwtIGFvdihsb2coQ2FyYXQpIH4gQ29sb3IsIGRhdGEgPSBEaWFtb25kczIpClR1a2V5SFNEKGRpYW1vbmRfYW92KQpgYGAK