---
title: "Tourist2023cn"
author: "Kai Wu"
date: "2023/5/5"
output: html_document
---

游客满意度调数据（模拟数据，42变量，373条记录）

这是一个随机数模拟的数据，可以练习使用，但是分析结论同学不要太认真

（1）电子问卷

http://kaiwu.city/openfiles/tourist_satisfaction_questionnaire_cn.docx
  

（2）csv格式调研数据

http://kaiwu.city/openfiles/tourist.csv

（3）Excel格式调研数据

http://kaiwu.city/openfiles/tourist_cn.xlsx

（4）SPSS格式调研数据
http://kaiwu.city/openfiles/data_tourist_cn.sav

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```

# 1.data preparation数据预处理

data preparation（数据准备）与data wrangling、data munging、data
cleaning的含义非常接近
https://theappsolutions.com/blog/development/data-wrangling-guide-to-data-preparation/

Data Wrangling in R 
(1)Dplyr - essential data-munging R package. Supreme
data framing tool. Especially useful for data management operating by
categories. 
(2)Purrr - good for list function operations and
error-checking. 
(3)Splitstackshape - an oldie but goldie. Good for
shaping complex data sets and simplifying the visualization.
(4)JSOnline - nice and easy parsing tool.
(5)Magrittr - good for
wrangling scattered sets and putting them into a more coherent form.



## 1.1 load packages载入拓展功能包

```{r}

# 最为稳健的方法，检测是否安装，没有安装就安装

if(!isTRUE(require("Hmisc"))){install.packages("Hmisc")}
#数据分析包
if(!isTRUE(require("psych"))){install.packages("psych")}
#数据分析包
if(!isTRUE(require("lavaan"))){install.packages("lavaan")} #结构方程分析

if(!isTRUE(require("broom"))){install.packages("broom")}# 分析结果整洁输出

if(!isTRUE(require("ggplot2"))){install.packages("ggplot2")}# 专业绘图，数据可视化
if(!isTRUE(require("ggrepel"))){install.packages("ggrepel")} # 绘图的标注（不遮挡）
if(!isTRUE(require("corrplot"))){install.packages("corrplot")} #相关分析绘图
if(!isTRUE(require("semPlot"))){install.packages("semPlot")} #结构方程分析绘图


if(!isTRUE(require("leaflet"))){install.packages("leaflet")}# 地图绘制，基于openstreet数据


#library(Hmisc) #数据分析包
#library(psych) #数据分析包
#library(broom) # 分析结果整洁输出
#library(leaflet) # 地图绘制，基于openstreet数据
#library(ggplot2) #专业绘图
#library(ggrepel) # 绘图的标注（不遮挡）
#library(lavaan) #结构方程分析
#library(semPlot) #结构方程分析绘图
#library(corrplot)# 相关系数图

```

## 1.2 import dataset导入数据

### 1.2a设定工作文件夹

```{r}

#设定工作文件夹, 根据自己电脑的情况进行修改
datafolder<-'D:/tdata/tourist_cn_R/'



#文件夹的设定，注意斜线的方向
#datafolder<-'D:\\tdata\\tourist_cn_R\\'


#数据可以通过网盘直链下载
# http://kaiwu.city/openfiles/tourist.csv 
#373条记录（被调查者），42个变量
#数据下载后存入datafolder

```

### 1.2b导入csv格式数据

数据可以通过网盘直链下载
http://kaiwu.city/openfiles/tourist.csv
373条记录（被调查者），42个变量 数据下载后存入datafolder

read.csv是常见的读入数据方法（csv格式，以逗号为分隔符的数据集）

```{r}
tourist<-read.csv('http://kaiwu.city/openfiles/tourist.csv',header=TRUE,encoding = 'UTF-8')



#tourist<-read.csv(paste(datafolder,'tourist.csv',sep=""),header=TRUE,encoding = 'UTF-8')

#tourist<-read.csv('D:/tdata/tourist_cn_R/tourist.csv',header=TRUE,encoding = 'UTF-8')


#3种相同效果的文件名及路径
#第1种：'D:/pdata/rdata/tourist_cn/tourist.csv'
#第2种：paste(datafolder,'tourist.csv',sep="")
#第3种：paste0(datafolder,'tourist.csv')
#注意第3种paste0,最后1位是数字0，不是字母。


#header表示表头——数据中包含变量名的，header==TRUE，否则header==FALSE
#encoding = 'UTF-8'是csv文件的编码格式，UTF-8支持中文等多种格式，但是注意excel另存为csv时，默认是ANSI格式，需要注意设置一下。

#head 表示呈现前6行数据，tail显示后6行数据
head(tourist)
#tail(tourist)



```

### 1.2c 练习1 导入csv格式数据

```{r}
# 练习1 ： 下载数据，试着用本地磁盘文件，修改datafolder目录，读入数据。
#tourist<-read.csv('http://kaiwu.city/openfiles/tourist.csv',header=TRUE,encoding = 'UTF-8')

#tourist<-read.csv(paste(datafolder,'tourist.csv',sep=""),header=TRUE,encoding = 'UTF-8')






```

### 1.2d导入SPSS格式数据

SPSS格式数据下载网址
http://kaiwu.city/openfiles/data_tourist_cn.sav
注意要安装haven拓展功能包，才可以导入

```{r}

#导入SPSS格式数据也可以，需要haven
#if(!isTRUE(require("haven"))){install.packages("haven")}
#tourist_spss<-read_sav(url("http://kaiwu.city/openfiles/data_tourist_cn.sav"))

#tourist_spss<-read_sav(paste0(datafolder,'data_tourist_satisfaction_cnlabels.sav'))
#head(tourist_spss)


#head(tourist_spss)


```




### 1.2e导入excel格式数据

如果想直接导入excel格式数据(xls或xlsx)，推荐使XLConnect
XLConnect可以实现Excel与R的数据双向交互

连接汇报用的excel文件 参考博文
https://www.r-bloggers.com/2016/03/few-steps-to-connect-r-with-excel-xlconnect-2/
http://www.sthda.com/english/articles/2-r/4-xlconnect-read-write-and-manipulate-microsoft-excel-files-from-within-r/

#导入Excel格式数据也可以，需要readxl

```{r}

#if(!isTRUE(require("readxl"))){install.packages("readxl")}

#tourist_excel<-read_excel(paste0(datafolder,"data_tourist_satisfaction_en.xlsx"),sheet="data")


```



#导入Excel格式数据也可以，需要openxlsx
```{r}


#if(!isTRUE(require("openxlsx"))){install.packages("openxlsx")}

#tourist_excel<- read.xlsx("http://kaiwu.city/openfiles/tourist_cn.xlsx",sheet="data")


#tourist_excel<- read.xlsx(paste0(datafolder,"data_tourist_satisfaction_en.xlsx"),sheet="data")


```


## 1.3 value labels变量值标签

### 1.3a 变量值标签：示例

```{r}
#value labels变量值标签，参考电子问卷
#   http://kaiwu.city/openfiles/tourist_satisfaction_questionnaire_cn.docx



#tourist是数据框dataframe的名字，gendr是变量的名字，tourist$gender就是调用tourist数据框的gender变量

#等同的效果是tourist[,"gender"]或tourist[,2]

#factor是R的一种数据形式，可以简单理解为分类变量，注意这个levels是可以有顺序的。

tourist$gender <- factor(tourist$gender,levels = c(1,2),labels = c("男","女"))


tourist$thotel <- factor(tourist$thotel,levels = c(1,2,3,4),labels = c("经济型酒店", "豪华型酒店", "民宿", "酒店式公寓"))

tourist$type3 <- factor(tourist$type3,levels = c(1,2,3),labels = c("自然风光", "历史文化", "自然与历史混合"))

tourist$type2 <- factor(tourist$type2,levels = c(1,2),labels = c("观赏型", "参与型"))

head(tourist)


```

### 1.3b 练习2 变量值标签

```{r}
##### 练习2
##### region（区域）变量的变量值标签。
##### 1	2	3	4	5	6	7
##### 华中	华东	华北	东北	西北	西南	华西


```

```{r}
##### 练习2的答案


tourist$region <- factor(tourist$region,levels = c(1,2,3,4,5,6,7),labels = c('华中',	'华东',	'华北',	'东北',	'西北',	'西南',	'华西'))

head(tourist)
```

## 1.4 compute variables 计算变量

### 1.4a 计算变量：示例

计算变量，compute variables与recode是两个思路，但是有时可以互相替换。

```{r}
#计算总体满意度sat
tourist<-transform(tourist,sat=(sat1+sat2+sat3+sat4+sat5+sat6)/6)

summary(tourist$sat)

```

### 1.4b 练习3 计算变量

```{r}
#compute variables
# 练习3 ： 基于byear（出生年份），计算age变量（年龄）。
# 2022- byear





```

```{r}
#练习3答案
#练习3 ： 基于byear（出生年份），计算age变量（年龄）。
# 方案1：transform
tourist<-transform(tourist,age=2023-byear)

# 方案2：直接计算
# tourist$age<-2022-tourist$byear

summary(tourist$age)


```

## 1.5 recode 重新编码

### 1.5a recode:收入段

attach和detach配合，可以减少tourist\$的使用 []里面是条件

```{r}
# recode 重新编码：income3收入段变量——分3段
attach(tourist)
tourist$income3[income  < 2000] <- 1
tourist$income3[income >= 2000 & income  <3000] <- 2
tourist$income3[income >= 3000] <- 3
detach(tourist)

# 加变量值标签
tourist$income3 <- factor(tourist$income3,levels = c(1,2,3),labels = c("2000以下", "2000-3000", "3000以上"))


summary(tourist$income3)
```

### 1.5b 练习4 recode

```{r}
#compute variables
# 练习4 ： 把收入变量income分为5段，2000以下，2000-2499，2500-2999，3000-3499，3500以上




```

```{r}
# 练习4答案
# 练习4 ： 把收入变量income分为5段，2000以下，2000-2499，2500-2999，3000-3499，3500以上

attach(tourist)
tourist$income5[income  < 2000] <- 1
tourist$income5[income >= 2000 & income  <2500] <- 2
tourist$income5[income >= 2500 & income  <3000] <- 3
tourist$income5[income >= 3000 & income  <3500] <- 4
tourist$income5[income >= 3500] <- 5
detach(tourist)

# 加变量值标签
tourist$income5 <- factor(tourist$income5,levels = c(1,2,3,4,5),labels = c("2000以下", "2000-2499","2500-2999","3000-3499", "3500以上"))


summary(tourist$income5)


```

### 1.5c recode:旅游支出分段

```{r}

# recode 重新编码：旅游花费变量——分3段
attach(tourist)
tourist$expense3[expense  < 300] <- 1
tourist$expense3[expense >= 300 & expense  < 400] <- 2
tourist$expense3[expense > 400] <- 3
detach(tourist)


# 加变量值标签
tourist$expense3 <- factor(tourist$expense3,levels = c(1,2,3),labels = c("300以下", "300-399", "400以上"))
summary(tourist$expense3)

```

## 1.6 summary()

快速浏览、汇总（最大值、最小值、数据类型、变量值标签） 有利于发现异常值

```{r}

summary(tourist)

```

## 1.7 数据的分布情况

数据分布对于统计方法选择，统计检验非常重要
为了提高运行效率，使用R默认的绘图工具，没有用ggplot2

### 1.7a 数据分布：示例

```{r}
#概率分布图

plot(density(tourist$sat1))
#main=""是图标题
#xlab='总体满意度'设置x轴标题
#ylab='频数'设置y轴标题

plot(density(tourist$sat1),main="", xlab='景区满意度', ylab='频数')

hist(tourist$sat1)
hist(tourist$sat1,main="", xlab='景区满意度', ylab='频数')

```

### 1.7b 练习5 数据分布

```{r}
#查看概率分布
# 练习5 ： 查看收入变量income的分布。


```

```{r}
# 练习5答案
# 练习5 ： 查看收入变量income的分布。
hist(tourist$income,main="", xlab='收入', ylab='频数')

```

### 1.7c 用循环语句查看多个变量的分布：密度图

程序设计可以提高批处理的效率，但是我们不做要求，感兴趣同学自己研究一下

```{r}
#10是变量所在的列的序号
#for (i in 10:40){
#   plot(density(tourist[,i]))
#   plot(density(tourist[,i]),main="", xlab=names(tourist[i]))
#}
```

### 1.7d 用循环语句查看多个变量的分布：柱形图

```{r}
#for (i in 10:40){
 #  hist(tourist[,i],main="", xlab=names(tourist[i]))
  
   #带频数值标签的柱形图
 #  h<-hist(tourist[,i],plot=FALSE)

 # hist(tourist[,i],main="", xlab=names(tourist[i]),labels = TRUE,ylim=c(0,1.1*max(h$counts)))

  # 关于柱形图的标注问题，参考如下网址
  # https://stackoverflow.com/questions/9317948/how-to-label-histogram-bars-with-data-values-or-percents-in-r  
#}
```

## 1.8删除一个变量

```{r}
#产生一个新变量
tourist$sat2023<-tourist$sat1/5



```

删除使用NULL

```{r}
#删除这个刚生成的变量
tourist$sat2023<-NULL
```

## 1.9 save dataset

保存为rda格式数据，以后可以反复调用，注意从2.开始，不用把1.部分重复一遍。

```{r}

save(tourist,file=paste0(datafolder,"tourist_cn.rda"))
#save(tourist,file=paste(datafolder,"tourist.rda",sep=""))



```

```{r}
#输出为csv文件
#write.csv(tourist,paste(datafolder,'tourist2023cn.csv',sep=""),fileEncoding = "GBK")
```

# 2. basic analysis基础分析

```{r}

#设定工作文件夹, 根据自己电脑的情况进行修改
datafolder<-'D:/tdata/tourist_cn_R/'


```

load是用于载入保存了的rda格式数据

```{r}
load(paste0(datafolder,'tourist_cn.rda'))
```

## 2.1 frequency analysis频数分析

频数分析 Freq_gender Freq_thotel_gender

FreqGender FreqThotelGender

```{r}
#性别的频数表
Freq_gender<-table(tourist$gender)
Freq_gender

#性别的百分比
prop.table(Freq_gender)

#保留小数点后2位数字
round(prop.table(Freq_gender)*100,2)



```

为了输出更整洁一点，可以转换为数据框data frame

```{r}
#性别的频数表


#频数表gender_count为数据框data frame
gender_count<-as.data.frame(table(tourist$gender))


#频数表gender_percent为数据框data frame
gender_percent<-as.data.frame(round(prop.table(table(tourist$gender)),4))

#合并gender_count和gender_percent
gender_Freq<-cbind(gender_count,gender_percent$Freq)

#更改变量名
colnames(gender_Freq)<-c('性别','频数','Percent')

print(gender_Freq)
#输出为csv文件
#write.csv(gender_Freq,paste(datafolder,'gender_Freq.csv',sep=""),fileEncoding = "GBK")


```

如果用xlconnect，可以写入到特定excel格式文件

```{r}
#住宿类型的频数表


#频数表thotel_count为数据框data frame
thotel_count<-as.data.frame(table(tourist$thotel))


#频数表thotel_percent为数据框data frame
thotel_percent<-as.data.frame(round(prop.table(table(tourist$thotel)),4))

#合并thotel_count和thotel_percent
thotel_Freq<-cbind(thotel_count,thotel_percent$Freq)

#更改变量名
colnames(thotel_Freq)<-c('thotel','Freq','Percent')

print(thotel_Freq)
#输出为csv文件
#write.csv(thotel_Freq,paste(datafolder,'thotel_Freq.csv',sep=""),fileEncoding = "GBK")


```

### 2.1a 练习6 :频数分析

```{r}
##### 练习6 ： 对region变量，进行频数分析。
#区域的频数表



```



```{r}
##### 练习6答案 ： 对region变量，进行频数分析。

#频数表region_count为数据框data frame
region_count<-as.data.frame(table(tourist$region))


#百分表region_percent为数据框data frame
region_percent<-as.data.frame(round(prop.table(table(tourist$region)),4))

#合并region_count和region_percent
region_Freq<-cbind(region_count,region_percent$Freq)

#更改变量名
colnames(region_Freq)<-c('区域','频数','百分比')

print(region_Freq)
#输出为csv文件
#write.csv(region_Freq,paste(datafolder,'region_Freq.csv',sep=""),fileEncoding = "GBK")

```

### 2.1b 连续变量的柱形图

```{r}
hist(tourist$byear)
  # 带频数值标签的柱形图
h2<-hist(tourist$byear,plot=FALSE)
hist(tourist$byear,main="", xlab='出生年份')
hist(tourist$byear,main="", xlab='出生年份',labels = TRUE)
hist(tourist$byear,main="", xlab='出生年份',labels = TRUE,ylim=c(0,50))
hist(tourist$byear,main="", xlab='出生年份',labels = TRUE,ylim=c(0,50),col='lightblue')
hist(tourist$byear,main="", xlab='出生年份',labels = TRUE,ylim=c(0,50),col='cyan4')
#html color codes
#https://html-color.codes/
# r colors names
# https://r-charts.com/colors/


```

### 2.1c 练习7 age4的频数分析

```{r}

#### 综合一点的练习

##### 练习7 ： 
##### step1 参考1.4 根据byear，计算生成年龄变量age
##### step2 参考1.7 查看年龄变量age的分布
##### step3 参考1.5 把连续变量age，转换为年龄段变量，例如age4（20岁以下，20-40，40-60，60岁以上）
##### step4 参考2.1 对age4变量进行频数分析


```

练习7的答案

```{r}
# 练习7 答案
age<-2023-tourist$byear
hist(age)

attach(tourist)
tourist$age4[age  < 20] <- 1
tourist$age4[age >= 20 & age < 40] <- 2
tourist$age4[age >= 40 & age < 60] <- 3
tourist$age4[age >= 60] <- 4
detach(tourist)

tourist$age4 <- factor(tourist$age4,levels = c(1,2,3,4),labels = c("20岁以下", "20-40", "40-60","60岁以上"))

t3<-table(tourist$age4)
t3
round(prop.table(t3)*100,2)
```

## 2.2 crosstable and chi-square test

列联表与卡方检验

```{r}
### crosstable 
### 推荐的命名是这样的ct_thotel_gender，但是可以简化使用ct1
### ct_thotel_gender的好处是：（1）了解列联表涉及的变量；（2）方便查找替换

ct1<-table(tourist$thotel,tourist$gender)
prop.table(ct1)
prop.table(ct1,1)
prop.table(ct1,2)
round(prop.table(ct1,1)*100,2)
# summary 输出卡方检验结果chi-square test
summary(ct1)
```

为了输出更整洁一点，可以转换为数据框data frame

```{r}
#
#频数表ct1为数据框data frame
### 推荐的命名是这样的ct_thotel_gender，但是可以简化使用ct_thotel_gender
### ct_thotel_gender的好处是：（1）了解列联表涉及的变量；（2）方便查找替换

Freq_thotel_gender<-as.data.frame(table(tourist$thotel,tourist$gender))


#频数表ct_thotel_gender为数据框data frame
Percent_thotel_gender<-as.data.frame(round(prop.table(table(tourist$thotel,tourist$gender),1)*100,2))


#合并Freq_thotel_gender和Percent_thotel_gender
ct_thotel_gender<-cbind(Freq_thotel_gender,Percent_thotel_gender$Freq)

#更改变量名
colnames(ct_thotel_gender)<-c('住宿类型','性别','频数','百分比')

print(ct_thotel_gender)

#输出为csv文件
#write.csv(ct_thotel_gender,paste(datafolder,'tourist_hotel.csv',sep=""))

```

### 练习8 卡方检验

```{r}

##### 练习8 ： 利用2.2的卡方检验，检验如下研究假设
##### H1 收入income3与酒店类型thotel有关



```

## 2.3 independent T-test

```{r}
t.test(tourist$sat1~tourist$gender)

```

更整洁的表格，需要broom包

```{r}
ttest_sat1_gender<-t.test(tourist$sat1~tourist$gender)
tidy(ttest_sat1_gender)
```

6个满意度与性别关系的T检验，合并一起。

```{r}
# 计算6个t-test
# 第10-15列是满意度变量sat1-sat6，所以tourist[,10:15]与tourist[,c('sat1','sat2','sat3','sat4','sat5','sat6')]效果相同，直接引用变量名，程序容易理解，但是写起来麻烦

ttest_sat <- purrr::map(tourist[,10:15], ~t.test(.x ~ tourist$gender))

#整洁输出
report_ttest_sat<-data.frame(tidy(ttest_sat[[1]]))[1,]

for(i in 2:6){
  report_ttest_sat<-rbind(report_ttest_sat,data.frame(tidy(ttest_sat[[i]]))[1,])
}

# 增加变量名

for(i in 1:6){
  report_ttest_sat[i,1]<-names(tourist[9+i])
}

colnames(report_ttest_sat)<-c('变量','男','女','t','p.value','parameter','conf.low','conf.high','method','alternative')
report_ttest_sat
					


# 输出分析结果
write.csv(report_ttest_sat,paste(datafolder,'report_ttest_sat.csv',sep=""),fileEncoding = "GBK")

```

### 练习9 卡方检验

```{r}

##### 练习9 ： 利用2.3的T检验，检验如下研究假设
##### H3 女士比男士（gender）对景区的总体满意度（sat）更高



```

## 2.4 ANOVA

单因素方差分析

```{r}
aov_income3<-aov(sat1~income3,tourist)
summary(aov_income3)
tidy(aov_income3)
```

## 2.5 correlation

相关分析

```{r}
#两个变量关系的散点图
plot(tourist$income,tourist$te)
```

```{r}

cor.test(tourist$income,tourist$expense)

cor(tourist$income,tourist$expense,method="spearman")

```

更加整洁的输出

```{r}



c1<-cor.test(tourist$income,tourist$expense)
tidy(c1)
```

相关系数矩阵

```{r}
# 第5列是income，第6列是expense，第10-15列是满意度变量sat1-sat6
cor_data<-tourist[,c(5,6,10:15)]
mydata.cor = cor(cor_data)
mydata.cor
```

rcorr来自Hmisc包

```{r}

mydata.rcorr <- rcorr(as.matrix(cor_data))
mydata.coeff <- mydata.rcorr$r
mydata.coeff 
mydata.p <- mydata.rcorr$P
mydata.p
```

tidy report
```{r}
#rcorr来自Hmisc包
tidy(mydata.rcorr)
```

相关系数矩阵图

```{r}

corrplot(mydata.cor)
```

# 3.ggplot2绘图
设定工作目录

```{r}
datafolder<-'D:/tdata/tourist_cn_R/'
```

加载数据
```{r}
load(paste0(datafolder,'tourist_cn.rda'))
```

## 3.1 qplot快捷做图

```{r}
qplot(tourist$income,tourist$expense)
```
## 3.2 标准的ggplot2 散点图

### 3.2a 散点图的基本框架geom_point
```{r}
ggplot(tourist,aes(x=income,y=expense))+
  geom_point(size=1)
#size是散点图总点的大小
```
### 3.2b 散点图的分组着色

基于gender给数据点着色
```{r}
ggplot(tourist,aes(x=income,y=expense))+
  geom_point(size=2,aes(colour=gender))
```

type3给数据点着色

```{r}
ggplot(tourist,aes(x=income,y=expense))+
  geom_point(size=2,aes(colour=type3))
```
### 3.2c 散点图的分面（facet）

使用type3进行分面
```{r}
ggplot(tourist,aes(x=income,y=expense))+
  geom_point(size=2)+
  facet_grid(type3 ~.)
```
使用gender进行分面
```{r}
ggplot(tourist,aes(x=income,y=expense))+
  geom_point(size=2,aes(colour=type3))+
  facet_grid(gender ~.)
```

使用type3,type2进行分面

### 3.2d 保存输出图片

```{r}
#png 用于保存图片
png(paste0(datafolder,'scatexpenser_income_expense.png'))
scatexpenser_income_expense<-ggplot(tourist,aes(x=income,y=expense))+
  geom_point(size=2,aes(colour=type3))+
  facet_grid(type3 ~type2)


scatexpenser_income_expense

#print 用于保存图片
print(scatexpenser_income_expense)
dev.off()

```

## 3.3 柱形图（histogram）

### 3.3a 柱形图的基本框架（geom_histogram）

```{r}
ggplot(tourist,aes(x=sat2,fill=gender))+
  geom_histogram(binwidth=1,colour="white")
  
```


### 3.3b 使用type3、type2进行分面
```{r}
ggplot(tourist,aes(x=sat2,fill=type3))+
  geom_histogram(binwidth=1,colour="white")+
  facet_grid(type2 ~type3)
```

# 4. map地图

```{r}
#设定工作文件夹, 根据自己电脑的情况进行修改
datafolder<-'D:/tdata/tourist_cn_R/'

```

load是用于载入保存了的rda格式数据

```{r}
load(paste0(datafolder,'tourist_cn.rda'))
```

leaflet地图
```{r}
# 中心位置，缩放等级
center_lon = median(tourist$longitude,na.rm = TRUE)
center_lat = median(tourist$latitude,na.rm = TRUE)

leaflet(tourist) %>% addProviderTiles("Esri.NatGeoWorldMap") %>%
  addCircles(lng = ~longitude, lat = ~latitude,radius = ~sqrt(income/10000000))  %>%
  setView(lng=center_lon, lat=center_lat,zoom = 10)
```

# 5 IPA图


IPA（importance-performance analysis）重要性-表现分析方法
IPA相对简单，但是在管理、营销领域广泛使用

Mejia, C., Bąk, M., Zientara, P., & Orlowski, M. (2022).
Importance-Performance Analysis of Socially Sustainable Practices in
U.s. Restaurants: A Consumer Perspective in the Quasi-Post-Pandemic
Context. International Journal of Hospitality Management, 103, 103209.
https://doi.org/10.1016/j.ijhm.2022.103209

```{r}

#设定工作文件夹, 根据自己电脑的情况进行修改
datafolder<-'D:/tdata/tourist_cn_R/'


```

load是用于载入保存了的rda格式数据

```{r}
load(paste0(datafolder,'tourist_cn.rda'))
```


## 5.1 data wrangling

```{r}
review_label<-as.data.frame(c('评论数量',	'评论发表日期',	'评论与酒店性能关联',	'正向评价',	'发布者资信度'))
colnames(review_label)<-c("label")
tourist1<-tourist[,c("ri1","ri2","ri3","ri4","ri5","rp1","rp2","rp3","rp4","rp5")]
hdata<-colMeans(tourist1)


im<-tourist[,c("ri1","ri2","ri3","ri4","ri5")]
pe<-tourist[,c("rp1","rp2","rp3","rp4","rp5")]


imdata<-as.data.frame(colMeans(im))
colnames(imdata)<-c("importance")

pedata<-as.data.frame(colMeans(pe))
colnames(pedata)<-c("performance")
ipa<-cbind(imdata,pedata,review_label)

```

## 5.2 计算x轴、y轴的取值范围

```{r}
print("------------importance----------")
print(paste("min :",round(min(ipa[1]),2)))

print(paste("max :",round(max(ipa[1]),2)))
print("------------performance----------")
print(paste("min :",round(min(ipa[2]),2)))
print(paste("max :",round(max(ipa[2]),2)))
print(colMeans(ipa[,1:2]))
```

## 5.3 IPA图

```{r}
# 基本图形架构
# xlim是x轴的区间——参考了5.2的计算结果
# ylim是x轴的区间——参考了5.2的计算结果
# labs是坐标轴的标注
# geom_vline是增加一条垂直的线，取值xintercept，颜色color，线型linetype
# geom_hline是增加一条水平的线，取值yintercept，颜色color，线型linetype
# geom_abline是增加一条斜线：y=ax+b，截距intercept，斜率slope，颜色color，粗细size
# geom_text_repel，数据点的标注问题，不重叠

ipa_chart<-ggplot(ipa,aes(x=performance,y=importance,))+
  geom_point(size=2)+ 
  xlim(4.2, 4.8)+
  ylim(4.2, 4.8)+
  labs(x = "表现",y="重要性")+
  geom_vline(xintercept=4.516,color = "red",linetype="dashed")+
  geom_hline(yintercept=4.555,color = "red",linetype="dashed")+
  geom_abline(intercept = 0, slope = 1, color="blue",size=0.5)+
  geom_text_repel(aes(label=label),size=4,hjust=0.3,vjust=0.3)

ipa_chart
#图片输出
#ggsave(paste(datafolder,'ipa_chart.png',sep=""),ipa_chart, width = 10)

```

# 6.EFA探索性因子分析


```{r}
#设定工作文件夹, 根据自己电脑的情况进行修改
datafolder<-'D:/tdata/tourist_cn_R/'

```

load是用于载入保存了的rda格式数据

```{r}
load(paste0(datafolder,'tourist_cn.rda'))
```


## 6.1 datasets数据集整理

选取变量的2种方法：
（1）用列号
zhdata<-tourist[,c(34:40)]

（2）用变量名称
zhdata<-tourist[,c('zh1','zh2','zh3','zh4','zh5','zh6','zh7')]


```{r}
#satdata<-tourist[,c(10:15)]
#tedata<-tourist[,c(26:33)]
zhdata<-tourist[,c(34:40)]
#zhdata<-tourist[,c('zh1','zh2','zh3','zh4','zh5','zh6','zh7')]


#head(satdata)
#head(tedata)
head(zhdata)
```

## 6.2 kmo test

KMO(Kaiser-Meyer-Olkin)
检验值0-1之间，值越大越适合探索性因子分析EFA，通常至少0.7以上

```{r}
kmo_test<-KMO(zhdata)
print(kmo_test, stats = c("both", "MSA", "KMO"), vars = "all", sort = FALSE, show = "all", digits = 3)
```

## 6.3 Bartlett's test

Bartlett球形检验 Bartlett's test of sphericity
p<0.05，统计上显著适合探索性因子分析EFA

```{r}
cortest.bartlett(zhdata)
```

## 6.4判断因子数量

```{r}
parallel <- fa.parallel(zhdata, fm = 'minres', fa = 'both',n.iter=100)

```

推荐提取2个因子

## 6.5 提取2个因子

```{r}
efa2 <- fa(zhdata,nfactors = 2,rotate = "oblimin",fm="minres")
print(efa2)

print(efa2$loadings,cutoff = 0.5)
```

## 6.6图示

EFA的图示1
```{r}
factor.plot(efa2,labels=rownames(efa2$loadings))
```

EFA的图示2
```{r}
fa.diagram(efa2)

```

# 7.CFA验证性因子分析Confirmatory Factor Analysis


```{r}
#设定工作文件夹, 根据自己电脑的情况进行修改
datafolder<-'D:/tdata/tourist_cn_R/'

```

load是用于载入保存了的rda格式数据

```{r}
load(paste0(datafolder,'tourist_cn.rda'))
```


     MR1        MR2
zh1 0.821
zh2 0.911
zh3 0.875
zh4 0.733
zh5            0.737
zh6            0.868
zh7            0.910

## 7.1 CFA模型设定

```{r}
cfa_model <- " D1 =~ zh1 + zh2 + zh3 + zh4 
 D2 =~ zh5 + zh6 + zh7 "
```

## 7.2 CFA模型估计

```{r}
cfa.fit <- cfa(cfa_model, data=zhdata)
summary(cfa.fit, fit.measures=TRUE)

```

## 7.3 CFA模型图示

```{r}
semPaths(cfa.fit, what="est", fade=FALSE, residuals=FALSE,edge.label.cex=0.75)
```

# 8.信度分析reliability analysis
信度分析参考资料

https://lhbikos.github.io/ReC_Psychometrics/rxx.html


```{r}
#设定工作文件夹, 根据自己电脑的情况进行修改
datafolder<-'D:/tdata/tourist_cn_R/'

```

load是用于载入保存了的rda格式数据

```{r}
load(paste0(datafolder,'tourist_cn.rda'))
```


```{r}
if(!require(psych)){install.packages("psych")}

```

维度1的信度分析

```{r}
zh_d1<-tourist[,c('zh1','zh2','zh3','zh4')]
#alpha(zh_d1)
splitHalf(zh_d1, raw = TRUE, brute = TRUE)

```




维度2的信度分析

```{r}
zh_d2<-tourist[,c('zh5','zh6','zh7')]
splitHalf(zh_d2, raw = TRUE, brute = TRUE)
```
