1 In-class
#templates
需先下載prettydoc Package, 而後使用cayman
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 1641-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-4 分析概要報表
小數點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\).
方差分析表
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模型擬合圖
plot(women, xlab = "Height (in)", ylab = "Weight (lb)")
abline(m0, lty = 2)
grid()
殘差圖
檢查殘差分配有沒有規律
plot(resid(m0) ~ fitted(m0), xlab = "Fitted values",
ylab = "Residuals", ylim = c(-3.5, 3.5))
grid()
abline(h = 0, lty = 2)
驗證殘差常態分佈
qqnorm(resid(m0))
qqline(resid(m0))
grid()
結束
顯示演練單元信息
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] compiler_4.1.1 magrittr_2.0.1 fastmap_1.1.0 tools_4.1.1
## [5] htmltools_0.5.2 prettydoc_0.4.1 yaml_2.2.1 stringi_1.7.4
## [9] rmarkdown_2.11 highr_0.9 knitr_1.34 stringr_1.4.0
## [13] xfun_0.26 digest_0.6.28 rlang_0.4.11 evaluate_0.142 In-class 使用TIMSS 2011 年資料檔
dta <- read.table(file = "TIMSSmath.txt", header = TRUE, fileEncoding="big5")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.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.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#手動下載lattice
install.packages("lattice")## Warning: 正在使用 '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')
圖2.3 進階版
#利用 lattice 的 splom 指令重畫兩兩變項散佈圖，算是進階版
splom(~ dta_scores, cex = 0.1, pch = '.', axis.text.cex = 0.5,
type = c('p', 'r', 'g'))
#數學與物理分數相關
round(cor(dta$數學,dta$物理), 3)## [1] 0.7252.4程式報表-所有學科分數相關
3為小數點後幾位，可自修改
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.72452.5程式報表-一次檢定所有相關
#載進 Hmist
require(Hmisc)## 載入需要的套件：Hmisc## 載入需要的套件：survival## 載入需要的套件：Formula## 載入需要的套件：ggplot2##
## 載入套件：'Hmisc'## 下列物件被遮斷自 'package:base':
##
## format.pval, unitsrcorr(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 02.6程式報表 檢驗數學與物理、數學與生物何者相關高
require(cocor)## 載入需要的套件：cocorcocor(~數學+物理 | 數學+生物, 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 連續變項與類別變項
##兩個類別 #看看家裡有無電腦跟數學間的關係 type = “g”代表有網格
densityplot(~數學, groups=電腦, data=dta, xlab='數學分數', lty=c(1.2), plot.points = F, type = "g", , main = '電腦 (無 = 虛線, 有 = 實線)')
圖2.5 QQ圖比較
type = c(‘p’, ‘g’) —to get both the grid and points. aspect = 1– modify the width, pch = ‘.’–用點來呈現 ,pch =1 表示用圓圈呈現，不同數字有對應的符號
qq(電腦 ~ 數學, data = dta, type = c('p', 'g'), pch ='.' , aspect = 1,
xlab = '數學分數 (有電腦)', ylab = '數學分數 (無電腦)')
#看看有無電腦的學生數學平均與數標差 aggregate類似SAS的group by做分組後的計算
aggregate(數學 ~ 電腦, data = dta, FUN = mean)## 電腦 數學
## 1 有 605.5
## 2 無 532.7aggregate(數學 ~ 電腦, 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.3aggregate(數學 ~ 國家, 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.0show(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")