为了考察一种新的经济学教学方法对学生成绩的影响，进行了调查，共得到了32个样本数据。数据见以下概况介绍。GRADE取1时表示新近学习成绩提高，取0时表示其他；GPA是平均积分点；TUCE是以往经济学成绩；PSI取1时表示受到新的经济学教学方法的指导，取0时表示其他。假如想要了解GPA、TUCE和PSI因素对学生成绩是否有影响，以及根据学生的GPA、TUCE和PSI预测学生成绩是否会提高，该如何建模分析？
mdata <- read.csv('http://data.galaxystatistics.com/blog_data/regression/grade_data.csv',header=T)
head(mdata)## obs GRADE GPA TUCE PSI
## 1 1 0 2.66 20 0
## 2 2 0 2.89 22 0
## 3 3 0 3.28 24 0
## 4 4 0 2.92 12 0
## 5 5 1 4.00 21 0
## 6 6 0 2.86 17 0modeldata <- subset(mdata, select=c(GRADE,GPA,TUCE,PSI))
#描述性统计计算
modeldata1 <- subset(modeldata, PSI==0)
modeldata2 <- subset(modeldata, PSI==1)
summary(modeldata1)## GRADE GPA TUCE PSI
## Min. :0.0000 Min. :2.630 Min. :12.00 Min. :0
## 1st Qu.:0.0000 1st Qu.:2.777 1st Qu.:19.25 1st Qu.:0
## Median :0.0000 Median :2.905 Median :21.50 Median :0
## Mean :0.1667 Mean :3.101 Mean :21.56 Mean :0
## 3rd Qu.:0.0000 3rd Qu.:3.310 3rd Qu.:24.75 3rd Qu.:0
## Max. :1.0000 Max. :4.000 Max. :29.00 Max. :0summary(modeldata2)## GRADE GPA TUCE PSI
## Min. :0.0000 Min. :2.060 Min. :14.00 Min. :1
## 1st Qu.:0.0000 1st Qu.:2.845 1st Qu.:21.00 1st Qu.:1
## Median :1.0000 Median :3.140 Median :23.00 Median :1
## Mean :0.5714 Mean :3.138 Mean :22.43 Mean :1
## 3rd Qu.:1.0000 3rd Qu.:3.533 3rd Qu.:24.75 3rd Qu.:1
## Max. :1.0000 Max. :4.000 Max. :28.00 Max. :1#数据集字段名称转换
renames <- function(x) {
ocol <- ncol(x)
if (ocol==2) {
xnam <- paste0("X")
names(x) <- c('Y', xnam)
}else if(ocol>2) {
xnam <- paste0("X", 1:(ocol-1))
names(x) <- c('Y', xnam)
}
return(x)
}
modeldata <- renames(modeldata)
head(modeldata)## Y X1 X2 X3
## 1 0 2.66 20 0
## 2 0 2.89 22 0
## 3 0 3.28 24 0
## 4 0 2.92 12 0
## 5 1 4.00 21 0
## 6 0 2.86 17 0class(modeldata)## [1] "data.frame"dim(modeldata)## [1] 32 4regFormula <- as.formula(paste("Y ~ ", paste(names(modeldata)[-1], collapse= "+")))
regFormula## Y ~ X1 + X2 + X3由于因变量是二元选择的结果，因此按传统线性回归模型所计算的判定系数R^2不再有实际的意义。可以定义为\[CountR^2=\frac{正确预测的个数}{总观测值个数}\]当\(Y\)的实际预测值大于\(0.5\)时，视其预测值为\(1\)；当小于\(0.5\)时，视其预测值为\(0\)。然后比较预测值和实际值是否存在差异，如果不存在差异，则认为是正确地预测。最后将正确的个数与总预测个数比较，得到一个新的拟合优度指标。
###线性回归
#前6行数
# head(modeldata)
#行列数
# dim(modeldata)
fit <- lm(regFormula, data=modeldata)
#展示拟合模型的详细结果
summary(fit)##
## Call:
## lm(formula = regFormula, data = modeldata)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.78153 -0.27731 0.00531 0.21089 0.81145
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1.49802 0.52389 -2.859 0.00793 **
## X1 0.46385 0.16196 2.864 0.00784 **
## X2 0.01050 0.01948 0.539 0.59436
## X3 0.37855 0.13917 2.720 0.01109 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.3881 on 28 degrees of freedom
## Multiple R-squared: 0.4159, Adjusted R-squared: 0.3533
## F-statistic: 6.646 on 3 and 28 DF, p-value: 0.001571#列出拟合模型的模型参数
coefficients(fit) ## (Intercept) X1 X2 X3
## -1.49801712 0.46385168 0.01049512 0.37855479head(predict(fit, modeldata[-1]),10)## 1 2 3 4 5 6
## -0.05426921 0.07340692 0.27529932 -0.01762875 0.57778716 0.00701576
## 7 8 9 10
## -0.03936941 0.05363477 0.16983152 0.62464001从分析的结果来看，在\(5%\)的显著性水平上，PSI对GRADE的影响是显著的。也就是说，在GPA和TUCE都一样的情况下，接受过新的教学方法的学生与没有接受过新的教学方法的学生相比，学习成绩提高的概率要多\(0.3786\)。此外，GPA对成绩提高的边际影响是\(0.46\)，也就是说，在其他条件相同的情况下，GPA每增加1，学习成绩提高的概率为\(46%\)。
Probit模型是一种广义的线性模型。服从正态分布。
最简单的Probit模型就是指被解释变量\(Y\)是一个\((0,1)\)变量，事件发生地概率是依赖于解释变量，即\(P（Y=1）= f(X)\)，也就是说,\(Y=1\)的概率是一个关于X的函数，其中\(f(.)\)服从标准正态分布。
# link设为probit
Probit.fit <- glm(regFormula, family = binomial(link="probit") , data=modeldata)
#展示拟合模型的详细结果
summary(Probit.fit)##
## Call:
## glm(formula = regFormula, family = binomial(link = "probit"),
## data = modeldata)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.9392 -0.6508 -0.2229 0.5934 2.0451
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -7.45231 2.57152 -2.898 0.00376 **
## X1 1.62581 0.68973 2.357 0.01841 *
## X2 0.05173 0.08119 0.637 0.52406
## X3 1.42633 0.58695 2.430 0.01510 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 41.183 on 31 degrees of freedom
## Residual deviance: 25.638 on 28 degrees of freedom
## AIC: 33.638
##
## Number of Fisher Scoring iterations: 6#列出拟合模型的模型参数
coefficients(Probit.fit) ## (Intercept) X1 X2 X3
## -7.45231349 1.62581226 0.05172846 1.42633069head(predict(Probit.fit, modeldata[-1], type='response'),10)## 1 2 3 4 5 6
## 0.01817085 0.05308069 0.18992678 0.01857098 0.55457680 0.02723334
## 7 8 9 10
## 0.01850346 0.04457162 0.10880826 0.66312101从分析结果看，Probit模型的设定形式是：Prob( \(GRADE_i\) = 1 | \(x_i\) ) = \(\phi ( X_l^\prime \beta)\)，其中\(\phi(∙)\)是标准正态分布的累积分布函数。将系数的估计结果代入，得到估计的模型为：
\(\widetilde{Prob}\) ( \(GRADE_i\) = 1 | \(x_i\) ) = \(\phi (-7.452320+1.625810*GPA+0.051729*TUCE+1.426331*PSI)\)
# link设为logit
Logit.fit <- glm(regFormula, family = binomial(link="logit") , data=modeldata)
#展示拟合模型的详细结果
summary(Logit.fit)##
## Call:
## glm(formula = regFormula, family = binomial(link = "logit"),
## data = modeldata)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.9551 -0.6453 -0.2570 0.5888 2.0966
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -13.02135 4.93127 -2.641 0.00828 **
## X1 2.82611 1.26293 2.238 0.02524 *
## X2 0.09516 0.14155 0.672 0.50143
## X3 2.37869 1.06456 2.234 0.02545 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 41.183 on 31 degrees of freedom
## Residual deviance: 25.779 on 28 degrees of freedom
## AIC: 33.779
##
## Number of Fisher Scoring iterations: 5#列出拟合模型的模型参数
coefficients(Logit.fit) ## (Intercept) X1 X2 X3
## -13.02134686 2.82611259 0.09515766 2.37868765head(predict(Logit.fit, modeldata[-1], type='response'),10)## 1 2 3 4 5 6
## 0.02657799 0.05950126 0.18725993 0.02590164 0.56989295 0.03485827
## 7 8 9 10
## 0.02650406 0.05155900 0.11112666 0.69351131Logit模型的估计结果也给出了与Probit模型估计结果相似的似然比和拟合优度指标。但对应的Logit模型为：
\(\widetilde{Prob}\) ( \(GRADE_i\) = 1 | \(x_i\) ) = \(\Lambda (-13.02135+2.8261133*GPA+0.095158*TUCE+2.378688*PSI)\)
Logit模型也叫Logistic模型，服从Logistic分布。Probit模型服从正态分布。两个模型都是离散选择模型的常用模型。但Logit模型简单直接，应用更广。
而且，当因变量是名义变量时，Logit和Probit没有本质的区别，一般情况下可以换用。区别在于采用的分布函数不同，前者假设随机变量服从逻辑概率分布，而后者假设随机变量服从正态分布。其实，这两种分布函数的公式很相似，函数值相差也并不大，唯一的区别在于逻辑概率分布函数的尾巴比正态分布粗一些。但是，如果因变量是序次变量，回归时只能用有序Probit模型。有序Probit可以看作是Logit的扩展。
似然比检验类似于检验模型整体显著性的\(F\)检验，原假设为全部解释变量的系数都为\(0\)，检验的统计量\(LR\)为： \[LR = 2(lnL-lnL_0)\] 其中，\(lnL\)为对概率模型进行\(MLE\)估计的对数似然函数值，\(lnL_0\)为估计只有截距项的模型的对数似然函数值。当原假设成立时，\(LR\)的渐近分布是自由度为\(k-1\)（即除截距项外的解释变量的个数）\(χ^2\)分布。 对于Probit和Logit模型，同样可以计算（上文中）\(CountR^2\)，以反映模型的拟合优度。此外，还可以计算类似于传统\(R^2\)的McFadden似然比指数（McFadden’s Likelihood Ratio Index）来度量拟合优度。似然比指数的定义为：
\[McFadden R^2 = 1 - \frac{\ln{L}}{\ln{L_0}}\]
McFadden \(R^2\) 总是介于\(0\)和\(1\)之间。当所有的斜率系数都为\(0\)时，\(McFadden R^2 = 0\)，但是，McFadden \(R^2\) 不会恰好等于\(1\)。McFadden \(R^2\) 越大，表明拟合得越好。
\(ln{L}\) 是指 Null deviance。
\(ln{L_0}\) 是指 Residual deviance。
library(lmtest)
lrtest(Probit.fit)## Likelihood ratio test
##
## Model 1: Y ~ X1 + X2 + X3
## Model 2: Y ~ 1
## #Df LogLik Df Chisq Pr(>Chisq)
## 1 4 -12.819
## 2 1 -20.592 -3 15.546 0.001405 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1lrtest(Logit.fit)## Likelihood ratio test
##
## Model 1: Y ~ X1 + X2 + X3
## Model 2: Y ~ 1
## #Df LogLik Df Chisq Pr(>Chisq)
## 1 4 -12.890
## 2 1 -20.592 -3 15.404 0.001502 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1上面的Probit.fit检验结果中还给出了有关模型的似然比检验和拟合优度的信息。根据似然比的公式得，L为Model1 LogLik的值，为\(-12.8188\)，L_0为Model2 LogLik的值，为\(-20.5917\)，LR值即为Chisq，为\(15.5458\)，它对应的p值为\(0.0014049\)，因此，它是显著的，表明模型整体是显著的。由于glm()函数并没有直接给出拟合优度的信息，因此需要自己编写一个函数来求解 拟合优度信息值。
上面Logit.fit的检验结果中还给出了有关模型的似然比检验和拟合优度的信息。根据似然比的公式得，L为Model1 LogLik的值，为\(-12.8896\)，L_0为Model2 LogLik的值，为\(-20.5917\)，LR值即为Chisq，为\(15.4042\)，它对应的p值为\(0.0015019\)，因此，它是显著的，表明模型整体是显著的。
Count.Rsquare <- function(glm.object, p){ #p值一般取0.5
num <- length(predict(glm.object, type="response")[predict(glm.object, type="response") > p])
sum_count <- length(predict(glm.object, type="response"))
Count_Rsquare <- num/sum_count
list(Count.Rsquare=Count_Rsquare)
}
Count.Rsquare(Probit.fit, 0.5)## $Count.Rsquare
## [1] 0.34375#Count.Rsquare(Probit.fit, 0.5)$Count.Rsquare
Count.Rsquare(Logit.fit, 0.5)## $Count.Rsquare
## [1] 0.34375McFa.Rsquare <- function(glm.object) {
deviance <- glm.object$deviance
null.deviance <- glm.object$null.deviance
McFa.Rsquare <- 1-(deviance/null.deviance)
list(McFadden.Rsquare = McFa.Rsquare)
#return (McFa.Rsquare)
}
McFa.Rsquare(Probit.fit)## $McFadden.Rsquare
## [1] 0.377478#McFa.Rsquare(Probit.fit)$McFadden.Rsquare
McFa.Rsquare(Logit.fit)## $McFadden.Rsquare
## [1] 0.3740383#连续型变量（非虚拟变量）
#具体每个“自变量”对Y值的影响=平均边际影响*回归系数值
#平均边际影响
Marginal.Effects <- function(glmCoef, modeldata, otype) {
modeldataMean <- apply(modeldata[-1], 2, mean)
n <- length(modeldataMean)
if (otype=='Probit') {
avg.Marginal.Effects <- dnorm(glmCoef[1]+sum(glmCoef[-1]*modeldataMean))
}else if (otype=='Logit') {
avg.Marginal.Effects <- dlogis(glmCoef[1]+sum(glmCoef[-1]*modeldataMean))
}
list(Marginal.Effects = avg.Marginal.Effects)
}
(Probit.effect <- Marginal.Effects(coefficients(Probit.fit), modeldata, 'Probit'))## $Marginal.Effects
## (Intercept)
## 0.3280504(Logit.effect <- Marginal.Effects(coefficients(Logit.fit), modeldata, 'Logit'))## $Marginal.Effects
## (Intercept)
## 0.1889022#“连续型”自变量边际影响=平均边际影响*回归系数值
coef(Probit.fit)[-1]*Probit.effect$Marginal.Effects## X1 X2 X3
## 0.53334835 0.01696954 0.46790834coef(Logit.fit)[-1]*Logit.effect$Marginal.Effects## X1 X2 X3
## 0.53385882 0.01797549 0.44933928#离散型变量（虚拟变量）
#分别计算虚拟变量取1和0时方程式的值，二者的差即为虚拟解释变量的边际影响。
###具体操作过程中，不知道有哪些是虚拟变量，虚拟变量有几个，所以具体情况具体分析！！！
#“离散型”自变量边际影响=平均边际影响(取1)-平均边际影响(取0)Sample.Count.Rsquare <- function(glm.object, testdata, p){
test_prob <- predict(glm.object, newdata=testdata, type="response")
num <- length(test_prob[test_prob > p])
sum_count <- length(test_prob)
Sample_Count_Rsquare <- num/sum_count
list(Sample.Count.Rsquare=Sample_Count_Rsquare)
}
Sample.Count.Rsquare(Probit.fit, modeldata[,-1], 0.5)## $Sample.Count.Rsquare
## [1] 0.34375#Sample.Count.Rsquare(Probit.fit, modeldata[,-1], 0.5)$Sample.Count.Rsquare
Sample.Count.Rsquare(Logit.fit, modeldata[,-1], 0.5)## $Sample.Count.Rsquare
## [1] 0.34375所谓事后预测，所要做的工作是在模型应用上线后，实际发生的结果和模型预测结果相比较，预测成功的比率。
glm.object <- Logit.fit
testdata <- modeldata[,-1] # 测试集x变量
default_test <- modeldata[,1] # 测试集y变量
test_prob <- predict(glm.object, newdata=testdata, type="response")
s <- 0
m <- 0
xulie <- NULL
jilu <- NULL
for(k in 1:32){
if(test_prob[k]>0.5){
a <- 1
}else{
a <- 0
}
if(a==default_test[k]){
s <- s+1
}else{
m <- m+1
xulie[m] <- k
}
xulie
s
}
xulie## [1] 14 19 24 26 31 32s## [1] 26