• lm()
  即linear model线性模型函数，用来建立OLS回归模型

m1<-lm(iris$Sepal.Length~iris$Sepal.Width+iris$Petal.Length+iris$Petal.Width)
summary(m1)#查看m1详细情况

• OLS线性回归

####OLS回归
# OLS（最小二乘法）主要用于线性回归的参数估计，
# 它的思路很简单，就是求一些使得实际值和模型估值之差的平方和达到最小的值，
# 将其作为参数估计值。就是说，通过最小化误差的平方和寻找数据的最佳函数匹配。
# 利用最小二乘法可以简便地求得未知的数据，
# 并使得这些求得的数据与实际数据之间误差的平方和为最小。
# 最小二乘法可用于曲线拟合，
# 其他一些优化问题也可通过最小化能量或最大化熵用最小二乘法来表达。
# OLS法通过一系列的预测变量来预测响应变量（也可以说是在预测变量上回归响应变量）。
# 线性回归是指对参数β为线性的一种回归（即参数只以一次方的形式出现）模型：
# Yt=α+βxt+μt （t=1……n）表示观测数
# 
# Yt  被称作因变量
# xt  被称作自变量
# α、β 为需要最小二乘法去确定的参数，或称回归系数
# μt  为随机误差项
# OLS线性回归的基本原则：最优拟合曲线应该使各点到直线的距离的平方和（即残差平方和，简称RSS）最小：
# 
# 
# 
# OLS线性回归的目标是通过减少响应变量的真实值与预测值的差值来获得模型参数（截距项和斜率），
# 就是使RSS最小。
# 
# 为了能够恰当地解释OLS模型的系数，数据必须满足以下统计假设：
# 
# 正态性：对于固定的自变量值，因变量值成正太分布
# 独立性：个体之间相互独立
# 线性相关：因变量和自变量之间为线性相关
# 同方差性：因变量的方差不随自变量的水平不同而变化，即因变量的方差是不变的
#
####用lm()拟合回归模型
# 在R中，拟合回归模型最基本的函数是lm()，格式为：
# lm(formula, data)
#
# formula中的符号注释：
# 
# ~ 分割符号，左边为因变量，右边为自变量，例如， z~x+y，表示通过x和y来预测z
# + 分割预测变量
# : 表示预测变量的交互项，例如，z~x+y+x:y
# * 表示所有可能的交互项，例如，z~x*y 展开为 z~x+y+x:y
# ^ 表示交互项的次数，例如，z ~ (x+y)^2，展开为z~x+y+x:y
# . 表示包含除因变量之外的所有变量，例如，如果只有三个变量x，y和z，那么代码 z~. 
#展开为z~x+y+x:y
# -1 删除截距项，强制回归的直线通过原点
# I() 从算术的角度来解释括号中的表达式，例如，z~y+I(x^2) 表示的拟合公式是 z=a+by+cx2
# function 可以在表达式中应用数学函数，例如，log(z) ~x+y
# 对于拟合后的模型（lm函数返回的对象），可以应用下面的函数，得到模型的更多额外的信息。
# 
# summary() 展示拟合模型的详细结果
# coefficients() 列出捏模型的参数（截距项intercept和斜率）
# confint() 提供模型参数的置信区间
# residuals() 列出拟合模型的残差值
# fitted() 列出拟合模型的预测值
# anova() 生成一个拟合模型的方差分析表
# predict() 用拟合模型对新的数据预测响应变量
# 例如，使用基础安装包中的数据集women（包含了年龄在30-39岁之间女性的身高和体重信息），通过身高来预测体重。
# 
#
#
# 1，线性拟合
# 
# 使用lm()函数来拟合模型，获得模型对象fit
fit <- lm(weight~height,data=women)
# 使用summary()函数来查看模型的信息：
summary(fit)

# Call:
#   lm(formula = weight ~ height, data = women)
# 
# Residuals:
#   Min      1Q  Median      3Q     Max 
# -1.7333 -1.1333 -0.3833  0.7417  3.1167 
# 
# Coefficients:
#   Estimate Std. Error t value Pr(>|t|)    
# (Intercept) -87.51667    5.93694  -14.74 1.71e-09 ***
#   height        3.45000    0.09114   37.85 1.09e-14 ***
# ---
#   Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
# 
# Residual standard error: 1.525 on 13 degrees of freedom
# Multiple R-squared:  0.991,    Adjusted R-squared:  0.9903 
# F-statistic:  1433 on 1 and 13 DF,  p-value: 1.091e-14
# 残差标准差（Residual standard error）：表示模型用身高来预测体重的平均误差
# 
# R的平方项（Multiple R-squared）：表明模型可以解释体重99.1%的方差，
#是实际和预测值之间相关系数的平方。
# 
# F统计量（F-statistic）：检验所有的预测变量预测响应变量是否都在某个概率水平之上，
# 
# 从Coefficients 组中，可以看到 Intercept（截距项）是 - 87.51667，height的系数是3.45，
#截距项和系数的标准差、t值和Pr(>|t|)，其中，Pr(>|t|) 表示双边检验的p值
# 
# 注，p值的表示方法通常有p-value，或Pr，p值是概率，表示某一事件发生的可能性大小。
#如果P值很小，说明原假设情况的发生的概率很小，而如果出现了，根据小概率原理，
#我们就有理由拒绝原假设，P值越小，我们拒绝原假设的理由越充分。
#总之，P值越小，表明结果越显著。
#
# 因此可以得到预测等式：
# Weight= - 87.51667 + 3.45 * Height
# 绘制拟合的直线：
plot(women$height, women$weight,xlab='height',ylab='weight')
abline(fit)
#
#
#
#2、多项式拟合
# 使用lm()，**在formula参数中使用I()函数来进行多项式拟合**
fit2 <- lm(weight~I(height^2),data=women)
plot(women$height, women$weight,xlab='height',ylab='weight')
lines(women$height,fitted(fit2))
summary(fit2)
# Call:
#   lm(formula = weight ~ I(height^2), data = women)
# 
# Residuals:
#   Min      1Q  Median      3Q     Max 
# -1.2562 -0.7636 -0.1837  0.4622  2.2654 
# 
# Coefficients:
#   Estimate Std. Error t value Pr(>|t|)    
# (Intercept) 2.390e+01  2.109e+00   11.34 4.12e-08 ***
#   I(height^2) 2.659e-02  4.926e-04   53.98  < 2e-16 ***
#   ---
#   Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
# 
# Residual standard error: 1.072 on 13 degrees of freedom
# Multiple R-squared:  0.9956,    Adjusted R-squared:  0.9952 
# F-statistic:  2913 on 1 and 13 DF,  p-value: < 2.2e-16
#
#非线性模型可用nls()函数进行拟合。
#
#
#
#####回归诊断（标准方法）
# 回归诊断用于评价回归模型的拟合程度，模型返回的参数多大程度上匹配原始数据。
# 
# 根据OLS回归的统计假设来评价模型的拟合情况，对于lm()拟合的模型对象，
# 使用plot()函数生成评价模型拟合情况的四幅图形，
weight<-c(80,84,85,87,88,89,90,78)
height<-c(180,175,174,156,165,176,177,178)
fit <- lm(weight~height,data=women)
par(mfrow=c(2,2))
plot(fit)
# 正态性：当预测变量值固定时，因变量成正态分布，则残差值也应该是一个均值为0的正态分布。
# 正态QQ图（Normal Q-Q）是在正态分布对应的值下，标准化残差的概率图。若满足正态假设，
# 那么图上的点应该落在呈45度角的直线上；若不是如此，那么久违反了正态假设。
# 
# 独立性：无法判断因变量和自变量的值是否相互独立，只能从收集的数据中来验证。
# # 
# 线性相关性：若因变量与自变量线性相关，那么残差值与预测（拟合）值除了系统误差之外，
# 就没有任何关联。在残差和拟合图（Residuals VS Fitted）中，可以看到一条曲线，
# 这暗示着可能需要对拟合模型加上一个二次项。
# 
# 同方差性：若满足不变方差假设，那么在位置尺度图（Scale-Location）中，水平线周围的点应该随机分布。
# 
# 残差和杠杆图（Residuals VS Leverage）提供了特殊的单个观测点（离群点、高杠杆点、强影响点）的信息。
# 
# 离群点：表明拟合回归模型对其预测效果不佳（产生巨大的残差）
# 高杠杆点：指自变量因子空间中的离群点（异常值），由许多异常的自变量值组合起来的，与因变量没有关系。
# 强影响点：表明它对模型参数的估计产生的影响过大，非常不成比例。
# 强影响点可以通过Cook距离（Cook distance）统计量来前别。
# 残差和杠杆图的可读性差，不够实用。
# 
# 什么是高杠杆点？离群点是指对于给定的预测值  来说，响应值  异常的点。
# 相反，高杠杆(high leverage) 表示观测点  是异常的。例如，图1(a)左下图中的观测点41具有高杠杆值，
# 因为它的预测变量值比其他观测点都要大。实线是对数据的最小二乘拟合，而虚线是删除观测点41后的拟合。
# 事实上，高杠杆的观测往往对回归直线的估计有很大的影响。如果一些观测对最小二乘线有重大影响，
# 那么它们值得特别关注，这些点出现任何问题都可能使整个拟合失效。因此找出高杠杆观测是十分重要的 。
# 
# 
# 
# ####回归诊断（car包）
# car包提供了大量的函数，大大增强了拟合和评价回归模型的能力
# 
# 1，正态性
# 
# qqplot()函数提供了精确的正态假设检验方法，
install.packages("carData")
install.packages("car")
library(carData)
library(car)
par(mfrow=c(1,2))
fit <- lm(weight~height,data=women)
qqPlot(fit,labels=row.names(women),id.method='identity',simulate=TRUE,main='qq-fit')
# [1]  1 15
fit2 <- lm(weight~height+I(height^2),data=women)
qqPlot(fit2,labels=row.names(women),id.method='identity',simulate=TRUE,main='qq-fit2')
# [1] 13 15
# 
# 2，误差的独立性
# car包提供了durbinWatsonTest()函数，用于做Durbin-Watson检验，检测误差的序列相关性。
durbinWatsonTest(fit)
# p值 （p=0）不显著，误差项之间独立
# 
# 
# 3，线性相关性
# 
# 通过成分残差图（component + residual plots）检查因变量和自变量之间是否呈线性关系。
crPlots(fit)
# 若图形存在非线性，则说明可能对预测变量的函数形式建模不够充分，那么需要添加一些曲线成分，比如多项式，对一个或多个变量进行变换（log(x)代替x），或用其他回归变体形式而不是线性回归。
# 
# 
# 4，同方差性
# 判断方差是否恒定，nvcTest()函数生成一个记分检验，原假设为误方差不变，
# 备择假设为误差方差随着拟合值水平的变化而变化。若检验显著，说明存在误方差不很定。
# 
# spreadLevelPlot()函数创建了一个添加了最佳拟合曲线的散点图，
# 展示了标准化残差绝对值与拟合值得关系。
ncvTest(fit)
spreadLevelPlot(fit)
# 记分检验不显著：p=0.36954，说明满足方差不变假设，也可以通过分布水平看到这一点，
# 点在水平的最佳拟合曲线周围呈随机分布。
# 
# 
# ####回归诊断（gvlma包）
# gvlma()函数，用于对线性模型假设进行综合验证，同事还能验证偏斜度，峰度和异方差的评价。
# 从输出项中可以看出，假设检验的显著性水平是5%。如果p<0.05，说明违反了假设条件。
# 从Global Stat 的Decision 栏中，可以看到数据满足OLS回归模型所有的统计假设（p=0.597）。
install.packages("gvlma")
library(gvlma)
gvmodel <- gvlma(fit)
summary(gvmodel)                  

• LASSO回归

LASSO是由1996年Robert Tibshirani首次提出，全称Least absolute shrinkage and selection operator。该方法是一种压缩估计。它通过构造一个惩罚函数得到一个较为精炼的模型，使得它压缩一些回归系数，即强制系数绝对值之和小于某个固定值；同时设定一些回归系数为零。因此保留了子集收缩的优点，是一种处理具有复共线性数据的有偏估计。

image.png
image.png

#Load the Data
#define response variable
y <- mtcars$hp

#define matrix of predictor variables
x <- data.matrix(mtcars[, c('mpg', 'wt', 'drat', 'qsec')])

#Fit the Lasso Regression Model

library(glmnet)

#perform k-fold cross-validation to find optimal lambda value
cv_model <- cv.glmnet(x, y, alpha = 1)

#find optimal lambda value that minimizes test MSE
best_lambda <- cv_model$lambda.min
best_lambda



#produce plot of test MSE by lambda value
plot(cv_model)


# Analyze Final Model
#find coefficients of best model
best_model <- glmnet(x, y, alpha = 1, lambda = best_lambda)
coef(best_model)

# 5 x 1 sparse Matrix of class "dgCMatrix"
# s0
# (Intercept) 484.20742
# mpg          -2.95796
# wt           21.37988
# drat          .      
# qsec        -19.43425



# predict the value for hp of this new observation:
new = matrix(c(24, 2.5, 3.5, 18.5), nrow=1, ncol=4) 

#use lasso regression model to predict response value
predict(best_model, s = best_lambda, newx = new)

# [1,] 109.0842



# calculate the R-squared of the model on the training data:
#use fitted best model to make predictions
y_predicted <- predict(best_model, s = best_lambda, newx = x)

#find SST and SSE
sst <- sum((y - mean(y))^2)
sse <- sum((y_predicted - y)^2)

#find R-Squared
rsq <- 1 - sse/sst
rsq

# [1] 0.8047064

• Ridge回归
  注：在存在共线性问题和病态数据偏多的研究中有较大的实用价值。
  岭回归(英文名：ridge regression, Tikhonov regularization)是一种专用于共线性数据分析的有偏估计回归方法，实质上是一种改良的最小二乘估计法，通过放弃最小二乘法的无偏性，以损失部分信息、降低精度为代价获得回归系数更为符合实际、更可靠的回归方法，对病态数据的拟合要强于最小二乘法。
  image.png
  image.png

特点：

通常岭回归方程的R平方值会稍低于普通回归分析，但回归系数的显著性往往明显高于普通回归，在存在共线性问题和病态数据偏多的研究中有较大的实用价值。

# 首先读取数据，可视化变量间的关系：
library(readxl);library(caret);library(glmnet);library(corrplot)
library(Metrics);library(ggplot2)
## 读取数据
#diabete <- read.csv("GSE84465_GBM_All_data.csv",sep = "\t")
str(iris)
diabete <- iris[,1:3]
## 可视化相关系数
diabete_cor <- cor(diabete)
corrplot.mixed(diabete_cor,tl.col="black",tl.pos = "d",number.cex = 0.8)
#在建立Ridge回归模型之前，先将数据切分为训练集和测试集，并对数据集进行标准化，
#使用标准化后的数据进行建模。
## 切分为训练集核测试集，70%训练，30%测试
set.seed(123)
d_index <- createDataPartition(diabete$Petal.Length,p = 0.7)
train_d <- diabete[d_index$Resample1,]
test_d <- diabete[-d_index$Resample1,]
## 数据标准化
scal <- preProcess(train_d,method = c("center","scale"))
train_ds <- predict(scal,train_d)
test_ds <- predict(scal,test_d)
## 查看标准化使用的均值和标准差
scal$mean
scal$std
#
# 需要指定一个合适的参数lambda，该参数为施加在回归系数上的惩罚系数，
#不同的参数lambda可以得到不同的Ridge回归模型。
# glmnet包中的cv.glmnet()可以通过交叉验证的方式，来分析在不同的参数lambda下回归模型的效果。
# 下面使用cv.glmnet()函数和训练数集分析模型参数的影响。
#
lambdas <- seq(0,5, length.out = 200)
X <- as.matrix(train_ds[,1:3])
Y <- train_ds[,2]
set.seed(1245)
ridge_model <- cv.glmnet(X,Y,alpha = 0,lambda = lambdas,nfolds =3)
plot(ridge_model)
plot(ridge_model$glmnet.fit, "lambda", label = T)
# 使用cv.glmnet()函数来进行参数lambda的影响分析。通过指定参数alpha = 0来建立Ridge回归，
#如果参数alpha = 1，则建立的是Lasso回归模型，nfolds =3表示使用3折交叉验证。
#使用plot(ridge_model)可视化lambda对模型均方误差的影响，而plot(ridge_model$glmnet.fit,…)
# 用来可视化不同的lambda下，各个自变量回归系数的变化，
# 在cv.glmnet()函数的输出结果中，包含一个lambda.min取值，且ridge_model$lambda.min = 0.15，
# 该值为在模型均方误差最小的情况下参数lambda的取值，可以使用该值建立Ridge回归模型。
# 下面使用该参数对全部的训练数据集建立Ridge模型。
ridge_min <- ridge_model$lambda.min
## 使用ridge_min 拟合ridge模型
ridge_best <- glmnet(X,Y,alpha = 0,lambda = ridge_min)
coef(ridge_best)
# 为了验证Ridge回归模型的预测效果，可在测试集上进行预测，并计算出平均绝对值误差的大小。
# 
test_pre <- predict(ridge_best,as.matrix(test_ds[,1:3]))
sprintf("标准化后平均绝对误差为: %f",mae(test_ds$Y,test_pre))
#
## 将预测值逆标准化和原始数据进行比较
test_pre_o <- as.vector(test_pre[,1] * scal$std[2] + scal$mean[3])
sprintf("标准化前平均绝对误差为: %f",mae(test_d$Y,test_pre_o))

• Logistic回归
  logistic 回归又叫对数几率回归，适合数值型的二值型输出的拟合，它其实是一个分类模型，比如根据患者的医疗数据判断它是否能被治愈。最大似然法来解决方程估计和检验问题。something that can take two values such as true/false, yes/no, and so on.
  原理：Logistic回归原理及公式推导https://blog.csdn.net/AriesSurfer/article/details/41310525

# How to Perform Logistic Regression in R (Step-by-Step)
# Step 1: Load the Data
#load dataset
install.packages("ISLR")
data <- ISLR::Default

#view summary of dataset
summary(data)

# default    student       balance           income     
# No :9667   No :7056   Min.   :   0.0   Min.   :  772  
# Yes: 333   Yes:2944   1st Qu.: 481.7   1st Qu.:21340  
# Median : 823.6   Median :34553  
# Mean   : 835.4   Mean   :33517  
# 3rd Qu.:1166.3   3rd Qu.:43808  
# Max.   :2654.3   Max.   :73554  

#find total observations in dataset
nrow(data)

# [1] 10000




# Step 2: Create Training and Test Samples
#make this example reproducible
set.seed(1)

#Use 70% of dataset as training set and remaining 30% as testing set
sample <- sample(c(TRUE, FALSE), nrow(data), replace=TRUE, prob=c(0.7,0.3))
train <- data[sample, ]
test <- data[!sample, ]  
# Step 3: Fit the Logistic Regression Model
#fit logistic regression model
model <- glm(default~student+balance+income, family="binomial", data=train)

#disable scientific notation for model summary
options(scipen=999)

#view model summary
summary(model)

# Call:
#   glm(formula = default ~ student + balance + income, family = "binomial", 
#       data = train)

# Deviance Residuals: 
#   Min       1Q   Median       3Q      Max  
# -2.5586  -0.1353  -0.0519  -0.0177   3.7973  

# Coefficients:
#   Estimate    Std. Error z value            Pr(>|z|)    
# (Intercept) -11.478101194   0.623409555 -18.412 <0.0000000000000002 ***
#   studentYes   -0.493292438   0.285735949  -1.726              0.0843 .  
# balance       0.005988059   0.000293765  20.384 <0.0000000000000002 ***
#   income        0.000007857   0.000009965   0.788              0.4304  
# ---
#   Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# 
# (Dispersion parameter for binomial family taken to be 1)
# 
# Null deviance: 2021.1  on 6963  degrees of freedom
# Residual deviance: 1065.4  on 6960  degrees of freedom
# AIC: 1073.4
# 
# Number of Fisher Scoring iterations: 8
# Assessing Model Fit:
caret::varImp(model)

# Overall
# studentYes  1.726393
# balance    20.383812
# income      0.788449
# VIF Values:
  #calculate VIF values for each predictor variable in our model
car::vif(model)
# student  balance   income 
# 2.754926 1.073785 2.694039




# Step 4: Use the Model to Make Predictions
#define two individuals
new <- data.frame(balance = 1400, income = 2000, student = c("Yes", "No"))

#predict probability of defaulting
predict(model, new, type="response")

# 1          2 
# 0.02732106 0.04397747
#calculate probability of default for each individual in test dataset
predicted <- predict(model, test, type="response")

# Step 5: Model Diagnostics
install.packages("InformationValue")
library(InformationValue)

#convert defaults from "Yes" and "No" to 1's and 0's
test$default <- ifelse(test$default=="Yes", 1, 0)

#find optimal cutoff probability to use to maximize accuracy
optimal <- optimalCutoff(test$default, predicted)[1]
optimal

# [1] 0.5451712
confusionMatrix(test$default, predicted)

# 0  1
# 0 2912 64
# 1   21 39
#calculate sensitivity
sensitivity(test$default, predicted)

# [1] 0.3786408

#calculate specificity
specificity(test$default, predicted)

# [1] 0.9928401

#calculate total misclassification error rate
misClassError(test$default, predicted, threshold=optimal)

# [1] 0.027

#plot the ROC curve
plotROC(test$default, predicted)

• Poission回归
  Poisson regression is a special type of regression in which the response variable consists of “count data.

Assumptions of Poisson Regression

（Before we can conduct a Poisson regression, we need to make sure the following assumptions are met so that our results from the Poisson regression are valid）:

Assumption 1: The response variable consists of count data.

Assumption 2: Observations are independent.

Assumption 3: The distribution of counts follows a Poisson distribution.

Assumption 4: The mean and variance of the model are equal.

实现代码暂略

• 生存分析和cox回归
  从采取的分析方法上看，生存分析有非参数法（如Wilcoxon法、Log-rank法）、参数法（如Weibull回归、lognormal回归等）和半参数分析（Cox回归）。 Cox回归要求满足比例风险假定(proportional-hazards assumption)的前提条件。Cox回归分为单因素cox回归和多因素cox回归。

实现代码暂略

• 多维标度MDS分析
  MDS(多维尺度变换)

多维尺度变换算法解决的问题是:当n个对象之间的相似性给定，确定这些对象在低维空间中的表示，并使其尽可能与原先的相似性大致匹配。

其中MDS又分为classical MDS和no-classical MDS。

Classical MDS(经典多维尺度变换):经典多维尺度变换的距离标准采用欧式距离。

No-classical MDS(非经典多维度尺度变换):非经典多维度尺度变换的距离标准采用非欧式距离。

实现代码暂略

• 负二项回归
  做lrtest 如果lrtest成立 就用possion，如果不成立，就用负二项。如果研究X对于Y的影响，Y是计数资料，一般可以使用Poisson回归进行研究。但是Poisson回归要求数据满足等离散现象（平均值与方差相等），如果说数据具有一定的聚焦性，此时很可能就会产生过离散现象，即数据平均值与方差明显不相等。此时使用负二项回归更为科学。

实现代码暂略

参考：
https://www.cnblogs.com/ljhdo/p/4807068.html
https://blog.csdn.net/xiaozhu_1024/article/details/80585151
https://zhuanlan.zhihu.com/p/166179138
https://zhuanlan.zhihu.com/p/377988205
https://baike.baidu.com/item/%E5%B2%AD%E5%9B%9E%E5%BD%92/554917
https://www.statology.org/lasso-regression-in-r/
https://zhuanlan.zhihu.com/p/95132284
https://baike.baidu.com/item/logistic%E5%9B%9E%E5%BD%92/2981575
https://www.sciencedirect.com/topics/computer-science/logistic-regression
https://www.zhihu.com/question/264754736/answer/1538471612
https://zhuanlan.zhihu.com/p/51441355