範例一:美國婦女的身高和體重
蘇睿寧
2021 10 16
套內數據
十五位三十到十九歲美國婦女的身高和體重.
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基本統計圖形
下面R程式碼畫出women數據集的散點圖.
plot(women, type = 'b', xlab = "Height (in)", ylab = "Weight (lb)")
grid()
線性模型分析
\[y_i = \beta_0 + \beta_1 x_i + \epsilon_i ,~~ \epsilon_i \sim N(0, \sigma^2) \]
體重 = 截距參數 + 斜率參數 x 身高 + 殘差(常態分佈)
分析概要報表
小數點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 17134)
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 yaml_2.2.1 jquerylib_0.1.4 stringi_1.7.5
[9] rmarkdown_2.11 highr_0.9 knitr_1.36 stringr_1.4.0
[13] xfun_0.26 digest_0.6.28 rlang_0.4.11 evaluate_0.14 #讀檔案，資料來自於 TIMSS 2011 年
dta <- read.table(file = "TIMSSmath.txt", header = TRUE)#看看前六筆 #程式報表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，準備畫圖。
library(lattice)
require(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
library(Hmisc)載入需要的套件：survival載入需要的套件：Formula載入需要的套件：ggplot2
載入套件：'Hmisc'下列物件被遮斷自 'package:base':
format.pval, unitslibrary(lattice)
library(survival)
library(Formula)
library(ggplot2)
require(Hmisc)
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
library(cocor)
require(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 = '電腦 (無 = 虛線, 有 = 實線)')
densityfunction (x, ...)
UseMethod("density")
<bytecode: 0x0000000012daf698>
<environment: namespace:stats>#也可用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.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 套件
library(latticeExtra)
載入套件：'latticeExtra'下列物件被遮斷自 'package:ggplot2':
layerrequire(latticeExtra)#在圖中加入了誤差 #圖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")