幸福指数相关影响因素研究

数据集

mdata <- read.csv('http://data.galaxystatistics.com/blog_data/regression/happiness_index_data.csv',header=T)
head(mdata)

##   happy money sex love work
## 1    10    36   0    3    4
## 2     8    47   1    3    1
## 3     8    53   0    3    5
## 4     8    35   1    3    3
## 5     4    88   1    1    2
## 6     9   175   1    3    4

modeldata <- mdata
# 描述性统计计算
summary(modeldata)

##      happy           money             sex            love     
##  Min.   : 4.00   Min.   : 35.00   Min.   :0.00   Min.   :1.00  
##  1st Qu.: 5.00   1st Qu.: 40.00   1st Qu.:1.00   1st Qu.:2.00  
##  Median : 8.00   Median : 45.00   Median :1.00   Median :3.00  
##  Mean   : 6.96   Mean   : 65.32   Mean   :0.76   Mean   :2.52  
##  3rd Qu.: 8.00   3rd Qu.: 81.00   3rd Qu.:1.00   3rd Qu.:3.00  
##  Max.   :10.00   Max.   :175.00   Max.   :1.00   Max.   :3.00  
##       work     
##  Min.   :1.00  
##  1st Qu.:3.00  
##  Median :4.00  
##  Mean   :3.32  
##  3rd Qu.:4.00  
##  Max.   :5.00

#数据集字段名称转换
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 X4
## 1 10  36  0  3  4
## 2  8  47  1  3  1
## 3  8  53  0  3  5
## 4  8  35  1  3  3
## 5  4  88  1  1  2
## 6  9 175  1  3  4

class(modeldata)

## [1] "data.frame"

dim(modeldata)

## [1] 25  5

regFormula <- as.formula(paste("Y ~ ", paste(names(modeldata)[-1], collapse= "+")))
regFormula

## Y ~ X1 + X2 + X3 + X4

# save.image()
getwd()

## [1] "D:/BaiduYunDisk/BaiduYunDisk/galaxystatistics/galaxystatistics/银河统计-Statistics_Method/report/regression/article8"

有序Logit模型1

#有序Logit模型和有序Probit模型可以用MASS包里的polr()函数拟合：
library(MASS)
# method设为logistic
Logit.fit <- polr(regFormula, data=modeldata, method="logistic")
#展示拟合模型的详细结果
summary(Logit.fit)
#列出拟合模型的模型参数
coefficients(Logit.fit) 
head(predict(Logit.fit, modeldata[-1], type='probs'),10)

有序Logit模型2

# load('.RData')

# 有序Logit模型也可以用rms包中的lrm()函数拟合：
# install.packages('rms')
library(rms)
Logit.fit1 <- lrm(regFormula, data=modeldata)
Logit.fit1

## 
## Logistic Regression Model
## 
## lrm(formula = regFormula, data = modeldata)
## Frequencies of Responses
## 
##  4  5  6  7  8  9 10 
##  4  3  2  2 10  3  1 
## 
##                       Model Likelihood     Discrimination    Rank Discrim.    
##                          Ratio Test            Indexes          Indexes       
## Obs            25    LR chi2      35.37    R2       0.783    C       0.911    
## max |deriv| 6e-09    d.f.             4    g        4.727    Dxy     0.822    
##                      Pr(> chi2) <0.0001    gr     112.947    gamma   0.839    
##                                            gp       0.450    tau-a   0.660    
##                                            Brier    0.064                     
## 
##       Coef     S.E.   Wald Z Pr(>|Z|)
## y>=5  -10.3732 3.3822 -3.07  0.0022  
## y>=6  -12.8730 3.7414 -3.44  0.0006  
## y>=7  -14.4137 4.0341 -3.57  0.0004  
## y>=8  -15.4344 4.2140 -3.66  0.0002  
## y>=9  -19.4267 4.9594 -3.92  <0.0001 
## y>=10 -21.2577 5.0879 -4.18  <0.0001 
## X1      0.0109 0.0120  0.91  0.3622  
## X2     -2.7295 1.2532 -2.18  0.0294  
## X3      5.5087 1.5002  3.67  0.0002  
## X4      0.8215 0.6044  1.36  0.1741

head(predict(Logit.fit1, modeldata[-1]),10)

##           1           2           3           4           5           6 
##  9.83089499  4.75665690 10.83753411  6.26899436 -4.99267149  8.61510558 
##           7           8           9          10 
##  8.61510558  3.59873235 -0.06117843 -3.70902439

# 可见对于有序的Logit模型，lrm()函数拟合的结果与polr()函数拟合的结果相同。并且lrm()函数还给出了似然比检验结果。

有序Probit模型

#method设为probit
library(MASS)
Probit.fit <- polr(regFormula, data=modeldata, method="probit")
#展示拟合模型的详细结果
summary(Probit.fit)
#列出拟合模型的模型参数
coefficients(Probit.fit) 
head(predict(Probit.fit, modeldata[-1], type='probs'),10)

幸福指数相关影响因素研究

baorui

摘要

数据集

有序Logit模型1

有序Logit模型2

有序Probit模型