HW1
fce
2022-03-01
fLoc <- "https://cpcheng.neocities.org/Rbook/01/data/Quality_of_Life.txt"dta <- read.table(file = fLoc, header = TRUE)class(dta)## [1] "data.frame"看看dta維度： 多少行、多少欄(變項)
dim(dta)## [1] 2077 7#利用names看看變項名稱
names(dta)## [1] "性別" "教育" "年齡" "視力" "聽力" "移動" "溝通"看前六筆
head(dta)## 性別 教育 年齡 視力 聽力 移動 溝通
## 1 女 專科 56到60歲 1.00 1.33 1.2 1.0
## 2 女 大學 66到70歲 1.00 1.33 1.2 1.5
## 3 男 大學 76到80歲 1.00 1.33 1.0 1.0
## 4 女 專科 81歲以上 2.67 2.00 1.0 1.0
## 5 男 專科 71到75歲 1.67 2.33 1.0 1.0
## 6 女 大學 56到60歲 1.00 1.00 1.0 1.0#看看第一列第一欄對應資料
dta[1, 1]## [1] "女"這是類別資料
#看看第九列，這是資料框架
class(dta[9, ])## [1] "data.frame"#看看第一欄，這「不是」資料框架，是類別資料
class(dta[,1])## [1] "character"#前者是類別變項，後者是資料框架
class(dta[,'教育'])## [1] "character"class(dta['教育'])## [1] "data.frame"#還是資料框架
class(dta[5:7, c('視力', '聽力')])## [1] "data.frame"檢視資料結構
str(dta)## 'data.frame': 2077 obs. of 7 variables:
## $ 性別: chr "女" "女" "男" "女" ...
## $ 教育: chr "專科" "大學" "大學" "專科" ...
## $ 年齡: chr "56到60歲" "66到70歲" "76到80歲" "81歲以上" ...
## $ 視力: num 1 1 1 2.67 1.67 1 1 3.67 1.33 1 ...
## $ 聽力: num 1.33 1.33 1.33 2 2.33 1 1.33 1.33 1 1 ...
## $ 移動: num 1.2 1.2 1 1 1 1 1 2.33 1.2 1 ...
## $ 溝通: num 1 1.5 1 1 1 1 1 1 1 1 ...看看類別變項的屬性，有五個level的categorical variable
attributes(dta[, "教育"])## NULL將變項的attribute改變
dta[,1:3] <- lapply(dta[, 1:3], factor)用table指令檢視類別變項分佈
table(dta[, "教育"])##
## 大學 小學 初中 高中 專科
## 466 107 163 188 1153將類別變項轉成數值變項，再利用table檢視之
table(as.numeric(dta[, '教育']))##
## 1 2 3 4 5
## 466 107 163 188 1153當summary用在data frame
summary(dta['教育'])## 教育
## 大學: 466
## 小學: 107
## 初中: 163
## 高中: 188
## 專科:1153當summary出現在類別變項前
summary(dta[, '教育'])## 大學 小學 初中 高中 專科
## 466 107 163 188 1153class指令
#這一個
class(dta[, '視力'])## [1] "numeric"#另一個
class(dta['視力'])## [1] "data.frame"#當class加上指令t(轉置矩陣)，於是資料框架被轉成矩陣了
class(t(dta['視力']))## [1] "matrix" "array"平均數的兩種寫法
#第一個很直接，只是如果有遺漏值就麻煩了
dta$功能 <- (dta$視力 + dta$聽力 + dta$移動 + dta$溝通)/4;
head(dta$功能)## [1] 1.1325 1.2575 1.0825 1.6675 1.5000 1.0000#第二個則可以警告是否有遺漏值
tail(dta$功能 <- rowMeans(dta[, 4:7]))## [1] 1.0000 1.1450 1.4425 1.2500 1.0000 1.0950基本統計量／描述統計 Descriptive statistic
#summary指令又出現了
summary(dta)## 性別 教育 年齡 視力 聽力
## 女:1201 大學: 466 56到60歲:393 Min. :1.000 Min. :1.000
## 男: 876 小學: 107 61到65歲:456 1st Qu.:1.000 1st Qu.:1.000
## 初中: 163 66到70歲:383 Median :1.000 Median :1.000
## 高中: 188 71到75歲:283 Mean :1.208 Mean :1.302
## 專科:1153 76到80歲:217 3rd Qu.:1.330 3rd Qu.:1.330
## 81歲以上:345 Max. :4.000 Max. :4.000
## 移動 溝通 功能
## Min. :1.000 Min. :1.000 Min. :1.000
## 1st Qu.:1.000 1st Qu.:1.000 1st Qu.:1.000
## Median :1.000 Median :1.000 Median :1.165
## Mean :1.515 Mean :1.078 Mean :1.276
## 3rd Qu.:1.800 3rd Qu.:1.000 3rd Qu.:1.417
## Max. :4.000 Max. :4.000 Max. :4.000標準差的計算
sd(dta$功能)## [1] 0.3302689#再來認識一個新的指令－with
with(dta, sd(功能))## [1] 0.3302689安裝與載入
install.packages("moments", repos='https://cran.rstudio.com/')## 將程式套件安載入 'C:/Users/user/Documents/R/win-library/4.1'
## (因為 'lib' 沒有被指定)## package 'moments' successfully unpacked and MD5 sums checked
##
## The downloaded binary packages are in
## C:\Users\user\AppData\Local\Temp\RtmpyQg33j\downloaded_packageslibrary(moments)## Warning: 套件 'moments' 是用 R 版本 4.1.1 來建造的計算偏態
with(dta, skewness(功能))## [1] 1.821716一次計算所有變項的平均數
sapply(dta[, -c(1:3)], mean)## 視力 聽力 移動 溝通 功能
## 1.208199 1.301815 1.515315 1.077756 1.275772Conditional means
aggregate(cbind(視力, 聽力, 移動, 溝通) ~ 教育, data = dta, mean)## 教育 視力 聽力 移動 溝通
## 1 大學 1.141953 1.243047 1.262339 1.036481
## 2 小學 1.401402 1.294299 1.913178 1.172897
## 3 初中 1.355153 1.480061 1.890982 1.180982
## 4 高中 1.242660 1.423351 1.704096 1.103723
## 5 專科 1.190650 1.281249 1.496748 1.066782定義一個函數 可以同時算出平均數、標準差、偏態與峰度 一次計算多個變項平均數、標準差、偏態與峰度
my_summary <- function(x) {
require(moments)
funs <- c(mean, sd, skewness, kurtosis)
sapply(funs, function(f) f(x, na.rm = TRUE))
}
sapply(dta[, c(4:8)], my_summary)## 視力 聽力 移動 溝通 功能
## [1,] 1.2081993 1.3018151 1.5153154 1.0777564 1.2757715
## [2,] 0.4254268 0.4603758 0.7489312 0.2909127 0.3302689
## [3,] 3.0324554 2.1140343 1.4311597 4.9739873 1.8217162
## [4,] 14.4583674 8.6589573 3.9949684 34.3036465 7.7898425計算不同教育程度健康功能標準差（兩種寫法）
#第一種
tapply(dta$功能, dta$年齡, sd)## 56到60歲 61到65歲 66到70歲 71到75歲 76到80歲 81歲以上
## 0.2915871 0.2925204 0.2703624 0.3217985 0.3057586 0.4122367#第二種
with(dta, tapply(功能, 年齡, sd))## 56到60歲 61到65歲 66到70歲 71到75歲 76到80歲 81歲以上
## 0.2915871 0.2925204 0.2703624 0.3217985 0.3057586 0.4122367計算不同性別、教育的平均數（兩種寫法）
#第一種
with(dta, aggregate(dta[, 4:8], by = list(性別, 年齡), FUN = mean))## Group.1 Group.2 視力 聽力 移動 溝通 功能
## 1 女 56到60歲 1.201143 1.164667 1.411714 1.052381 1.207476
## 2 男 56到60歲 1.156393 1.243607 1.244754 1.054645 1.174850
## 3 女 61到65歲 1.197425 1.204664 1.466567 1.072761 1.235354
## 4 男 61到65歲 1.169894 1.260160 1.286915 1.098404 1.203843
## 5 女 66到70歲 1.167885 1.210192 1.569375 1.050481 1.249483
## 6 男 66到70歲 1.147314 1.370914 1.329371 1.065714 1.228329
## 7 女 71到75歲 1.238194 1.236065 1.604387 1.077419 1.289016
## 8 男 71到75歲 1.184531 1.352422 1.461328 1.093750 1.273008
## 9 女 76到80歲 1.236484 1.346016 1.606484 1.050781 1.309941
## 10 男 76到80歲 1.142022 1.263933 1.369663 1.134831 1.227612
## 11 女 81歲以上 1.330259 1.480043 1.947802 1.099138 1.464310
## 12 男 81歲以上 1.318142 1.659027 1.892566 1.132743 1.500619#第二種
hlth <- aggregate(cbind(視力, 聽力, 移動, 溝通, 功能) ~ 性別 + 年齡,
data = dta, mean)#命名 hlth 的前兩欄
names(hlth)[1:2] <- c('性別', '年齡')不同性別與年齡的健康功能圖
#圖1.1
with(hlth, dotchart(t(功能), group = 性別, labels = 年齡, xlab = '功能'))
依移動排序，可以看到男性移動困難退化的較慢
hlth[order(hlth$移動), c(1:2, 5)]## 性別 年齡 移動
## 2 男 56到60歲 1.244754
## 4 男 61到65歲 1.286915
## 6 男 66到70歲 1.329371
## 10 男 76到80歲 1.369663
## 1 女 56到60歲 1.411714
## 8 男 71到75歲 1.461328
## 3 女 61到65歲 1.466567
## 5 女 66到70歲 1.569375
## 7 女 71到75歲 1.604387
## 9 女 76到80歲 1.606484
## 12 男 81歲以上 1.892566
## 11 女 81歲以上 1.947802繪圖顯示
#先將教育程度水準重排
dta$教育 <- factor(dta$教育, levels = c('小學', '初中', '高中', '專科', '大學'))先繪出不同教育水準的健康功能盒鬚圖 把資料顯示改成水平，較容易做比較 這是原始資料，而非統計量的圖像
plot(功能 ~ 教育, data = dta, horizontal = T)
#加上格線，有助比較
grid()
library(lattice)#取出四個要畫圖的變項
v4 <- dta[, 4:7]p <- lapply(names(v4), function(x) histogram( ~ v4[x] | 年齡 + 性別, data = dta,
xlab = x )) 利用gridExtra排一下四張圖
install.packages("gridExtra", repos='https://cran.rstudio.com/')## 將程式套件安載入 'C:/Users/user/Documents/R/win-library/4.1'
## (因為 'lib' 沒有被指定)## package 'gridExtra' successfully unpacked and MD5 sums checked
##
## The downloaded binary packages are in
## C:\Users\user\AppData\Local\Temp\RtmpyQg33j\downloaded_packageslibrary(gridExtra)## Warning: 套件 'gridExtra' 是用 R 版本 4.1.2 來建造的grid.arrange(p[[1]], p[[2]], p[[3]], p[[4]], nrow = 2, ncol = 2)
繪製不同年齡在生理健康的平均數 加上利用拔靴法得到的平均數信賴區間
#先檢查再安裝
#載進 ggplot2 in tidyverse，準備畫圖
#Hmisc提供額外功能
if (!require(pacman)) install.packages("pacman")## 載入需要的套件：pacman## Warning: 套件 'pacman' 是用 R 版本 4.1.2 來建造的# 載入套件
pacman::p_load(tidyverse, Hmisc, here)
#調整圖形，避免重疊
pd <- position_dodge(width = .3)
ggplot(data = dta, aes(x = 年齡, y = 功能, color = 性別) ) +
stat_summary(fun.data = 'mean_cl_boot', position = pd) +
facet_grid(教育 ~ . ) +
labs(x = '年齡', y = '功能分數') +
coord_cartesian(ylim = range(dta['功能'])) +
coord_flip() +
scale_color_manual(values = c('black', 'gray'),
guide = guide_legend(title = '性別', reverse = TRUE)) +
theme_bw() +
theme(legend.justification = c(1, 0), legend.position = c(1, 0))## Coordinate system already present. Adding new coordinate system, which will replace the existing one.
示範資料重新結構 把資料由寬形變長形
dtal <- gather(dta, key = '功能項目', value = '分數', 4:7)
head(dtal)## 性別 教育 年齡 功能 功能項目 分數
## 1 女 專科 56到60歲 1.1325 視力 1.00
## 2 女 大學 66到70歲 1.2575 視力 1.00
## 3 男 大學 76到80歲 1.0825 視力 1.00
## 4 女 專科 81歲以上 1.6675 視力 2.67
## 5 男 專科 71到75歲 1.5000 視力 1.67
## 6 女 大學 56到60歲 1.0000 視力 1.00p <- ggplot(data = dtal, aes(x = 年齡, y = 分數,
shape = reorder(功能項目, 分數, mean), group = reorder(功能項目, 分數))) +
stat_summary(fun.data = mean_se, geom = 'errorbar', position = pd) +
stat_summary(fun.y = mean, geom = 'point', position = pd, size = 3) +
stat_summary(fun.y = mean, geom = 'line', position = pd, linetype = 'dotted') +
facet_grid( . ~ 性別) +
labs(x = '年齡', y = '分數') +
scale_shape_manual(values = c(2, 8, 1, 6),
guide = guide_legend(title = '功能項目', reverse = TRUE)) +
theme_bw() +
theme(legend.justification = c(0, 1), legend.position = c(0, 1))## Warning: `fun.y` is deprecated. Use `fun` instead.
## `fun.y` is deprecated. Use `fun` instead.print(p)
存檔、圖
dir.create("W3")
dir.create(here("W3", "data"))
write.csv(dta, here("W3", "data", "quality_of_life.csv"), quote = FALSE, row.names = FALSE)
dir.create(here("W3", "figs"))
ggsave(plot = p, here("W3", "figs", "Fig1_5.png"), device = "png")## Saving 7 x 5 in image