R语言中用于计算Rsquare的包rsq

た入场券 2022-09-12 06:58 19阅读 0赞

### 文章目录 ###

*  理论介绍
 *   *  线性模型情形
     *  广义线性模型情形
 *  函数介绍
 *   *  rsq()的介绍
     *   *  实例
     *  rsq.partial()的介绍
     *   *  相关文献
     *  pcor()函数介绍
     *  vresidual()函数介绍
     *   *  实例

# 理论介绍 #

## 线性模型情形 ##

R-squared（值范围 0-1）描述的 输入变量对输出变量的解释程度。在单变量线性回归中R-squared 越大，说明拟合程度越好。

数学表达式：

R 2 = S S R / T S S = 1 − R S S / T S S R^2=SSR/TSS=1-RSS/TSS R2=SSR/TSS=1−RSS/TSS  
**其中：**

*  TSS 是执行回归分析前，响应变量固有的方差
 *  RSS 是残差平方和（就是回归模型不能解释的方差）
 *  SSR 回归模型可以解释的方差

然而，只要增加了更多的变量，无论增加的变量是否和输出变量存在关系，则R-squared 要么保持不变，要么增加。

所以需要 adjusted R-squared （范围 ( − ∞ , 1 \] (-\\infty,1\] (−∞,1\]）,它会对那些增加的且不会改善模型效果的变量增加一个惩罚项。

数学表达式：

R 2 ( a d j ) = 1 − ( 1 − R 2 ) × n − 1 n − p − 1 R^2(adj)=1-(1-R^2)\\times\\frac\{n-1\}\{n-p-1\} R2(adj)=1−(1−R2)×n−p−1n−1

(p 为变量个数，n 为样本个数)

另外一个表达是：

R 2 ( a d j ) = 1 − R S S / ( n − p − 1 ) T S S / ( n − 1 ) R^2(adj)=1-\\frac\{RSS/(n-p-1)\}\{TSS/(n-1)\} R2(adj)=1−TSS/(n−1)RSS/(n−p−1)

结论，如果单变量线性回归，则使用 R-squared评估，多变量，则使用adjusted R-squared。

在单变量线性回归中，R-squared和adjusted R-squared是一致的。

另外，如果增加更多无意义的变量，则 R-squared 和adjusted R-squared之间的差距会越来越大，Adjusted R-squared会下降。但是如果加入的特征值是显著的，则adjusted R-squared 也会上升。  
参考来源：[https://zhuanlan.zhihu.com/p/340105912][https_zhuanlan.zhihu.com_p_340105912]

R语言中的req包可以对线性模型和广义线性模型的多种形式的 R 2 R^2 R2进行计算

## 广义线性模型情形 ##

针对广义线性模型R-squared通常定义为解释的不确定性的比例（the proportionate reduction in uncertainty, measured by KL divergence, due to the inclusion of regressors.）  
一种基于KL散度的R-squared定义是：  
 R K L 2 = 1 − K ( y , μ ^ ) K ( y , μ ^ 0 ) R\_\{KL\}^2=1-\\frac\{K(y,\\hat\{\\mu\})\}\{K(y,\\hat\{\\mu\}\_0)\} RKL2=1−K(y,μ^0)K(y,μ^),  
其中，  
 K ( y , μ ^ ) = ∑ i = 1 n \[ log ⁡ f y ( y i ) − log ⁡ f μ ^ ( y i ) \] K(y,\\hat\{\\mu\})=\\sum\_\{i=1\}^n\[\\log f\_y(y\_i)-\\log f\_\{\\hat\{\\mu\}\}(y\_i)\] K(y,μ^)=∑i=1n\[logfy(yi)−logfμ^(yi)\],  
 μ ^ i = exp ⁡ ( x i β ^ ) \\hat\{\\mu\}\_i=\\exp(x\_i \\hat\{\\beta\}) μ^i=exp(xiβ^)表示模型估计的均值  
 μ ^ 0 \\hat\{\\mu\}\_0 μ^0是对 μ \\mu μ直接估计的极大似然估计，往往是 y y y的均值

**参考文献：**  
A. C. Cameron, F. A. G. Windmeijer. An R-squared measure of goodness of fit for some common nonlinear regression models. Journal of Econometrics. 1997, 77: 329-342.

# 函数介绍 #

下面的函数介绍来源于：[https://cran.r-project.org/web/packages/rsq/rsq.pdf][https_cran.r-project.org_web_packages_rsq_rsq.pdf]

## rsq()的介绍 ##

rsq(fitObj,adj=FALSE,type=c('v','kl','sse','lr','n'))

**参数介绍：**

*  fitObj：一个“lm”,“glm”, “merMod”, “lmerMod” 或者类 “lme”; 通常是 lm, glm, glm.nb, lmer , glmer , lme4中的glmer ，nlme中 lme的一个结果。
 *  adj：这是一个逻辑值，如果是TRUE，则计算adjusted R^2
 *  type: 表示R-squared的类型 (仅对广义线性模型有用):  
    ’v’ (default) – variance-function-based (Zhang, 2016), calling rsq.v;  
    ’kl’ – KL-divergence-based (Cameron and Windmeijer, 1997), calling rsq.kl;  
    ’sse’ – SSE-based (Efron, 1978), calling rsq.sse;  
    ’lr’ – likelihood-ratio-based (Maddala, 1983; Cox and Snell, 1989; Magee,  
    1990), calling rsq.lr;  
    ’n’ – corrected version of ’lr’ (Nagelkerke, 1991), calling rsq.n

**输出参数：**  
该函数除了返回R-square的值，针对（广义）线性混合模型还有下列输出：

*  R\_M^2: 模型总共可以解释的变异性的比例，包含固定效应以及随机效应.
 *  R\_F^2：可以被固定因子解释的变异性所占的比例
 *  R\_R^2：可以被随机效应因子解释的变异性所占的比例

### 实例 ###

library(rsq)
    data(hcrabs)
    attach(hcrabs)
    y <- ifelse(num.satellites>0,1,0)
    bnfit <- glm(y~color+spine+width+weight,family=binomial)
    rsq(bnfit)
    # [1] 0.2171238
    rsq(bnfit,adj=TRUE)
    # [1] 0.1839109
    quasibn <- glm(y~color+spine+width+weight,family=quasibinomial)
    rsq(quasibn)
    # [1] 0.2171238
    rsq(quasibn,adj=TRUE)
    # [1] 0.1839109
    psfit <- glm(num.satellites~color+spine+width+weight,family=poisson)
    rsq(psfit)
    # [1] 0.1172267
    rsq(psfit,adj=TRUE)
    # [1] 0.07977572
    quasips <- glm(num.satellites~color+spine+width+weight,family=quasipoisson)
    rsq(quasips)
    # [1] 0.1172267
    rsq(quasips,adj=TRUE)
    # [1] 0.07977572
    
    
    
    # Linear mixed models
    require(lme4)
    lmm1 <- lmer(Reaction~Days+(Days|Subject),data=sleepstudy)
    rsq(lmm1)
    # $model
    # [1] 0.8003832
    # 
    # $fixed
    # [1] 0.2864714
    # 
    # $random
    # [1] 0.5139119
    rsq.lmm(lmm1)
    # $model
    # [1] 0.8003832
    # 
    # $fixed
    # [1] 0.2864714
    # 
    # $random
    # [1] 0.5139119

## rsq.partial()的介绍 ##

该函数用于计算aka partial  R 2 R^2 R2, 可用于线性模型和广义线性模型

rsq.partial(objF,objR=NULL,adj=FALSE,type=c('v','kl','sse','lr','n'))

**参数介绍：**

*  objF：是一个 “lm” 或 "glm"对象, 也是lm, glm, 或者glm.nb 的用于拟合**全模型**的结果
 *  objR：是一个 “lm” 或 "glm"对象, 也是lm, glm, 或者glm.nb 的用于拟合**缩减模型**的结果
 *  adj: logical; if TRUE, calculate the adjusted partial R^2.  
    type: R-squared的类型:  
    ’v’ (default) – variance-function-based (Zhang, 2016), calling rsq.v;  
    ’kl’ – KL-divergence-based (Cameron and Windmeijer, 1997), calling rsq.kl;  
    ’sse’ – SSE-based (Efron, 1978), calling rsq.sse;  
    ’lr’ – likelihood-ratio-based (Maddala, 1983; Cox and Snell, 1989; Magee,1990), calling rsq.lr;  
    ’n’ – corrected version of ’lr’ (Nagelkerke, 1991), calling rsq.n.

**输出参数：**  
返回值包含adjustment 和partial.rsq. 当objR 为 NULL时, variable.full和 variable.reduced将被返回; 否则返回的是variable。

*  adjustment：逻辑值; 如果为TRUE， 则计算adjusted partial R^2.
 *  variable.full ：全模型中全部的协变量
 *  variable.reduced：缩减模型中的全部协变量.
 *  variable： 全模型中全部的协变量.
 *  partial.rsq： partial  R 2 R^2 R2 或者 adjusted partial  R 2 R^2 R2.

### 相关文献 ###

1.  Cameron, A. C. and Windmeijer, A. G. (1997) An R-squared measure of goodness of fit for some common nonlinear regression models. Journal of Econometrics, 77: 329-342.
2.  Cox, D. R. and Snell, E. J. (1989) The Analysis of Binary Data, 2nd ed. London: Chapman and Hall.
3.  Efron, B.(1978) Regression and ANOVA with zero-one data: measures of residual variation. Journal of the American Statistical Association, 73: 113-121.
4.  Maddala, G. S. (1983) Limited-Dependent and Qualitative Variables in Econometrics. Cambridge University.
5.  Magee, L. (1990) R^2 measures based on Wald and likelihood ratio joint significance tests. The American Statistician, 44: 250-253.
6.  Nagelkerke, N. J. D. (1991) A note on a general definition of the coefficient of determination. Biometrika, 78: 691-692.
7.  Zhang, D. (2017). A coefficient of determination for generalized linear models. The American Statistician, 71(4): 310-316.

## pcor()函数介绍 ##

该函数用于计算线性模型和广义线性模型的偏相关系数

pcor(objF,objR=NULL,adj=FALSE,type=c('v','kl','sse','lr','n'))

**参数介绍：**

## vresidual()函数介绍 ##

用于计算variance-function-based的残差, 这是用来计算variance-function-based R-squared.

vresidual(y,yfit,family=binomial(),variance=NULL)

**参数介绍：**

*  y：观测值的向量
 *  yfit：拟合值的向量
 *  family：分布族
 *  variance：方差函数 (specified by family by default).

**注意：**  
残差的计算依赖于方差函数, 对于quasi 模型需要很好的定义. 当方差函数是常数或者线性函数时，这会导出经典的残差。注意只有方差函数需要被设定，通过 “family”" 或者 “variance”。

### 实例 ###

data(hcrabs)
    attach(hcrabs)
    y <- ifelse(num.satellites>0,1,0)
    bnfit <- glm(y~color+spine+width+weight,family="binomial")
    vresidual(y,bnfit$fitted.values,family="binomial")
    # Effectiveness of Bycycle Safety Helmets in Thompson et al. (1989)
    y <- matrix(c(17,218,233,758),2,2)
    x <- factor(c("yes","no"))
    tbn <- glm(y~x,family="binomial")
    yfit <- cbind(tbn$fitted.values, 1-tbn$fitted.values)
    vr0 <- vresidual(matrix(0,2,1),yfit[,1],family="binomial")
    vr1 <- vresidual(matrix(1,2,1),yfit[,2],family="binomial")
    y[,1]*vr0+y[,2]*vr1

[https_zhuanlan.zhihu.com_p_340105912]: https://zhuanlan.zhihu.com/p/340105912
[https_cran.r-project.org_web_packages_rsq_rsq.pdf]: https://cran.r-project.org/web/packages/rsq/rsq.pdf