##讀取檔案、提取資料與製造變項
#這是一般 TXT 檔，檔頭有變項名稱
#資料來自於 NHIS 2010 調查,取 1955 年以前出生者（55歲以上）
#視力困難、聽力困難、移動困難與溝通困難變項，分數越高越有困難
dta <- read.table(file = 'data/Quality_of_Life.txt', header = TRUE)

#dta 的類型是資料框架（data frame）
class(dta)
## [1] "data.frame"
#看看dta維度
dim(dta)
## [1] 2077    7
#利用names看看變項名稱
names(dta)
## [1] "性別" "教育" "年齡" "視力" "聽力" "移動" "溝通"
#看前六筆
#程式報表1.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] 女
## Levels: 女 男
#看看第九列，這是資料框架
#看看第一欄，這「不是」資料框架，是類別
dta[9, ]
##   性別 教育     年齡 視力 聽力 移動 溝通
## 9   女 小學 56到60歲 1.33    1  1.2    1

dta[,1]

前者是類別變項，後者是資料框架

dta[,‘教育’] dta[‘教育’]

#還是資料框架
dta[5:7, c('視力', '聽力')]
##   視力 聽力
## 5 1.67 2.33
## 6 1.00 1.00
## 7 1.00 1.33
#檢視資料結構
#程式報表1.2
str(dta)
## 'data.frame':    2077 obs. of  7 variables:
##  $ 性別: Factor w/ 2 levels "女","男": 1 1 2 1 2 1 2 1 1 1 ...
##  $ 教育: Factor w/ 5 levels "大學","小學",..: 5 1 1 5 5 1 5 5 2 5 ...
##  $ 年齡: Factor w/ 6 levels "56到60歲","61到65歲",..: 1 3 5 6 4 1 5 6 1 1 ...
##  $ 視力: 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 ...
#看看類別變項的屬性，這是五個水準類別變項
attributes(dta[, '教育'])
## $levels
## [1] "大學" "小學" "初中" "高中" "專科"
## 
## $class
## [1] "factor"
#利用table指令看看類別變項分佈
#程式報表1.3
table(dta[, '教育'])
## 
## 大學 小學 初中 高中 專科 
##  466  107  163  188 1153
#將類別變項轉成數值變項，再利用table
#程式報表1.3
table(as.numeric(dta[, '教育']))
## 
##    1    2    3    4    5 
##  466  107  163  188 1153
#指令summary是常用指令，可以用在不同物件上
#指令summary作用結果視作用物件而定
#這是用在資料框架上
#程式報表1.4
summary(dta['教育'])
##    教育     
##  大學: 466  
##  小學: 107  
##  初中: 163  
##  高中: 188  
##  專科:1153
#這是用在類別變項上
#程式報表1.4
summary(dta[, '教育'])
## 大學 小學 初中 高中 專科 
##  466  107  163  188 1153
#第一個是數值，第二個是資料框架
#指令t用來轉置矩陣，資料框架被轉成矩陣了
class(dta[, '視力'])
## [1] "numeric"
class(dta['視力'])
## [1] "data.frame"
class(t(dta['視力']))
## [1] "matrix"
#製造平均數兩種寫法
#很直覺，不過如果有遺漏值就麻煩了
dta$功能 <-(dta$視力 + dta$聽力 + dta$移動 + dta$溝通)/4
                          
#有遺漏值會警告我們
tail(dta$功能 <- rowMeans(dta[, 4:7]))
## [1] 1.0000 1.1450 1.4425 1.2500 1.0000 1.0950
###基本統計量（示範應用函數）
#看看資料基本統計
#程式報表1.5
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
#計算健康情形標準差
#利用 with 會讓指令比較好寫好讀
sd(dta$功能)
## [1] 0.3302689
with(dta, sd(功能))
## [1] 0.3302689
#計算偏態
#先安裝再載入套件 moments
#注意，安裝只需要一次，但載入則是每次
install.packages('moments', dependencies = TRUE)
require(moments)
skewness(dta$功能)
#一次計算所有平均數
sapply(dta[, -c(1:3)], mean)
##     視力     聽力     移動     溝通     功能 
## 1.208199 1.301815 1.515315 1.077756 1.275772
#定義一個函數，可以同時算出平均數、標準差、偏態與峰度
#一次計算多個變項平均數、標準差、偏態與峰度
#程式報表1.6
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)
## Loading required package: moments
## Warning: package 'moments' was built under R version 3.2.2
##            視力      聽力      移動       溝通      功能
## [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
#把不同性別、教育者平均數算出來（兩種寫法）
#第二種寫法把輸出放到 hlth 物件了，沒有輸出到螢幕
#程式報表1.7
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 = '功能'))

#依移動排序，可以看到男性移動困難退化的較慢
#程式報表1.8
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('小學', '初中', '高中', '專科', '大學'))

#不同教育水準健康功能盒鬚圖
#把資料改成水平，更容易比較
#這是直接針對原始資料，而非統計量畫圖
#圖1.2
plot(功能 ~ 教育, data = dta, horizontal = T)

#加上格線，有助比較
grid()

#載入 lattice，準備畫圖
#取出四個要畫圖的變項
require(lattice)
## Loading required package: lattice
## Warning: package 'lattice' was built under R version 3.2.2
v4 <- dta[, 4:7]

#把變項年齡與性別上的直方圖記下來
p <- lapply(names(v4), function(x) histogram( ~ v4[x] | 年齡 + 性別, data = dta,
            xlab = x )) 
                            
#利用 gridExtra 排一下四張圖
#圖1.3
library(gridExtra)
## Warning: package 'gridExtra' was built under R version 3.2.2
grid.arrange(p[[1]], p[[2]], p[[3]], p[[4]], nrow = 2, ncol = 2)

#載進 ggplot2，準備畫圖
#Hmisc提供 ggplot2 額外功能
require(ggplot2)
## Loading required package: ggplot2
## Warning: package 'ggplot2' was built under R version 3.2.2
require(Hmisc)
## Loading required package: Hmisc
## Warning: package 'Hmisc' was built under R version 3.2.2
## Loading required package: grid
## Loading required package: survival
## Warning: package 'survival' was built under R version 3.2.2
## Loading required package: Formula
## Warning: package 'Formula' was built under R version 3.2.2
## 
## Attaching package: 'Hmisc'
## 
## The following objects are masked from 'package:base':
## 
##     format.pval, round.POSIXt, trunc.POSIXt, units
#調整圖形，避免重疊
pd <- position_dodge(width = .3)

#繪製不同年齡生理健康的平均數
#加上利用拔靴法得到的平均數信賴區間
#圖1.4
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))

#載入reshape2、把資料由寬形變長形
require(reshape2)
## Loading required package: reshape2
## Warning: package 'reshape2' was built under R version 3.2.2
dtal <- melt(dta[, -8], variable.name = '功能項目', value.name = '分數')
## Using 性別, 教育, 年齡 as id variables
#看看資料
#程式報表1.9
head(dtal)
##   性別 教育     年齡 功能項目 分數
## 1   女 專科 56到60歲     視力 1.00
## 2   女 大學 66到70歲     視力 1.00
## 3   男 大學 76到80歲     視力 1.00
## 4   女 專科 81歲以上     視力 2.67
## 5   男 專科 71到75歲     視力 1.67
## 6   女 大學 56到60歲     視力 1.00
#不同年齡性別各功能項目平均數
#圖1.5
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))
## ymax not defined: adjusting position using y instead
## ymax not defined: adjusting position using y instead
## ymax not defined: adjusting position using y instead
## ymax not defined: adjusting position using y instead

#不同性別在功能上平均數檢定（兩種寫法）
#程式報表1.10
t.test(功能 ~ 性別, data = dta)
## 
##  Welch Two Sample t-test
## 
## data:  功能 by 性別
## t = 2.6445, df = 1920.4, p-value = 0.008249
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  0.009959744 0.067135903
## sample estimates:
## mean in group 女 mean in group 男 
##         1.292030         1.253482
with(dta, t.test(功能 ~ 性別))
## 
##  Welch Two Sample t-test
## 
## data:  功能 by 性別
## t = 2.6445, df = 1920.4, p-value = 0.008249
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  0.009959744 0.067135903
## sample estimates:
## mean in group 女 mean in group 男 
##         1.292030         1.253482