linear_reg.Rmd

---
title: "linear_reg"
author: "liuc"
date: '2022-06-01'
output: pdf_document
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```


## linear regression

线性回归作为应用很多，注意事项也很多的算法，在医学研究中有很多应用，此处做一个详细的示例。
主要关注点在于 计算标准化回归系数。
同时还会加上：Nomogram，森林图，校准曲线，DCA曲线等预测模型常见的内容。

回归模型的一大特点便是其具有a one-unit change in x has the same effect on Y regardless of what the other terms in the model are.
contrast plots和conditional plots的异同点在于？conditional plot是在考虑其他所有term的情况下展示某一term的方法， 一般其他term set to their median if the term is numberic or the most common category if the term is a factor. contrast plot 是对单位变化的绘制。

所用数据为`lasso_model.Rmd`中的office数据集。



### 标准化回归系数

标准化回归系数，是指消除了因变量和自变量所取单位的影响之后的回归系数，其绝对值的大小直接反映了自变量对因变量的影响程度。标准化回归系数的比较结果只是适用于某一特定环境的，而不是绝对正确的，它可能因时因地而变化。
通常我们在构建多因素回归模型时，方程中呈现的是未标准化回归系数，它是方程中不同自变量对应的原始的回归系数。它反映了在其他因素不变的情况下，该自变量每变化一个单位对因变量的作用大小。通过未标准化回归系数和常数项构建的方程，便可以对因变量进行预测，并得出结论。
而对于标准化回归系数，它是在对自变量和因变量同时进行标准化处理后所得到的回归系数，数据经过标准化处理后消除了量纲、数量级等差异的影响，使得不同变量之间具有可比性，因此可以用标准化回归系数来比较不同自变量对因变量的作用大小。
*通常我们主要关注的是标准化回归系数的绝对值大小，绝对值越大，可认为它对因变量的影响就越大。*


```{r, include=FALSE}
library(car)
library(lmtest)
library(MASS)
library(glmnet)
library(rms)
library(sjPlot) # 绘制标准化回归系数
# library(QuantPsyc) # 计算标准化回归系数
library(visreg)
library(regplot) # A function to plot a regression nomogram of regression objects
library(finalfit) #
```


#### 一个lm的示例

最小二乘法，最小实际response变量和预测response值间的difference，即minimize the sum of squared residuals。OLS also assumes that the independent variables are fixed and measured without error, but this assumption is typically relaxed in practice. 既是讲OLS将自变量Predictor视为一个固定值，通常在GLM中自变量被视为随机变量。

除了`lm()`外，R中还提供了一系列的函数帮助提取lm对象的内容：
- `summary()` Displays detailed results for the fitted model
- `coefficients()` Lists the model parameters (intercept and slopes) for the fitted model
- `confint()` Provides confidence intervals for the model parameters (95% by default)
- `fitted()` Lists the predicted values in a fitted model
- `residuals()` Lists the residual values in a fitted model
- `anova()` Generates an ANOVA table for a fitted model, or an ANOVA table comparing two or more fitted models
- `vcov()` Lists the covariance matrix for model parameters
- `AIC()` Prints Akaike’s Information Criterion
- `plot()` Generates diagnostic plots for evaluating the fit of a model 
- `predict()` Uses a fitted model to predict response values for a new dataset


##### 简单线性回归

如下所示的简单线性回归，其在实际的应用中等同于Pearson's product-moment correlation. 可以理解为简单线性回归所能提供的
数据分析结果为两变量间的相关性。

```{r}
# 
fit1 <- lm(weight ~ height, data=women)

summary(fit1)
```
*interpret: *The residual standard error (1.53 pounds) can be thought of as the average error in predicting weight from height using this model. 
在只有一个自变量的情况下，f-statistic值和t-test的检验P值一致。

```{r}
# 和pearson相关性的相关系数进行对比
# 系数的平方=R2
cor(women$weight, fitted(fit1))^2
```


```{r}
# 得到每一个观测的预测值
broom::augment(fit1)
```


如何理解线性关系呢？可以通过散点图等观察X自变量和Y因变量间的关系，似乎没有某种计算方式验证是否线性。
```{r}
plot(women$height,women$weight, 
xlab="Height (in inches)", ylab="Weight (in pounds)")

abline(fit1)
```

##### 多项式回归

Fitting Polynomial Regression

多项式回归或者通过log等对自变量进行转化后的模型依旧是线性模型，似乎对非线性模型有所误解呀我。。或者不能说是线性关系，因为多项式是明显的非线形关系。

在对数据进行log处理的时候，一般选择对自变量进行log处理，如果对因变量进行log处理，则其和其他自变量的关系都会发生改变。

```{r}
# 这种公式的写法并不推荐
fit2 <- lm(weight ~ height + I(height^2), data=women)

# summary(fit2)

# an nth-degree polynomial produces a curve with n-1 bends
fit3 <- lm(weight ~ height + I(height^2) +I(height^3), data=women)

summary(fit3)
```

对于多项式回归，同时要考虑到变量间的共线性问题，所以一般采用`poly`函数的结果更为可靠，因其提供了正交多项式以消除共线性问题。

> https://stats.stackexchange.com/questions/95939/how-to-interpret-coefficients-from-a-polynomial-model-fit

对于临床上应用的多项式回归在进行线性检验时，有额外的考量，临床标准的线性检验中使用了两个统计量，ADL（偏离直线平均差异，average deviation from linearity）与PctBnd（百分区界 ，percent bound） ，对于大多数分析物PctBnd取5%。如ADL小于所要求的临界判断值，则可认为数据组具有临床可接受的线性，所拟合出的最适非线性多项式无临床意义。

```{r}
# 和上面的fit3的结果差别还是蛮大的
# 主要原因在于poly的应用上，poly(height, degree = 3, raw = TRUE) 将得到一致的结果
fit3_2 <- lm(weight ~ poly(height, degree = 3, raw = FALSE), data=women)

summary(fit3_2)
```

```{r}
ggplot(data = women, aes(x = weight, y = height)) +
  geom_point() +
  geom_smooth(method = 'lm', 
              formula = y ~ ploy(x, 2)
              )
```

```{r}
sjPlot::tab_model(fit3_2, show.se = T, show.stat = T, show.df = T, 
                  show.fstat = T, show.dev = T)
```

```{r}
modelsummary::modelsummary(fit3_2)

modelsummary::modelplot(fit3_2)
```

*interpret the Result: *
We can conclude from this plot that the slope is (significantly, because the confidence interval does not cover 0) negative until 0, and then becomes positive. 0 corresponds indeed to the point of inversion of the curve.


```{r}
# estime_slopes 提供对连续变量间的关系
deriv <- modelbased::estimate_slopes(fit3_2, trend = "x", at = "x")

plot(deriv) +
  geom_hline(yintercept = 0, linetype = "dashed")

summary(deriv)
```




##### multiple regression

多元回归，自变量不止一个的回归模型。

自变量间的共线性问题需要时时关注，同时在实际应用中，各种协变量、混淆变量、调节变量等等的分析总是很重要的。

虽则`先单后多`的因素筛选方法有值得借鉴的地方，但同时应该谨记这只是用于参考的方法，具体的变量筛选还是要结合多种方面进行选择。

```{r}
states <- as.data.frame(state.x77[,c("Murder", "Population", "Illiteracy", "Income", "Frost")])
```


`car`包所提供的scatter plot，方便整体检查变量间的关系。
```{r}
car::scatterplotMatrix(states, smooth=FALSE, main="Scatter Plot Matrix")
```

```{r}
fit4 <- lm(Murder ~ Population + Illiteracy + Income + Frost, data=states)

summary(fit4)
```
*interpret: *When there’s more than one predictor variable, the regression coefficients indicate the increase in the dependent variable for a unit change in a predictor variable, holding all other predictor variables constant. For example, the regression coefficient for Illiteracy is 4.14, suggesting that an increase of 1% in illiteracy is associated with a 4.14% increase in the murder rate, controlling for population, income, and temperature. 


##### 交互作用的示例

以一个只含有交互作用的简单例子进行解释，如果增加额外的变量要视变量类型进行分析。

```{r}
fit5 <- fit <- lm(mpg ~ hp + wt + hp:wt, data=mtcars)

summary(fit5)
```
*interpret results: *mpg = 49.81 – 0.12 × hp – 8.22 × wt + 0.03 × hp × wt
上述公式很好的说明了交互作用的考量。


```{r}
# list里面需要指定具体的数值，不过对于分类变量该如何呢
plot(effects::effect("hp:wt", fit5,, list(wt = c(2.2, 3.2, 4.2))),
  lines = c(1, 2, 3), multiline = TRUE
)

```

*进一步的交互作用interpret的理清：*
首先看待一个三个变量交互作用的情况。三个变量或多个变量间的交互在阐释时变的复杂，过多的交互作用，不论是从样本量还是实际意义来说都是要避免的，不过三个还算是常见的。
三个分类变量交互作用的一大特点就是，谨记在dummy编码的情况下，模型的分类变量的系数是可以随着变量level的不同而不存在的，即为0时，此项系数就为0了。最为极端的例子就是所有的分类变量取intercept中的level值时，所有的回归系数在方程中都为0.

一个简单的理解就是在不同的level组合上数据具有不同的表现。




##### 模型诊断（model diagnostics）

模型诊断在`base`R中通过`plot`函数提供了常见的四幅图。不过更为常用的还是基于`car`包中所提供的函数，显然做回归相关的事情是离不开此包的。

以下为`car`包中所提供的常用的模型诊断函数。

qqPlot()   Quantile comparisons plot
durbinWatsonTest()    Durbin–Watson test for autocorrelated errors
crPlots()    Component plus residual plots
ncvTest()    Score test for nonconstant error variance
spreadLevelPlot()     Spread-level plots
outlierTest()   Bonferroni outlier test
avPlots()      Added variable plots
influencePlot()    Regression influence plots
scatterplot()   Enhanced scatter plots
scatterplotMatrix() Enhanced scatter plot matrixes
vif()       Variance inflation factors

```{r}
# NORMALITY
car::qqPlot(fit5)

# INDEPENDENCE OF ERRORS
car::durbinWatsonTest(fit)

# LINEARITY
car::crPlot(fit)

# HOMOSCEDASTICITY
car::ncvTest(fit)
car::spreadLevelPlot(fit)

# Multicollinearity
car::vif(fit4)
car::vif(fit5, type ='predictor')
```


在展开回归分析时，同时应该考虑异常样本的情况，比如outlier等，有些样本的存在对回归系数的影响颇大，这不是一个好的方向，往往也会影响模型的泛化能力。

```{r}
car::outlierTest(fit4)
# Nevada 这一观测 为outlier样本

# 和check_outlier对比
performance::check_outliers(fit4)
```

此函数整合了outlier, leverage, and influence plots into one informative plot
```{r}
car::influencePlot(fit4, 
                   id=TRUE, 
                   main="Influence Plot",
                   sub="Circle size is proportional to Cook's distance")
```
*interpret results:* 此图超过|2|的为outlier，States above 0.2 or 0.3 on the horizontal axis have high leverage。


在`car`包所提供的丰富的信息中，还有提供在回归模型不满足条件时对Y变量的转换，包括删除outlier、transformation Y等。
什么时候可以对变量进行转换，什么时候不可以呢？比如在进行log转换时，结果的interpret要在转换后的基础上进行阐释。

```{r}
# 函数powerTransform提供的结果可供参考
# 其提供Box-Cox转换，即正态性转换
summary(f2 <- car::powerTransform(states$Murder))

```

```{r}
# 最佳变换
hist(bcPower(states$Murder, f2$lambda))

# 对数变换
hist(bcPower(states$Murder, f2$roundlam))

```

除了上述的正态性转换外，对不满足线形关系的还可以进行Box-Tidwell转换
```{r}
car::boxTidwell(Murder ~ Population, data=states)

```


*还可以利用lmtest包提供的函数进行模型的检验：*
```{r}
# 残差检验, 检验异方差性
# gqtest()：用于检验线性回归模型的异方差性。它可以使用Godfrey-Quinn检验和Newey-West检验等方法。
lmtest::gqtest(fit4)
lmtest::bptest(fit4)

# 用于检验线性回归模型系数的显著性。它可以使用t检验、F检验和Wald检验等方法。
coeftest(fit4, vcov = vcovHC(fit4)) # 异方差性稳健标准误差
```



##### 模型选择和对比

`anova`只能针对nested models进行？？？
Note that although the ANOVA approach requires nested models, the AIC approach doesn’t.

```{r}
fit6 <- lm(Murder ~ Population + Illiteracy, data=states)

anova(fit4, fit6)

AIC(fit4, fit6)
```

*stepwise regression:*
```{r}
stats::step(fit4, direction="backward")
```


*ALL SUBSETS REGRESSION*

In all subsets regression, every possible model is inspected. 在这一点上可以看出其比step方法要好一些。
All subsets regression is performed using the regsubsets() function from the leaps package.

```{r}
leaps <- leaps::regsubsets(Murder ~ Population + Illiteracy + Income + Frost, data=states, nbest=4)

leaps
```


```{r}
subsTable <- function(obj,scale){
  x <- summary(leaps)
  m <- cbind(round(x[[scale]],3), x$which[,-1]) 
  colnames(m)[1] <- scale
  m[order(m[,1]), ]
}

subsTable(leaps, scale="adjr2")
```


##### 模型泛化与预测

以预测为目的的模型和以interpret为目的模型在模型拟合上具有不同之处，OLS/ML等都是在尽量让数据符合模型而得到的参数，而预测模型同时还需要考虑其的泛化能力。

这些知识点在机器学习的部分已经有更为充分的函数和脚本，此处暂且不论。


*变量对outcome的重要性：*

除了标准化回归系数以外还有一些其他的方法：

Several possible approaches to relative importance are captured in the `relaimpo` package.

```{r}
vip::vip(fit4) +
  ggthemes::theme_economist()
```

```{r}

sjPlot::tab_model(fit4, show.std = 'std')
```
可以看到二者的结果还是相似的。

```{r}
relaimpo::boot.relimp(fit4)
```

R in action 中所提供的函数：
```{r}
relweights <- function(fit, ...) {
  R <- cor(fit$model)
  nvar <- ncol(R)
  rxx <- R[2:nvar, 2:nvar]
  rxy <- R[2:nvar, 1]
  svd <- eigen(rxx)
  evec <- svd$vectors
  ev <- svd$values
  delta <- diag(sqrt(ev))
  lambda <- evec %*% delta %*% t(evec)
  lambdasq <- lambda^2
  beta <- solve(lambda) %*% rxy
  rsquare <- colSums(beta^2)
  rawwgt <- lambdasq %*% beta^2
  import <- (rawwgt / rsquare) * 100
  import <- as.data.frame(import)
  row.names(import) <- names(fit$model[2:nvar])
  names(import) <- "Weights"
  import <- import[order(import), 1, drop = FALSE]
  dotchart(import$Weights,
    labels = row.names(import),
    xlab = "% of R-Square", pch = 19,
    main = "Relative Importance of Predictor Variables", sub = paste("Total R-Square=", round(rsquare, digits = 3)), ...
  )
  return(import)
}
```


书中所提供的这种更为合理的方法，所得到的结果和上文的结果差异蛮大呀。。。
```{r}
relweights(fit4, col="blue")

```






#### 一些线性模型应用中常用到的指标


像在`glm.Rmd`中演示的那样，用glm或lm所得到的线性模型结果是一致的。
使用的数据为lasso中的imdb数据集。
数据输入格式数据探索的内容在此先不考虑。
```{r}
lreg <- glm(imdb_rating ~ ., data=dat, family = 'gaussian')
```


绘制某变量的 校正其它变量后，再作图观察某一个自变量与应变量之间的关系, 亦可同时绘制全部。既是矫正后的散点图。
此处绘制`kelly`变量的散点图。

```{r}
visreg::visreg(lreg, 'kelly')
```


标准化回归系数, 也可以将数据scale后在进行建模，得到的结果是一致的。
```{r}
QuantPsyc::lm.beta(mymodel)
```


利用`sjPlot`包可以同时得到标准化回归系数和绘制图型
```{r}
sjPlot::plot_model(lreg,
           type = "std",      # 计算标准化回归系数
           sort.est = TRUE,    # 进行排序
           show.values = TRUE,  # 显示数值
           show.p = TRUE,     # 显示p值
           value.offset = 0.2,  # 调整数值的位置
           vline.color = "gray",    # 修改穿越0的线条颜色
           title = "Summary of the regression model")    # 添加题目


# 给出各个变量的标准化回归系数
sjPlot::tab_model(lreg, show.std = 'std')

```

### 模型变量的重要性。

```{r}
vip::vi(lreg)
```


### *通过线性模型理解dummy变量和contrast*

在R中分类变量会被R默认设置为dummy变量。


```{r}
df2 <- readRDS('./datasets/df2.rds')


lm_fit <- lm(
  score ~ time,
  data = df2
)
```

```{r}
tapply(df2$score, df2$time, mean)
```


_手动设置分类变量的编码形式_

> https://stats.oarc.ucla.edu/r/library/r-library-contrast-coding-systems-for-categorical-variables/

`stats::contr.XX`提供了很多contrast的方法。

不同的编码方式会影响分类变量的回归系数值。不过contrast应该不会变的吧？对于一个有分类变量的回归模型而言，其的方程式为相较于reference的变化，连续变量为每改变一个单位因变量的改变。

```{r, eval=FALSE}
# Don't run!!!
contrasts(df2$time)

# 设置完成后df2数据并不会发生变化，但是在后续的模型构建中会采用不同的编码方式
contrasts(df2$time) <- contr.treatment(3) # 依据因子的level数量设置


# 设置简单编码simple code
# 简单编码的Intercept为总的平均值tapply(df2$score, df2$time, mean) |> mean()
sim.code <- contr.treatment(3) - matrix(rep(1/3, 6), ncol=2)
# contrasts(df2$time) <- sim.code


# Deviation coding
contrasts(df2$time) <- contr.sum(3)

# Orthogonal Polynomial Coding
# 正交多项式编码一般用于有序分类变量，R默认对有序分类采用此编码方式
is.ordered(df2$time)
as.ordered(df2$time)
contrasts(df2$time) <- contr.poly(3)


# Helmert Coding, Reverse Helmert Coding
# Helmert Coding比较的是当前类别下的因变量平均值和后面的类别的因变量平均值。
# 和Helmert Coding刚好完全相反，每一个类别和它前面的类别比较。
# R中提供了contr.helmert，可以进行Reverse Helmert Coding：
contr.helmert(3)


# Forward Difference Coding/Backward Difference Coding
# 这种编码方式比较的是当前类别的因变量均值和它相邻的下一个类别的因变量均值；这种编码方式和前一种刚好相反，比较的是当前类别和它相邻的前一个类别的因变量均值。
# 在R中没有对应的函数，但是上面链接有提供

```

```{r}
lm_fit <- lm(
  score ~ time,
  data = df2
)
```


```{r}
summary(lm_fit)
```
*interpret* 在一个简单的线性模型中，time这一分类变量会自动变成dummy变量，Intercept的Estimate值为612.26，和因变量按照分类变量计算的均值中的Day0（即reference）均值是一样的，且-147.26=465.0000-612.2597，-64.70=547.5584-612.2597，以此类推。既是讲对于分类变量而言模型的系数为分组变量减去reference后的值。


*对于纳入了一个协变量的线性模型而言：*

```{r}
lm_fit2 <- lm(
  score ~ time + AGE,
  data = df2
)

summary(lm_fit2)
```

*interpret* 模型分类变量的Estimate值和lm_fit1的是保持一致的，但是Intercept截距的值变成了529.597，这是在纳入AGE协变量后的取值。



### regplot绘制列线图

```{r}
regplot::regplot(lreg)
```



### 利用`finalfit`一步法得到回归模型。

`finalfit`确实包含在R中做几个常见模型的整体流程，做的多了确实应该也需要一个框架似的整合。

finalfit支持的模型包括：linear lm(), logistic glm(), hierarchical logistic lme4::glmer() and Cox proportional hazards survival::coxph() regression models. 需要测试`lme4::lmer`。

```{r}
library(finalfit)


# 指定单因素的自变量
explanatory = c("jim", "kelly", "kevin", "michael")
# 指定多因素的自变量
explanatory_multi = c("kevin", "michael")
# 指定结局变量
dependent = "imdb_rating"
# 输出结果
dat %>%
  finalfit(dependent, explanatory, explanatory_multi)

# 使用ff_plot()函数将结果数据绘制森林图
dat %>%
  ff_plot(dependent, explanatory)
```