1 In-class 1

Theme= simplex

範例一:美國婦女的身高和體重

1.1 套內數據

十五位三十到十九歲美國婦女的身高和體重.

women
   height weight
1      58    115
2      59    117
3      60    120
4      61    123
5      62    126
6      63    129
7      64    132
8      65    135
9      66    139
10     67    142
11     68    146
12     69    150
13     70    154
14     71    159
15     72    164

1.2 基本統計圖形

下面R程式碼畫出women數據集的散點圖.

plot(women, type = 'b', xlab = "Height (in)", ylab = "Weight (lb)")
grid()

1.3 線性模型分析

\[y_i = \beta_0 + \beta_1 x_i + \epsilon_i ,~~ \epsilon_i \sim N(0, \sigma^2) \]

體重 = 截距參數 + 斜率參數 x 身高 + 殘差(常態分佈)

1.3.1 分析概要報表

小數點4位,去掉星星.

options(digits = 4, show.signif.stars = FALSE)
summary(m0 <- lm(weight ~ height, data = women))

Call:
lm(formula = weight ~ height, data = women)

Residuals:
   Min     1Q Median     3Q    Max 
-1.733 -1.133 -0.383  0.742  3.117 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)
(Intercept) -87.5167     5.9369   -14.7  1.7e-09
height        3.4500     0.0911    37.9  1.1e-14

Residual standard error: 1.53 on 13 degrees of freedom
Multiple R-squared:  0.991, Adjusted R-squared:  0.99 
F-statistic: 1.43e+03 on 1 and 13 DF,  p-value: 1.09e-14

根據這份數據, 三十到十九歲美國的婦女們,平均身高多出一英寸時, 平均體重大約增加3.45英磅(誤差為0.09). 殘差估計為\(\hat{\sigma} = 1.53\).

1.3.2 方差分析表

anova(m0)
Analysis of Variance Table

Response: weight
          Df Sum Sq Mean Sq F value  Pr(>F)
height     1   3333    3333    1433 1.1e-14
Residuals 13     30       2

1.3.3 模型擬合圖

plot(women, xlab = "Height (in)", ylab = "Weight (lb)")
abline(m0, lty = 2)
grid()

1.3.4 殘差圖

檢查殘差分配有沒有規律

plot(resid(m0) ~ fitted(m0), xlab = "Fitted values", 
     ylab = "Residuals", ylim = c(-3.5, 3.5))
grid()
abline(h = 0, lty = 2)

1.3.5 驗證殘差常態分佈

qqnorm(resid(m0))
qqline(resid(m0))
grid()

1.3.6 結束

顯示演練單元信息

sessionInfo()
R version 4.1.1 (2021-08-10)
Platform: x86_64-w64-mingw32/x64 (64-bit)
Running under: Windows 10 x64 (build 19043)

Matrix products: default

locale:
[1] LC_COLLATE=Chinese (Traditional)_Taiwan.950 
[2] LC_CTYPE=Chinese (Traditional)_Taiwan.950   
[3] LC_MONETARY=Chinese (Traditional)_Taiwan.950
[4] LC_NUMERIC=C                                
[5] LC_TIME=Chinese (Traditional)_Taiwan.950    

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

loaded via a namespace (and not attached):
 [1] digest_0.6.27    R6_2.5.1         jsonlite_1.7.2   magrittr_2.0.1  
 [5] evaluate_0.14    highr_0.9        rlang_0.4.11     stringi_1.7.4   
 [9] jquerylib_0.1.4  bslib_0.3.0      rmarkdown_2.11.2 tools_4.1.1     
[13] stringr_1.4.0    xfun_0.26        yaml_2.2.1       fastmap_1.1.0   
[17] compiler_4.1.1   htmltools_0.5.2  knitr_1.34       sass_0.4.0

2 In-class 2

ch02 file

# 讀檔案，資料來自於 TIMSS 2011 年
dta <- read.table(file = "TIMSSmath.txt", header = TRUE, fileEncoding="big5")
# 檢視資料結構
# 程式報表2.2
str(dta)
'data.frame':   15549 obs. of  14 variables:
 $ 國家    : chr  "台灣" "台灣" "台灣" "台灣" ...
 $ 性別    : chr  "女生" "女生" "女生" "女生" ...
 $ 父母教育: chr  "高中" "高中" "國小以下" "初中" ...
 $ 數學    : num  729 776 718 607 658 ...
 $ 化學    : num  625 662 709 633 618 ...
 $ 地科    : num  608 666 629 620 614 ...
 $ 生物    : num  651 628 617 592 560 ...
 $ 物理    : num  618 617 663 616 591 ...
 $ 電腦    : chr  "有" "有" "有" "有" ...
 $ 書桌    : chr  "無" "有" "有" "無" ...
 $ 書      : chr  "有" "有" "無" "無" ...
 $ 房間    : chr  "無" "無" "無" "無" ...
 $ 網路    : chr  "有" "有" "有" "有" ...
 $ 書量    : chr  "26-100本" "11-25本" "10本以下" "200本以上" ...
# 看看前六筆
# 程式報表2.1
head(dta)
  國家 性別 父母教育  數學  化學  地科  生物  物理 電腦 書桌 書 房間 網路
1 台灣 女生     高中 729.4 624.8 607.8 650.7 618.3   有   無 有   無   有
2 台灣 女生     高中 776.2 661.8 665.7 627.9 616.7   有   有 有   無   有
3 台灣 女生 國小以下 718.2 708.8 629.3 617.5 663.3   有   有 無   無   有
4 台灣 女生     初中 607.2 633.0 620.3 592.1 616.5   有   無 無   無   有
5 台灣 女生 大學以上 658.2 617.8 614.1 560.2 591.0   有   無 無   無   有
6 台灣 女生     專科 478.6 601.2 591.0 598.8 554.0   有   有 有   有   有
       書量
1  26-100本
2   11-25本
3  10本以下
4 200本以上
5  26-100本
6  26-100本
#看看資料基本統計
#程式報表2.3
summary(dta)
     國家               性別             父母教育              數學    
 Length:15549       Length:15549       Length:15549       Min.   :136  
 Class :character   Class :character   Class :character   1st Qu.:546  
 Mode  :character   Mode  :character   Mode  :character   Median :609  
                                                          Mean   :603  
                                                          3rd Qu.:666  
                                                          Max.   :918  
      化學          地科          生物          物理         電腦          
 Min.   :228   Min.   :172   Min.   :196   Min.   :152   Length:15549      
 1st Qu.:510   1st Qu.:508   1st Qu.:510   1st Qu.:511   Class :character  
 Median :566   Median :563   Median :566   Median :568   Mode  :character  
 Mean   :564   Mean   :559   Mean   :560   Mean   :564                     
 3rd Qu.:621   3rd Qu.:614   3rd Qu.:614   3rd Qu.:621                     
 Max.   :910   Max.   :848   Max.   :829   Max.   :852                     
     書桌                書                房間               網路          
 Length:15549       Length:15549       Length:15549       Length:15549      
 Class :character   Class :character   Class :character   Class :character  
 Mode  :character   Mode  :character   Mode  :character   Mode  :character  
                                                                            
                                                                            
                                                                            
     書量          
 Length:15549      
 Class :character  
 Mode  :character  
                   
                   
                   
#載進 lattice，準備畫圖。
require(lattice)
載入需要的套件：lattice
#看看數學分數的直方圖
#圖2.1
histogram(~ 數學, data = dta, xlab = '數學分數', ylab='機率',type = "density")

##連續變項間關係
#把學科分數取出來
dta_scores <- dta[, c('數學', '化學', '地科', '生物', '物理')]
#兩兩變項畫散佈圖
#圖2.2
pairs(dta_scores, pch = '.', upper.panel = panel.smooth, lower.panel = NULL, 
      col = 'gray')

#利用 lattice 的 splom 指令重畫兩兩變項散佈圖，算是進階版
#圖2.3
splom(~ dta_scores, cex = 0.1, pch = '.', axis.text.cex = 0.5, 
      type = c('p', 'r', 'g'))

#數學與物理分數相關
round(cor(dta$數學,dta$物理), 3)
[1] 0.725
#所有學科分數相關
#程式報表2.4
round(cor(dta_scores), 3)
      數學  化學  地科  生物  物理
數學 1.000 0.723 0.694 0.711 0.725
化學 0.723 1.000 0.857 0.868 0.859
地科 0.694 0.857 1.000 0.876 0.861
生物 0.711 0.868 0.876 1.000 0.897
物理 0.725 0.859 0.861 0.897 1.000
#檢定相關是否顯著，也可以看到信賴區間
cor.test( ~ 數學 + 物理, data = dta_scores)

    Pearson's product-moment correlation

data:  數學 and 物理
t = 131, df = 15547, p-value <2e-16
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 0.7170 0.7319
sample estimates:
   cor 
0.7245 
#載進 Hmist，一次檢定所有相關
#程式報表2.5
require(Hmisc)
載入需要的套件：Hmisc
載入需要的套件：survival
載入需要的套件：Formula
載入需要的套件：ggplot2

載入套件：'Hmisc'
下列物件被遮斷自 'package:base':

    format.pval, units
rcorr(as.matrix(dta_scores), type="pearson")
     數學 化學 地科 生物 物理
數學 1.00 0.72 0.69 0.71 0.72
化學 0.72 1.00 0.86 0.87 0.86
地科 0.69 0.86 1.00 0.88 0.86
生物 0.71 0.87 0.88 1.00 0.90
物理 0.72 0.86 0.86 0.90 1.00

n= 15549 


P
     數學 化學 地科 生物 物理
數學       0    0    0    0  
化學  0         0    0    0  
地科  0    0         0    0  
生物  0    0    0         0  
物理  0    0    0    0       
#檢驗數學與物理、數學與生物何者相關高
#程式報表2.6
require(cocor)
載入需要的套件：cocor
cocor(~數學 + 物理 | 數學 + 生物, dta)

  Results of a comparison of two overlapping correlations based on dependent groups

Comparison between r.jk (數學, 物理) = 0.7245 and r.jh (數學, 生物) = 0.7108
Difference: r.jk - r.jh = 0.0137
Related correlation: r.kh = 0.8972
Data: dta: j = 數學, k = 物理, h = 生物
Group size: n = 15549
Null hypothesis: r.jk is equal to r.jh
Alternative hypothesis: r.jk is not equal to r.jh (two-sided)
Alpha: 0.05

pearson1898: Pearson and Filon's z (1898)
  z = 5.5418, p-value = 0.0000
  Null hypothesis rejected

hotelling1940: Hotelling's t (1940)
  t = 5.5631, df = 15546, p-value = 0.0000
  Null hypothesis rejected

williams1959: Williams' t (1959)
  t = 5.5544, df = 15546, p-value = 0.0000
  Null hypothesis rejected

olkin1967: Olkin's z (1967)
  z = 5.5418, p-value = 0.0000
  Null hypothesis rejected

dunn1969: Dunn and Clark's z (1969)
  z = 5.5518, p-value = 0.0000
  Null hypothesis rejected

hendrickson1970: Hendrickson, Stanley, and Hills' (1970) modification of Williams' t (1959)
  t = 5.5631, df = 15546, p-value = 0.0000
  Null hypothesis rejected

steiger1980: Steiger's (1980) modification of Dunn and Clark's z (1969) using average correlations
  z = 5.5498, p-value = 0.0000
  Null hypothesis rejected

meng1992: Meng, Rosenthal, and Rubin's z (1992)
  z = 5.5496, p-value = 0.0000
  Null hypothesis rejected
  95% confidence interval for r.jk - r.jh: 0.0182 0.0381
  Null hypothesis rejected (Interval does not include 0)

hittner2003: Hittner, May, and Silver's (2003) modification of Dunn and Clark's z (1969) using a backtransformed average Fisher's (1921) Z procedure
  z = 5.5493, p-value = 0.0000
  Null hypothesis rejected

zou2007: Zou's (2007) confidence interval
  95% confidence interval for r.jk - r.jh: 0.0088 0.0185
  Null hypothesis rejected (Interval does not include 0)
cocor(~數學 + 物理 | 地科 + 生物, dta)

  Results of a comparison of two nonoverlapping correlations based on dependent groups

Comparison between r.jk (數學, 物理) = 0.7245 and r.hm (地科, 生物) = 0.8764
Difference: r.jk - r.hm = -0.1519
Related correlations: r.jh = 0.6937, r.jm = 0.7108, r.kh = 0.8612, r.km = 0.8972
Data: dta: j = 數學, k = 物理, h = 地科, m = 生物
Group size: n = 15549
Null hypothesis: r.jk is equal to r.hm
Alternative hypothesis: r.jk is not equal to r.hm (two-sided)
Alpha: 0.05

pearson1898: Pearson and Filon's z (1898)
  z = -42.2316, p-value = 0.0000
  Null hypothesis rejected

dunn1969: Dunn and Clark's z (1969)
  z = -48.6493, p-value = 0.0000
  Null hypothesis rejected

steiger1980: Steiger's (1980) modification of Dunn and Clark's z (1969) using average correlations
  z = -48.1480, p-value = 0.0000
  Null hypothesis rejected

raghunathan1996: Raghunathan, Rosenthal, and Rubin's (1996) modification of Pearson and Filon's z (1898)
  z = -48.6493, p-value = 0.0000
  Null hypothesis rejected

silver2004: Silver, Hittner, and May's (2004) modification of Dunn and Clark's z (1969) using a backtransformed average Fisher's (1921) Z procedure
  z = -47.8431, p-value = 0.0000
  Null hypothesis rejected

zou2007: Zou's (2007) confidence interval
  95% confidence interval for r.jk - r.hm: -0.1590 -0.1449
  Null hypothesis rejected (Interval does not include 0)
##連續變項與類別變項間的關係
##兩個類別
#看看家裡有無電腦跟數學間的關係
#圖2.4
densityplot(~ 數學, groups = 電腦, data = dta, xlab = '數學分數', lty = c(1,2),
            plot.points = F, type = "g", , main = '電腦 (無 = 虛線, 有 = 實線)')

#也可用QQ圖比較
#圖2.5
qq(電腦 ~ 數學, data = dta, type = c('p', 'g'), pch = '.', aspect = 1, 
   xlab = '數學分數 (有電腦)', ylab = '數學分數 (無電腦)')

#看看有無電腦的學生數學平均與數標差
aggregate(數學 ~ 電腦, data = dta, FUN = mean)
  電腦  數學
1   有 605.5
2   無 532.7
aggregate(數學 ~ 電腦, data = dta, FUN = sd)
  電腦  數學
1   有 90.27
2   無 93.07
##多個類別
#看看不同國家學生的資料分數直方圖
#圖2.6
histogram(~ 數學 | 國家, data = dta, xlab = '數學分數', ylab='機率',
          type = 'density', layout = c(4, 1))

#以盒鬚圖看不同國家學生的差異
#圖2.7
bwplot(國家 ~ 數學, data = dta, xlab = '數學分數',
  panel = function(x, y, ...) {
    panel.bwplot(x, y, ...)
    panel.grid(v = -1, h = 0, ...) })

#看看不同國家的學生數學平均與數標差
aggregate(數學 ~ 國家, data = dta, FUN = mean)
  國家  數學
1 日本 573.7
2 台灣 616.4
3 香港 592.9
4 韓國 620.3
aggregate(數學 ~ 國家, data = dta, FUN = sd)
  國家   數學
1 日本  83.43
2 台灣 101.97
3 香港  79.97
4 韓國  86.36
#看看不同國家間，物理與數學間的關係是否類似
#圖2.8
xyplot(物理 ~ 數學 |  國家, data = dta, xlab = '數學分數', ylab = '物理分數',
       type = c("g", "p", "r"), cex = 0.1, layout = c(4, 1))

##連續變項與兩個類別變項間的關係
#看看不同國家、有無電腦的學生數學平均與對應的平均數標準誤
#程式報表2.7
show(m0 <- aggregate(數學 ~ 國家 + 電腦, data = dta, FUN = mean))
  國家 電腦  數學
1 日本   有 578.7
2 台灣   有 618.6
3 香港   有 593.2
4 韓國   有 620.7
5 日本   無 524.7
6 台灣   無 539.4
7 香港   無 542.0
8 韓國   無 577.0
show(m1 <- aggregate(數學 ~ 國家 + 電腦, data = dta, function(x) sd(x)/sqrt(length(x))))
  國家 電腦   數學
1 日本   有  1.470
2 台灣   有  1.519
3 香港   有  1.403
4 韓國   有  1.313
5 日本   無  4.645
6 台灣   無  9.955
7 香港   無 23.474
8 韓國   無 16.027
#把資料集中，並重排國家
m <- data.frame(m0, l = m0$數學 - m1$數學, u = m0$數學 + m1$數學)
m$國家 <- factor(m$國家, levels=c('日本', '香港', '台灣', '韓國'))
#載入 latticeExtra 套件
require(latticeExtra)
載入需要的套件：latticeExtra

載入套件：'latticeExtra'
下列物件被遮斷自 'package:ggplot2':

    layer
#在圖中加入了誤差
#圖2.9
segplot( 國家 ~ l + u | 電腦, data = m, centers = 數學, 
        draw.bands = F, xlab = '數學平均分數', ylab = '國家',
        main = '電腦', layout = c(1, 2),
        segments.fun = panel.arrows, ends = "both", angle = 90, length = 1, 
        unit = "mm")

3 In-class 3

xfun::pkg_load2(c("htmltools", "mime"))
xfun::embed_file("Animals.zip")

Download Animals.zip

4 In-class 4

xfun::embed_file("report_mos.zip")

Download report_mos.zip