Analysis of Variance 方差分析
library(dplyr)
library(ggplot2)
library(knitr)
library(palmerpenguins)
library(tidyr)
library(ISLR2)
Analysis of Variance 方差分析
Assumptions
- 每个样本服从正态分布。
- 每个样本方差\(\sigma^2\)相同。
- 样本之间相互独立 levels：a factor has many levels.
要做的事情：比较各总体（每个水平下的总体）的均值是否相等
One Way ANOVA
举例：比较三种疗法的好坏 如果是两种，可以使用t test
这里，可以把三种疗法看作3个因子，（因子变量/分类变量），响应变量：治疗效果
Treatment is a between-groups factor with 3 levels.
This model is called one-way ANOVA, bacause there's only one
classification variable.
So this is a one-way between-groups ANOVA.
这里可以看作是独立的，对于其中两组可以做独立T检验。
对于某一种疗法，下面引入一个新的变量------时间。例如：在5周之后做测量，在6个月之后做测量。
Time is within-groups factor with 2 levels.
This model is called one-way within-groups ANOVA. （因为疗法都是一样的）
这里不是独立的，对于其中两组做T检验时，paired=TRUE
如果同时考虑疗法和时间： This model is a two-factor ANOVA.
对于两个因子，可以考虑它们时间是否有交互作用。
Eg. 5种药物
library(tidyverse)
data(cholesterol, package = "multcomp")
plotdata <- cholesterol %>%
group_by(trt) %>%
summarize(n = n(), #n()表示出现系数。
mean = mean(response),
sd = sd(response))
plotdata
ggplot(plotdata, aes(x = trt, y = mean, group = 1)) +
geom_point(size = 3, color = "red") +
geom_line(linetype = "dashed", color = "darkgrey") +
geom_errorbar(aes(ymin = mean - qt(0.975, df = n - 1) * sd / sqrt(n),
ymax = mean + qt(0.975, df = n - 1) * sd / sqrt(n)), width = .1) +
labs(x = "Treatment", y = "Response", title = "Mean plot with 95% confidence interval")
# geom_errorbar()是在算95%的置信区间
One-way ANOVA
fit <- aov(response ~ trt, data = cholesterol)
summary(fit)
原假设：所有方法都一样 备择假设：至少有一个方法不一样
所以此处可以拒绝原假设
于是，我们对这样一个模型做linear regression：
当Y属于第j类时，\(X_j\) = 1，否则\(X_j\) = 0
注意在这里\(\mu\)可以直接用总体来估计。 并且\(\sum{\alpha_{ij}} = 0\)
\(\hat b = ( X^T X) X^T Y\)
mean(cholesterol$response)
fit2 <- lm(response - mean(response) ~ 0 + trt, data = cholesterol) #这里的结果表示
fit2 <- lm(response ~ 0 + trt, data = cholesterol) #这里的结果表示：trt1time的效果比平均效果差6.95
fit2 <- lm(response ~ trt, data = cholesterol) #这里的结果表示：截距就是trt1time的疗效。trt1time的结果是5.78，trt2time的结果比trt2time的高3.443
fit2
要验证\(H_0: \alpha_1 = \alpha_2 = \alpha_3 = \alpha_4 = \alpha_5 = 0\)
\(H_1\) = 至少有一个\(\alpha_j\)不为0。
这里可以选择验证\(\sum {\alpha_i^2} = 0\)
它服从卡方分布。 d.f. = K - 1
这里方差未知，需要估计$$\hat \sigma^2 = s^2 = \frac{1}{n-K} \sum \sum (Y_{ij} - \bar Y_{ij})^2 $$它服从卡方分布。
d.f. = n - K
用这两个卡方分布去构造一个F分布的统计量。
说白了就是一开始就想知道$$\frac{1}{\sigma^2}\sum_{j = 1}^{K} n_j (\bar Y.,j - \bar Y ..)^2$$是不是接近0，所以去找了它的分布，然后发现它是一个卡方分布，然后就可以构造求p值
$$$$
我们发现在$$Y = \mu + \alpha_j + \epsilon_{ij}$$
中前两项是已知的，最后一项可以生成，所以我们可以生成Y，
\(sdf\)
Eg. two-way factorial ANOVA
data(ToothGrowth)
ToothGrowth$dose <- factor(ToothGrowth$dose) # what happens if we remove this line?
stats <- ToothGrowth %>%
group_by(supp, dose) %>%
summarize(n = n(),
mean = mean(len),
sd = sd(len),
std.error = sd/sqrt(n))
stats
fit1 <- aov(len ~ supp * dose, data = ToothGrowth)
summary(fit1)
fit2 <- aov(len ~ supp + dose, data = ToothGrowth)
summary(fit2)
注意：把dose看作double和把dose看作一个因子，得到的是有差异的。
单因子H_0：所有的alpha都是0 多因子H_0：所有的alpha都是0，
\(y\) supp: \alpha_1 -> OJ, \alpha_2 -> VC dose: \beta_1 -> 0.5,
\beta_2 -> 1 \beta_3 -> 1
\(y_{ijk}\), i-th supp, j-th dose, k-th sample
\(y_{ijk} = \mu + \alpha_i + \beta_j + \epsilon_{ijk}\)
\(\mu + \alpha_1\) 表示OJ维生素的样本的均值
同样假设 \(\epsilon_{ijk} \sim N(0, \sigma^2)\)
\(n_{ij} = 10, i \in (0,2], j \in (0,3]\) the number of samples that take
i-th vitamin and j-th dose
\(H_0: \alpha_1 = \alpha_2 = 0\) \(\alpha_1^2 + \alpha_2^2\) 小即可
即\(\sum_{i=1}^2 (y_i - \bar y)^2\) 就是上面表里的205，sum of square, d.f.
= 2-1 = 1 (\(\alpha_1 + \alpha_2 = 0\)) Residual sum of square:
\(\sum (y_{ijk} - \bar y_{ij.})^2}\)
total sum of square \(\sum_{ijk}(y_{ijk} - \bar y_{...})^2\) d.f. = 60 - 1
= 59，最后residual的那个d.f. = 59 - 1 - 2 = 56.
# validate 交互效应
#y ~ supp + dose + supp:dose
#也就是乘法
同样的剂量和不同的维生素之间是有差异的，
如果没有交互效应，两条线应该是平行的。
data(ToothGrowth)
lm(len ~ supp + dose, data = ToothGrowth)
ToothGrowth$dose = as.double(ToothGrowth$dose)
lm(len ~ supp + dose, data = ToothGrowth)
intercept表示当不用suppVC也不用dose就是9.272
suppVC表示，使用suppVC效果会差3.7
dose表示，固定其他的东西的情况下，多使用一个单位的剂量，效果会好9，764
这里因为factor类型的只有一个supp, 所有条件 if OJ y = 9.272 + 9.764dose
if VC y = 9.272 - 3.700 + 9.764dose
注意如果当dose是分类变量的时候会不一样。会出现dose1，dose2，相应的所有变量的解释会变得不一样。
同时，d.f.也会不一样。
summary： perdictor-response 1. quan-quan :regression 2. cate-quan :
anova 3. mix-quan: anova
Classification
when response is categorical.
\(E(Y|X) = \beta_0 + \beta_1 X\) \(Var(Y|X) = \sigma^2\)
P(Y = 1|X) = \(\beta_0 + \beta_1 X\) P(Y = 0|X)
用线性模型来拟合概率，但是线性模型求出的概率可能为负数。
这时，我们使用logit。 利用logistic function\(y = \frac{e^x}{1+e^x}\)
此时\(P(Y = 1|X) = \frac{e^{\beta_0+\beta_1X}}{1 + e^{\beta_0 + \beta_1 X}}\)
Eg. Default dataset
data(Default)
b <- ggplot(Default[1:1000,],aes(x=balance,y=income))
b + geom_point(aes(color=student))
Eg. IRIS dataset
b <- ggplot(iris,aes(x=Sepal.Length,y=Sepal.Width))
b + geom_point(aes(color=Species))
response variable Y = 1,2,3 它们之间实际上没有顺序关系，不能用线性回归，
Logistic regression
odds = $\frac{p(X)}{1 - p(X)} $
表示划分为一类的概率和划分为另一类的概率之比。
X每增加一个单位，odds就增加\(e^{\beta_1}\)倍。
这里既然涉及概率，用Maximum Likelihood来求解。
联合概率是：\(P(Y_1 = y_1,... Y_n = y_n | X_1,...,X_n) = \Pi_{i=1}^n P(Y_i = y_i | X_i = x_i)\)
其中 if \(y_i = 0\)
\(P(Y_i = y_i | X_i = x_i) = P(Y_i = 0 | X_i = x_i) = 1 - P(Y_i = 1 | X_i = x_i) = 1 - \frac{e^{\beta_0+\beta_1 x_i}}{1+e^{\beta_0+\beta_1 x_i}} = \frac{1}{1+e^{\beta_0+\beta_1 x_i}}\)
if \(y_i = 1\)
\(P(Y_i = y_i | X_i = x_i) = P(Y_i = 0 | X_i = x_i) = 1 - P(Y_i = 1 | X_i = x_i) = 1 - \frac{e^{\beta_0+\beta_1 x_i}}{1+e^{\beta_0+\beta_1 x_i}} = \frac{1}{1+e^{\beta_0+\beta_1 x_i}}\)
没有解析解，只能求数值解，
set.seed(1)
n <- nrow(Default)
relabel <- sample(n,n)
train.index <- relabel[1:9000]
test.index <- relabel[9001:10000]
train.data <- Default[train.index,]
glm.fit <- glm(student~income,data=Default,family="binomial")
glm.fit$coefficients
如果family = "gaussian"就是linear regression
Hypothesis tests
\(H_0: \beta_1 = 0\)
z-score = \(\frac{\hat \beta_1}{SE(\hat \beta_1)}\)
还是想知道\(\beta_1\)的分布，但是只能知道它的渐进分布为正态分布。
summary(glm.fit)
x.new=data.frame(c(10000,40000))
colnames(x.new) = 'income' #要把列名强调好
predict(glm.fit,newdata = x.new,type = "response")
test.data <- Default[test.index,]
glm.prob <- predict(glm.fit,test.data,type="response")
glm.pred <- rep("No",1000)
glm.pred[glm.prob > 0.5] <- "Yes"
table(glm.pred,test.data$student)
mean(glm.pred != test.data$student)
Multiple Logistic Regression
\(\beta_1\)表示固定其他的不变，在X_1增加一个单位时，左边会增大\(e^\{\beta_1}\)倍。
glm.fit2 <- glm(student~balance+income,data=train.data,family='binomial')
glm.fit2$coefficients
summary(glm.fit2)$coef
test.data <- Default[test.index,]
glm.prob2 <- predict(glm.fit2,test.data,type="response")
glm.pred2 <- rep("No",1000)
glm.pred2[glm.prob2 > 0.5] <- "Yes"
table(glm.pred2,test.data$student)
mean(glm.pred2 != test.data$student)
发现增加了一个predictor效果还变差了。需要找个一个评估模型的方法。例如：cross-validation
Multinomial Logistic Regression
library(nnet)
multinom.fit <- multinom(Species ~ ., data=iris,trace = FALSE)
summary(multinom.fit)$coefficients
table(iris$Species,predict(multinom.fit))
