R Markdown

This is an R Markdown document. Markdown is a simple formatting syntax for authoring HTML, PDF, and MS Word documents. For more details on using R Markdown see http://rmarkdown.rstudio.com.

When you click the Knit button a document will be generated that includes both content as well as the output of any embedded R code chunks within the document. You can embed an R code chunk like this:

summary(cars)
##      speed           dist       
##  Min.   : 4.0   Min.   :  2.00  
##  1st Qu.:12.0   1st Qu.: 26.00  
##  Median :15.0   Median : 36.00  
##  Mean   :15.4   Mean   : 42.98  
##  3rd Qu.:19.0   3rd Qu.: 56.00  
##  Max.   :25.0   Max.   :120.00

Including Plots

You can also embed plots, for example:

Note that the echo = FALSE parameter was added to the code chunk to prevent printing of the R code that generated the plot.

library(haven)
library(sjmisc)
library(sjlabelled)
## 
## Attaching package: 'sjlabelled'
## The following objects are masked from 'package:sjmisc':
## 
##     to_character, to_factor, to_label, to_numeric
## The following objects are masked from 'package:haven':
## 
##     as_factor, read_sas, read_spss, read_stata, write_sas, zap_labels
library(tidyverse)
## -- Attaching packages ------------------------------------------ tidyverse 1.3.0 --
## √ ggplot2 3.3.2     √ purrr   0.3.4
## √ tibble  3.0.3     √ dplyr   1.0.2
## √ tidyr   1.1.2     √ stringr 1.4.0
## √ readr   1.3.1     √ forcats 0.5.0
## -- Conflicts --------------------------------------------- tidyverse_conflicts() --
## x tibble::add_case()       masks sjmisc::add_case()
## x forcats::as_factor()     masks sjlabelled::as_factor(), haven::as_factor()
## x dplyr::as_label()        masks sjlabelled::as_label()
## x dplyr::filter()          masks stats::filter()
## x purrr::is_empty()        masks sjmisc::is_empty()
## x dplyr::lag()             masks stats::lag()
## x sjlabelled::read_sas()   masks haven::read_sas()
## x sjlabelled::read_spss()  masks haven::read_spss()
## x sjlabelled::read_stata() masks haven::read_stata()
## x tidyr::replace_na()      masks sjmisc::replace_na()
## x sjlabelled::write_sas()  masks haven::write_sas()
## x sjlabelled::zap_labels() masks haven::zap_labels()
cgss2017 <- read_dta("C:/Users/huawei/Desktop/R/cgss/cgss2017.dta")
cgss2017backup <- cgss2017
table(cgss2017$a15)
## 
##    1    2    3    4    5   98   99 
##  593 2014 3261 4409 2300    3    2
descr(cgss2017$a8a)
## 
## ## Basic descriptive statistics
## 
##  var    type                        label     n NA.prc   mean      sd       se
##   dd numeric 您个人去年(2016年)全年总收入 12580   0.02 580188 2265003 20194.29
##     md trimmed               range   iqr skew
##  24000  580188 9999999 (0-9999999) 46400 3.91
cgss2017[cgss2017 >= 9999997] <- NA
cgss2017[cgss2017 == 9999996] <- 999999
cgss2017[cgss2017 == 98] <- NA
cgss2017[cgss2017 == 99] <- NA
cgss2017 <- sjlabelled::drop_labels(cgss2017)
cgss2017 <- sjmisc::to_label(cgss2017)
descr(cgss2017$a8a)
## 
## ## Basic descriptive statistics
## 
##  var    type                        label     n NA.prc     mean       sd
##   dd numeric 您个人去年(2016年)全年总收入 11900   5.42 41156.86 233509.9
##       se    md  trimmed               range   iqr  skew
##  2140.58 20000 41156.86 9993600 (0-9993600) 42000 37.51
table(cgss2017$a15)
## 
##   很不健康 比较不健康       一般   比较健康     很健康 
##        593       2014       3261       4409       2300
cgss2017 <- droplevels(cgss2017)
library(tidyverse)
cgssx <-
cgss2017 %>% select(
gender = a2,
independence = a59g,
health = a15,
edu = a7a,
income = a8a,
happiness = a36
)
#cgss2017 %>% select(contains("a"))
#cgss2017 %>% select(a2,a59g,a15,a7a,a8a,a36)
#table(cgssx$edu)
cgssx$education <- as.numeric(cgssx$edu)
cgssx$educ <-
car::recode(cgssx$education, "9:13='3';5:8='2';1;4='1';'14'=NA")
table(cgssx$educ)
## 
##    1    2    3 
## 3511 2326 5196
cgssx$educ <-
factor(cgssx$educ,
levels = c(1, 2, 3),
labels = c("pri", "middle", "high"))
table(cgssx$educ)
## 
##    pri middle   high 
##   3511   2326   5196
cgssx$health <- as.numeric(cgssx$health)
cgssx$health <-
car::recode(cgssx$health, "4:5='3';3='2';1:2='1'")
table(cgssx$health)
## 
##    1    2    3 
## 2607 3261 6709
cgssx$health <-
factor(cgssx$health,
levels = c(1, 2, 3),
labels = c("pri", "middle", "high"))
table(cgssx$health)
## 
##    pri middle   high 
##   2607   3261   6709
cgssx$independence <- as.numeric(cgssx$independence)
cgssx$independence <-
car::recode(cgssx$independence, "1:2='3';3='2';4='1'")
table(cgssx$independence)
## 
##    1    2    3 
##  424 1002 3347
cgssx$independence <-
factor(cgssx$independence,
levels = c(1, 2, 3),
labels = c("pri", "middle", "high"))
table(cgssx$independence)
## 
##    pri middle   high 
##    424   1002   3347
library(gtsummary)
## #BlackLivesMatter
cgssx$income <- as.numeric(cgssx$income)
tbl_summary(cgssx)
Characteristic N = 12,5821
gender
男 5,935 (47%)
女 6,647 (53%)
independence
pri 424 (8.9%)
middle 1,002 (21%)
high 3,347 (70%)
Unknown 7,809
health
pri 2,607 (21%)
middle 3,261 (26%)
high 6,709 (53%)
Unknown 5
edu
没有受过任何教育 1,528 (12%)
私塾、扫盲班 91 (0.7%)
小学 2,726 (22%)
初中 3,511 (28%)
职业高中 163 (1.3%)
普通高中 1,470 (12%)
中专 530 (4.2%)
技校 72 (0.6%)
大学专科（成人高等教育） 377 (3.0%)
大学专科（正规高等教育） 674 (5.4%)
大学本科（成人高等教育） 299 (2.4%)
大学本科（正规高等教育） 948 (7.5%)
研究生及以上 172 (1.4%)
其他 21 (0.2%)
income 20,000 (3,000, 45,000)
Unknown 682
happiness
非常不幸福 216 (1.7%)
比较不幸福 858 (6.8%)
说不上幸福不幸福 1,719 (14%)
比较幸福 7,502 (60%)
非常幸福 2,266 (18%)
Unknown 21
education 4.0 (3.0, 6.0)
educ
pri 3,511 (32%)
middle 2,326 (21%)
high 5,196 (47%)
Unknown 1,549

1 Statistics presented: n (%); Median (IQR) 

tbl_summary(cgssx,statistic = list(all_continuous() ~ "{mean}({sd}) {median} ({min}, {p25}, {p75},{max})"),label = list(income = "income"),digits = list(c(income) ~ 1))
Characteristic N = 12,5821
gender
男 5,935 (47%)
女 6,647 (53%)
independence
pri 424 (8.9%)
middle 1,002 (21%)
high 3,347 (70%)
Unknown 7,809
health
pri 2,607 (21%)
middle 3,261 (26%)
high 6,709 (53%)
Unknown 5
edu
没有受过任何教育 1,528 (12%)
私塾、扫盲班 91 (0.7%)
小学 2,726 (22%)
初中 3,511 (28%)
职业高中 163 (1.3%)
普通高中 1,470 (12%)
中专 530 (4.2%)
技校 72 (0.6%)
大学专科（成人高等教育） 377 (3.0%)
大学专科（正规高等教育） 674 (5.4%)
大学本科（成人高等教育） 299 (2.4%)
大学本科（正规高等教育） 948 (7.5%)
研究生及以上 172 (1.4%)
其他 21 (0.2%)
income 41,156.9(233,509.9) 20,000.0 (0.0, 3,000.0, 45,000.0,9,993,600.0)
Unknown 682
happiness
非常不幸福 216 (1.7%)
比较不幸福 858 (6.8%)
说不上幸福不幸福 1,719 (14%)
比较幸福 7,502 (60%)
非常幸福 2,266 (18%)
Unknown 21
education 5.2(3.3) 4.0 (1.0, 3.0, 6.0,14.0)
educ
pri 3,511 (32%)
middle 2,326 (21%)
high 5,196 (47%)
Unknown 1,549

1 Statistics presented: n (%); Mean(SD) Median (Minimum, IQR,Maximum) 

cgssx %>%
select(gender,independence, income, health, educ) %>% tbl_summary(
by = gender, missing = "no",
type = all_continuous() ~ "continuous2",
statistic = list(
all_continuous() ~ c("{N_nonmiss}",
"{median} ({p25}, {p75})",
"{min}, {max}")
),
label = list(income = "income"),
digits = list( c(income) ~ 1)
) %>%
add_n() %>% # add column with total number of non-missing observations
add_p() %>% # test for a difference between groups
modify_header(label = "**Variable**") %>% # update the column header
bold_labels()
Variable N 男, N = 5,935 女, N = 6,647 p-value1
independence 4,773 0.6
pri 240 (9.0%) 184 (8.7%)
middle 545 (20%) 457 (22%)
high 1,881 (71%) 1,466 (70%)
income 11,900.0 <0.001
N 5,651.0 6,249.0
Median (IQR) 30,000.0 (7,000.0, 50,000.0) 15,000.0 (1,000.0, 36,000.0)
Range 0.0, 9,993,600.0 0.0, 9,930,000.0
health 12,577 <0.001
pri 1,112 (19%) 1,495 (22%)
middle 1,483 (25%) 1,778 (27%)
high 3,336 (56%) 3,373 (51%)
educ 11,033 0.004
pri 1,818 (33%) 1,693 (31%)
middle 1,202 (22%) 1,124 (20%)
high 2,524 (46%) 2,672 (49%)

1 Statistical tests performed: chi-square test of independence; Wilcoxon rank-sum test 

library(table1)
## 
## Attaching package: 'table1'
## The following objects are masked from 'package:base':
## 
##     units, units<-
table1(~ gender+edu+income+independence+health+educ, data=cgssx)
Overall
(N=12582)
gender
男 5935 (47.2%)
女 6647 (52.8%)
edu
没有受过任何教育 1528 (12.1%)
私塾、扫盲班 91 (0.7%)
小学 2726 (21.7%)
初中 3511 (27.9%)
职业高中 163 (1.3%)
普通高中 1470 (11.7%)
中专 530 (4.2%)
技校 72 (0.6%)
大学专科（成人高等教育） 377 (3.0%)
大学专科（正规高等教育） 674 (5.4%)
大学本科（成人高等教育） 299 (2.4%)
大学本科（正规高等教育） 948 (7.5%)
研究生及以上 172 (1.4%)
其他 21 (0.2%)
income
Mean (SD) 41200 (234000)
Median [Min, Max] 20000 [0, 9990000]
Missing 682 (5.4%)
independence
pri 424 (3.4%)
middle 1002 (8.0%)
high 3347 (26.6%)
Missing 7809 (62.1%)
health
pri 2607 (20.7%)
middle 3261 (25.9%)
high 6709 (53.3%)
Missing 5 (0.0%)
educ
pri 3511 (27.9%)
middle 2326 (18.5%)
high 5196 (41.3%)
Missing 1549 (12.3%) 
table1(~ edu+income+educ|gender, data=cgssx)
男
(N=5935)		 女
(N=6647)		 Overall
(N=12582)
edu
没有受过任何教育 383 (6.5%) 1145 (17.2%) 1528 (12.1%)
私塾、扫盲班 38 (0.6%) 53 (0.8%) 91 (0.7%)
小学 1258 (21.2%) 1468 (22.1%) 2726 (21.7%)
初中 1818 (30.6%) 1693 (25.5%) 3511 (27.9%)
职业高中 85 (1.4%) 78 (1.2%) 163 (1.3%)
普通高中 768 (12.9%) 702 (10.6%) 1470 (11.7%)
中专 264 (4.4%) 266 (4.0%) 530 (4.2%)
技校 47 (0.8%) 25 (0.4%) 72 (0.6%)
大学专科（成人高等教育） 205 (3.5%) 172 (2.6%) 377 (3.0%)
大学专科（正规高等教育） 338 (5.7%) 336 (5.1%) 674 (5.4%)
大学本科（成人高等教育） 147 (2.5%) 152 (2.3%) 299 (2.4%)
大学本科（正规高等教育） 486 (8.2%) 462 (7.0%) 948 (7.5%)
研究生及以上 90 (1.5%) 82 (1.2%) 172 (1.4%)
其他 8 (0.1%) 13 (0.2%) 21 (0.2%)
income
Mean (SD) 48200 (204000) 34800 (257000) 41200 (234000)
Median [Min, Max] 30000 [0, 9990000] 15000 [0, 9930000] 20000 [0, 9990000]
Missing 284 (4.8%) 398 (6.0%) 682 (5.4%)
educ
pri 1818 (30.6%) 1693 (25.5%) 3511 (27.9%)
middle 1202 (20.3%) 1124 (16.9%) 2326 (18.5%)
high 2524 (42.5%) 2672 (40.2%) 5196 (41.3%)
Missing 391 (6.6%) 1158 (17.4%) 1549 (12.3%) 
table1(~ income|gender*educ, data=cgssx)
男
女
Overall
pri
(N=1818)		 middle
(N=1202)		 high
(N=2524)		 pri
(N=1693)		 middle
(N=1124)		 high
(N=2672)		 pri
(N=3511)		 middle
(N=2326)		 high
(N=5196)
income
Mean (SD) 32000 (43400) 55100 (299000) 62400 (232000) 24200 (113000) 43700 (311000) 48900 (339000) 28300 (84400) 49600 (305000) 55500 (292000)
Median [Min, Max] 25000 [0, 960000] 36000 [0, 9930000] 30000 [0, 9990000] 20000 [0, 4420000] 30000 [0, 9930000] 20000 [0, 9930000] 20000 [0, 4420000] 30000 [0, 9930000] 24000 [0, 9990000]
Missing 59 (3.2%) 62 (5.2%) 143 (5.7%) 90 (5.3%) 75 (6.7%) 162 (6.1%) 149 (4.2%) 137 (5.9%) 305 (5.9%) 
table1(~ income|gender*educ, data=cgssx,overall=FALSE)
男
女
pri
(N=1818)		 middle
(N=1202)		 high
(N=2524)		 pri
(N=1693)		 middle
(N=1124)		 high
(N=2672)
income
Mean (SD) 32000 (43400) 55100 (299000) 62400 (232000) 24200 (113000) 43700 (311000) 48900 (339000)
Median [Min, Max] 25000 [0, 960000] 36000 [0, 9930000] 30000 [0, 9990000] 20000 [0, 4420000] 30000 [0, 9930000] 20000 [0, 9930000]
Missing 59 (3.2%) 62 (5.2%) 143 (5.7%) 90 (5.3%) 75 (6.7%) 162 (6.1%) 
library(ggstatsplot)
## Registered S3 method overwritten by 'broom.mixed':
##   method      from 
##   tidy.gamlss broom
## Registered S3 methods overwritten by 'lme4':
##   method                          from
##   cooks.distance.influence.merMod car 
##   influence.merMod                car 
##   dfbeta.influence.merMod         car 
##   dfbetas.influence.merMod        car
## In case you would like cite this package, cite it as:
##      Patil, I. (2018). ggstatsplot: "ggplot2" Based Plots with Statistical Details. CRAN.
##      Retrieved from https://cran.r-project.org/web/packages/ggstatsplot/index.html
ggstatsplot::ggpiestats(
data = cgssx,
x = educ,
y = gender,
title = "教育与性别",
legend.title = "性别",
caption = substitute(paste(italic("Source"), ": 1974 Motor Trend US magazine"))
)

ggstatsplot::ggbarstats(
data = cgssx,
x = educ,
y = gender,
title = "Education",
xlab = "性别",
legend.title = "gender"
#ggtheme = hrbrthemes::theme_ipsum_pub(),
)

ggstatsplot::ggbarstats(
data = cgssx,
x = gender,
y = health,
title = "Health",
xlab = "性别",
legend.title = "health"
#ggtheme = hrbrthemes::theme_ipsum_pub(),
)

ggstatsplot::ggpiestats(
data = cgssx,
x = income,
y = health,
title = "收入与健康",
legend.title = "收入",
caption = substitute(paste(italic("Source"), ": 1974 Motor Trend US magazine"))
)
## Warning in stats::chisq.test(x = x_arg, correct = FALSE): Chi-squared
## approximation may be incorrect
## Warning:  Number of labels is greater than default palette color count.
##  Try using another color `palette` (and/or `package`).
## 

descr(cgssx$income)
## 
## ## Basic descriptive statistics
## 
##  var    type label     n NA.prc     mean       sd      se    md  trimmed
##   dd numeric    dd 11900   5.42 41156.86 233509.9 2140.58 20000 24005.94
##                range   iqr  skew
##  9993600 (0-9993600) 42000 37.51
cgssx$income1<-log(cgssx$income+1)
descr(cgssx$income1)
## 
## ## Basic descriptive statistics
## 
##  var    type label     n NA.prc mean   sd   se  md trimmed           range
##   dd numeric    dd 11900   5.42 8.35 3.85 0.04 9.9    8.96 16.12 (0-16.12)
##       iqr  skew
##  2.707739 -1.47
c + geom_dotplot(binwidth = 500)+scale_x_continuous(limits = c(0, 999999))
## NULL
c <- ggplot(cgssx[is.na(cgssx$income) == F,], aes(income))
c + geom_dotplot(binwidth = 500)+scale_x_continuous(limits = c(0, 999999))
## Warning: Removed 19 rows containing non-finite values (stat_bindot).

c + geom_histogram()+ylim(0,1500)+xlim(0, 999999)
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## Warning: Removed 19 rows containing non-finite values (stat_bin).
## Warning: Removed 3 rows containing missing values (geom_bar).

cgssx$income1<-log(cgssx$income+1)
descr(cgssx$income1)
## 
## ## Basic descriptive statistics
## 
##  var    type label     n NA.prc mean   sd   se  md trimmed           range
##   dd numeric    dd 11900   5.42 8.35 3.85 0.04 9.9    8.96 16.12 (0-16.12)
##       iqr  skew
##  2.707739 -1.47
d <- ggplot(cgssx[is.na(cgssx$income1) == F,], aes(income1))
d + geom_dotplot()+scale_x_continuous(limits = c(0, 14))
## `stat_bindot()` using `bins = 30`. Pick better value with `binwidth`.
## Warning: Removed 10 rows containing non-finite values (stat_bindot).

d + geom_histogram()+ylim(0,2000)+xlim(0, 14)+ scale_x_continuous(
breaks = c(2, 4, 6,8,10,12,14),
label = c("two", "four", "six","eight","ten","twelve","fourteen")
)
## Scale for 'x' is already present. Adding another scale for 'x', which will
## replace the existing scale.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## Warning: Removed 1 rows containing missing values (geom_bar).

ggplot(cgssx[is.na(cgssx$happiness) == F,], aes(x=happiness)) +
geom_bar(stat="count",position = "dodge",na.rm = TRUE)+geom_text(stat='count',aes(label=..count..),vjust=-1,position=position_dodge(width=0.9))

ggplot(cgssx[is.na(cgssx$happiness) == F,], aes(x=happiness,fill=educ,na.rm = TRUE)) +
geom_bar(stat="count",na.rm = TRUE)+geom_text(stat='count',aes(label=..count..),position=position_stack(vjust=0.5))+
facet_wrap(~gender)

ggplot(cgssx[is.na(cgssx$happiness) == F,], aes(x=happiness,fill=income1)) +geom_bar(stat="count",position = "dodge")+geom_text(stat='count',aes(label=..count..),vjust=-1,position=position_dodge(width=0.9))

ggplot(cgssx[is.na(cgssx$happiness) == F,], aes(x=happiness,fill=independence)) +geom_bar(stat="count",position = "dodge")+geom_text(stat='count',aes(label=..count..),vjust=-1,position=position_dodge(width=0.9))

c <- ggplot(mpg, aes(hwy))

descr(cgssx$income)
## 
## ## Basic descriptive statistics
## 
##  var    type label     n NA.prc     mean       sd      se    md  trimmed
##   dd numeric    dd 11900   5.42 41156.86 233509.9 2140.58 20000 24005.94
##                range   iqr  skew
##  9993600 (0-9993600) 42000 37.51
cgssx%>%filter(!is.na(educ)&!is.na(happiness))%>%
ggplot(.,aes(x=happiness,fill=independence)) +
geom_bar(stat="count",position = "dodge")+geom_text(stat='count',aes(label=..count..),vjust=-1,position=position_dodge(width=0.9))+theme(axis.text.x = element_text(size = 8, color = "blue", face = "bold", vjust = 0.5, hjust = 0.5, angle = 15))

cgssx%>%filter(!is.na(educ)&!is.na(happiness))%>%
ggplot(.,aes(x=happiness,fill=health)) +
geom_bar(stat="count",position = "dodge")+geom_text(stat='count',aes(label=..count..),vjust=-1,position=position_dodge(width=0.9))+theme(axis.text.x = element_text(size = 8, color = "blue", face = "bold", vjust = 0.5, hjust = 0.5, angle = 15))

library(emmeans)
library(effects)
## Loading required package: carData
## lattice theme set by effectsTheme()
## See ?effectsTheme for details.
library(ggeffects)
names(cgssx)
## [1] "gender"       "independence" "health"       "edu"          "income"      
## [6] "happiness"    "education"    "educ"         "income1"
sat_mod1 <- lm(happiness ~ gender+educ, # regression formula
data = cgssx) # data
## Warning in model.response(mf, "numeric"): using type = "numeric" with a factor
## response will be ignored
## Warning in Ops.factor(y, z$residuals): '-' not meaningful for factors
sat_mod2<-lm(happiness ~ gender+educ +health , data = cgssx)
## Warning in model.response(mf, "numeric"): using type = "numeric" with a factor
## response will be ignored

## Warning in model.response(mf, "numeric"): '-' not meaningful for factors
sat_mod3<-lm(happiness ~ gender+educ +health +educ*gender, data = cgssx)
## Warning in model.response(mf, "numeric"): using type = "numeric" with a factor
## response will be ignored

## Warning in model.response(mf, "numeric"): '-' not meaningful for factors
mydf<-ggpredict(sat_mod3, terms = c("health", "educ"),data = cgssx)
## Warning in Ops.factor(r, 2): '^' not meaningful for factors
ggplot(mydf, aes(x = x, y = predicted, colour = group)) +
geom_point() +
stat_smooth(method = "lm", se = FALSE)
## `geom_smooth()` using formula 'y ~ x'

mydf1<-ggpredict(sat_mod3, terms = c("educ", "health"),data = cgssx)
## Warning in Ops.factor(r, 2): '^' not meaningful for factors
ggplot(mydf1, aes(x = x, y = predicted, colour = group)) +
geom_point() +
stat_smooth(method = "lm", se = FALSE)
## `geom_smooth()` using formula 'y ~ x'

mydf2<-ggpredict(sat_mod3, terms = c("health","educ"),data = cgssx)
## Warning in Ops.factor(r, 2): '^' not meaningful for factors
ggplot(mydf2, aes(x = x, y = predicted, colour = group)) +
geom_line() +
stat_smooth(method = "lm", se = FALSE)
## `geom_smooth()` using formula 'y ~ x'
## geom_path: Each group consists of only one observation. Do you need to adjust
## the group aesthetic?

library(modelsummary)
models <- list(
sat_mod1,
sat_mod2 ,
sat_mod3
)
modelsummary(models,stars = TRUE, stars_note = TRUE)
## Warning in Ops.factor(r, 2): '^' not meaningful for factors
## Warning in Ops.factor(r, 2): '^' not meaningful for factors

## Warning in Ops.factor(r, 2): '^' not meaningful for factors

## Warning in Ops.factor(r, 2): '^' not meaningful for factors

## Warning in Ops.factor(r, 2): '^' not meaningful for factors
## Warning in Ops.factor(res, 2): '^' not meaningful for factors

## Warning in Ops.factor(res, 2): '^' not meaningful for factors

## Warning in Ops.factor(res, 2): '^' not meaningful for factors
## Warning in Ops.factor(weighted.residuals(object), 2): '^' not meaningful for
## factors
## Warning in Ops.factor(r, 2): '^' not meaningful for factors

## Warning in Ops.factor(r, 2): '^' not meaningful for factors
## Warning in Ops.factor(res, 2): '^' not meaningful for factors

## Warning in Ops.factor(res, 2): '^' not meaningful for factors

## Warning in Ops.factor(res, 2): '^' not meaningful for factors
## Warning in Ops.factor(weighted.residuals(object), 2): '^' not meaningful for
## factors
## Warning in Ops.factor(r, 2): '^' not meaningful for factors

## Warning in Ops.factor(r, 2): '^' not meaningful for factors
## Warning in Ops.factor(res, 2): '^' not meaningful for factors

## Warning in Ops.factor(res, 2): '^' not meaningful for factors

## Warning in Ops.factor(res, 2): '^' not meaningful for factors
## Warning in Ops.factor(weighted.residuals(object), 2): '^' not meaningful for
## factors
Model 1 Model 2 Model 3
(Intercept) 3.818NA 3.505NA 3.503NA
gender女 0.094NA 0.101NA 0.105NA
educmiddle 0.053NA 0.031NA 0.014NA
educhigh 0.019NA 0.019NA 0.031NA
healthmiddle 0.225NA 0.225NA
healthhigh 0.452NA 0.452NA
gender女 × educmiddle 0.034NA
gender女 × educhigh -0.023NA
Num.Obs. 11017 11013 11013
R2
R2 Adj.
AIC
BIC
Log.Lik.
F
* p < 0.1, ** p < 0.05, *** p < 0.01 
m1 <- glm(happiness ~ health+educ + gender+income1+independence,
data = cgssx,
family = binomial(link = "probit"))
m2 <- glm(happiness ~ health+educ + gender+income1+independence+income1*educ ,
data = cgssx,
family = binomial(link = "probit"))
models1 <- list(
m1,
m2
)
modelsummary(models1,stars = TRUE, stars_note = TRUE)
Model 1 Model 2
(Intercept) 0.849*** 0.911**
(0.282) (0.430)
healthmiddle 0.661*** 0.660***
(0.183) (0.183)
healthhigh 0.823*** 0.822***
(0.157) (0.157)
educmiddle 0.114 0.189
(0.168) (0.619)
educhigh 0.092 -0.072
(0.141) (0.504)
gender女 0.093 0.094
(0.127) (0.128)
income1 0.042* 0.035
(0.022) (0.040)
independencemiddle 0.381* 0.383*
(0.202) (0.201)
independencehigh 0.384** 0.385**
(0.164) (0.164)
educmiddle × income1 -0.008
(0.062)
educhigh × income1 0.017
(0.050)
Num.Obs. 4397 4397
AIC 435.1 438.9
BIC 492.6 509.1
Log.Lik. -208.542 -208.431
* p < 0.1, ** p < 0.05, *** p < 0.01 
modelsummary(models1,stars = TRUE, stars_note = TRUE,exponentiate = TRUE)
Model 1 Model 2
(Intercept) 2.337*** 2.486**
(0.282) (0.430)
healthmiddle 1.937*** 1.934***
(0.183) (0.183)
healthhigh 2.278*** 2.275***
(0.157) (0.157)
educmiddle 1.121 1.208
(0.168) (0.619)
educhigh 1.096 0.931
(0.141) (0.504)
gender女 1.097 1.098
(0.127) (0.128)
income1 1.042* 1.036
(0.022) (0.040)
independencemiddle 1.464* 1.466*
(0.202) (0.201)
independencehigh 1.468** 1.469**
(0.164) (0.164)
educmiddle × income1 0.992
(0.062)
educhigh × income1 1.017
(0.050)
Num.Obs. 4397 4397
AIC 435.1 438.9
BIC 492.6 509.1
Log.Lik. -208.542 -208.431
* p < 0.1, ** p < 0.05, *** p < 0.01 
library(ordinal)
## 
## Attaching package: 'ordinal'
## The following object is masked from 'package:dplyr':
## 
##     slice
m3 <- clm(happiness ~ educ + gender+income,
          data = cgssx,na.action=na.omit)
## Warning: (2) Model is nearly unidentifiable: very large eigenvalue
##  - Rescale variables? 
## In addition: Absolute and relative convergence criteria were met
m4 <- clm(happiness ~educ + gender+income+income*educ ,
          data = cgssx,na.action=na.omit)
## Warning: (2) Model is nearly unidentifiable: very large eigenvalue
##  - Rescale variables? 
## In addition: Absolute and relative convergence criteria were met
modelsc <- list( m3,m4)
modelsummary(modelsc,stars = TRUE, coef_omit = "幸福",stars_note = TRUE)
Model 1 Model 2
educmiddle 0.096* 0.193***
(0.055) (0.068)
educhigh 0.021 0.156***
(0.044) (0.053)
gender女 0.215*** 0.231***
(0.039) (0.039)
income 0.000*** 0.000***
(0.000) (0.000)
educmiddle × income 0.000***
(0.000)
educhigh × income 0.000***
(0.000)
Num.Obs. 10429 10429
AIC 23222.3 23200.9
BIC 23280.3 23273.4
Log.Lik. -11603.1606717287 -11590.4387995499
edf 8 10
* p < 0.1, ** p < 0.05, *** p < 0.01