Bài tập về nhà
Nguyễn Thanh Hùng
MSSV: 2121013044
1 Tuần 6
1.1 Scraping dữ liệu
1.1.1 Kiểm tra trang web có cho phép cạo dữ liệu hay không.
Để kiểm tra chúng ta có được phép lấy dữ liệu từ 1 trang web hay
không chúng ta dùng package robotxtst
library(robotstxt)
paths_allowed(paths="https://www.investing.com/indices/vn")##
www.investing.com## [1] TRUEKết quả trả về TRUE chứng tỏ trang web investing.com cho
phép chúng ta scraping dữ liệu.
1.1.2 Scraping dữ liệu của chỉ số VNINDEX.
Để thực hiện scraping dữ liệu trên một trang web bất kỳ, ta cần phải
gọi package rvest xuất hiện như sau:
library(rvest)
library(stringr)
library(DT)
library(tidyverse)## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr 1.1.1 ✔ purrr 1.0.1
## ✔ forcats 1.0.0 ✔ readr 2.1.4
## ✔ ggplot2 3.4.2 ✔ tibble 3.2.1
## ✔ lubridate 1.9.2 ✔ tidyr 1.3.0
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ readr::guess_encoding() masks rvest::guess_encoding()
## ✖ dplyr::lag() masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errorsLấy dữ liệu giá chứng khoán của chỉ số VNI:
url <- "https://www.investing.com/indices/vn-historical-data"
webpage <- read_html(url)
data <- webpage%>%html_table(fill = TRUE)
vni = data[[2]]
vni <- as.data.frame(vni)Sau khi cào xong dữ liệu chỉ số VNI, ta có các biến như sau:
Date: Ngày chỉ số VNI được thu thập
Price: Giá đóng cửa của từng ngày
Open: Giá mở cửa của từng phiên
High: Giá cao nhất trong ngày
Giá thấp nhất trong ngày
Vol: khối lượng giao dịch trong ngày
Change: Thay đổi về tỷ suất sinh lợi (%)
Xử lí dữ liệu: Cột Date chưa có định dạng ngày nên ta
thêm 1 biến mới làdate chuyển về định dạng ngày nhằm có thể
vẽ biểu đồ thể hiện xu hướng giá chứng khoán.
vni$date <- as.Date(vni$Date,format="%m/%d/%Y")Cột price chứa dấu (,) nên cấu trúc của cột Price đang là
chuỗi. Do đó chúng ta chuyển hết dữ liệu về số
như sau
pvni <- vni$Price
pricevni <- as.numeric(gsub(",", "", pvni))
pricevni## [1] 1173.13 1168.40 1165.42 1154.20 1151.77 1149.02 1138.07 1126.22 1134.62
## [10] 1132.00 1125.50 1120.18 1125.39 1138.35 1134.33 1132.03 1129.38 1125.30
## [19] 1118.46 1111.72 1105.40vni$Price <- pricevni1.1.3 Biểu đồ thể hiện xu hướng.
Trước hết chúng ta vẽ đồ thị thể hiện đường giá của chuỗi VNI
library(ggplot2)
ggplot(data=vni,aes(x=date,y=Price))+geom_line(color="blue")
Kết quả cho chúng ta thấy trend của chỉ số VNI trong ngắn hạn, chúng ta có thể đoán chỉ số VNI đang trong xu hướng tăng. Điều này cũng dễ hiểu vì gần đây ngân hàng trung ương đã hạ lãi suất nên thanh khoản chủa thị trường chứng khoán cụ thể là chỉ số VNI đã được cải thiện rõ rệt.
1.1.4 Kiểm định phân phối chuẩn (Jarque-Bera)
Để kiểm định một chuỗi dữ liệu có phải là phân phối chuẩn hay không, chúng ta có rất nhiều phương pháp. Trong phần này, ta sử dụng kiểm định Jarque-Bera để kiểm tra điều đó.
Ta có bài toán kiểm định
H0: Chuỗi VNI có phân phối chuẩn
H1: Chuỗi VNI không có phân phối chuẩn
library(tseries)## Registered S3 method overwritten by 'quantmod':
## method from
## as.zoo.data.frame zoojarque.bera.test(vni$Price)##
## Jarque Bera Test
##
## data: vni$Price
## X-squared = 1.1725, df = 2, p-value = 0.5564Ta thấy giá trị p_value =0.55 > mức ý nghĩa. Do đó chấp nhận H0 hay chuỗi VNI có phân phối chuẩn.
ggplot(data=vni,aes(x=Price))+geom_histogram(fill="blue")## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
Đây là biểu đồ thể hiện phân phối của chuỗi giá VNI.
1.1.5 Kiểm định tính dừng(Augmented Dickey–Fuller)
Ta có bài toán kiểm định như sau:
Giá định:
- H0: Chuỗi VNI không dừng
- H1: Chuỗi VNI dừng
Thực hiện kiểm định đối với chuỗi lợi suất VNI như sau:
library(tseries)
adf.test(vni$Price, alternative = "stationary")##
## Augmented Dickey-Fuller Test
##
## data: vni$Price
## Dickey-Fuller = -2.0676, Lag order = 2, p-value = 0.5466
## alternative hypothesis: stationaryTa thấy p_value= 0.54 > mức ý nghĩa. Do đó chấp nhận H0. Vậy chuỗi chỉ số VNI không dừng. Vậy hiện tại dữ liệu VNI không được phép sử dụng để chạy mô hình. Ta phải thực hiện một số thủ thuật để đưa về chuỗi dừng được đề cập ở bên dưới.
1.1.6 Kiểm định tự tương quan chuỗi Ljung_box
Ta có bài toán kiểm định như sau:
H0: Không có tương quan chuỗi
H1: Có tương quan chuỗi
Thực hiện kiểm định chuỗi chỉ số VNI như sau:
library(tseries)
Box.test(vni$Price,lag=20, type="Ljung-Box")##
## Box-Ljung test
##
## data: vni$Price
## X-squared = 111.52, df = 20, p-value = 1.033e-14Ta thấy p_value xấp xỉ 0% do đó bác bỏ H0, nghĩa là chỉ số VNI có tương quan chuỗi với nhau ở độ trễ 20 ngày.
1.1.7 Đưa dữ liệu không có tính dừng sang chuỗi dữ liệu dừng
Khi dữ liệu không có tính dừng, chúng ta sẽ không đưa được vào các mô hình để phân tích. Do đó, chúng ta phải biến dữ liệu gốc về dạng return hoặc log return. Mục đích là để cho sự biến động của chuỗi VNI được kiểm soát và có tính dừng. Ta thực hiện như sau:
library(tidyverse)
log_returns <- vni$Price %>%
lag() %>% # Lấy giá trị trước đó
log() - log(vni$Price) # Tính toán log return
log_returns <- na.omit(log_returns)Sau một loạt các câu lệnh, ta đã có một vecto
log_returns không còn chứa dữ liệu giá mà đã được chuyển
sang dữ liệu suất sinh lợi. Tiếp theo ta chuyển vector
log_returns thành dataframe và tạo bộ dữ liệu mới
retvni gồm 2 biến date chỉ ngày và
rvni chứa log_retureturns
rvni <- as.data.frame(log_returns)
date <- vni$date[c(1:20)]
retvni <- data.frame(date=date,rvni=log_returns)Bây giờ chúng ta thử vẽ biểu đồ của dataset retvni xem liệu rằng nó đẫ thay đổi hay chưa.
ggplot(data=retvni,aes(x=date,y=rvni))+geom_line(col="red")
Biểu đồ trên thể hiện suất sinh lợi của chỉ số VNI. Do có ít quan sát nên biểu đồ chưa thể dự đoán được retvni đã có tính dừng.
1.1.8 Thống kê mô tả
Thống kê mô tả chuỗi lợi suất retvni bằng gói
fBasics
library("fBasics")
basicStats(vni$Price)## X..vni.Price
## nobs 21.000000
## NAs 0.000000
## Minimum 1105.400000
## Maximum 1173.130000
## 1. Quartile 1125.390000
## 3. Quartile 1149.020000
## Mean 1136.137619
## Median 1132.030000
## Sum 23858.890000
## SE Mean 3.979064
## LCL Mean 1127.837438
## UCL Mean 1144.437800
## Variance 332.491889
## Stdev 18.234360
## Skewness 0.493574
## Kurtosis -0.696426Chúng ta thấy chuỗi giá vni bao gồm 20 quan sát và giá trị lớn nhất là 1173.13, giá trị nhỏ nhất là 1105.40. Tứ phân vị thứ nhất là 1125,39 cho biết 25% quan sát có giá trị dưới 1123,39 và 75% quan sát có giá trị trên 1123,39. Tứ phân vị thứ 3 là 1149.02 cho biết có 75% quan sát có giá trị dưới 1149,02 và 25% trên ngưỡng đó. Trung bình đạt 1136,13, trung vị chỉ số VNi đạt 1132.03…. Độ lệch chuẩn đạt 18.23 cho biết mức đôj phân tán của chuỗi VNI quanh giá trị trung bình.
1.2 Viết function thực hiện các công việc
Ta sử dụng dữ liệu tỷ suất sinh lợi của 5 thị trường tiền ảo bao gồm Bicoin, Ethereum, EOS, XRP và BitcoinCash đã được tôi thu thập từ trước. Các biến bao gồm:
DATE: Ngày nghiên cứu: từ 1/1/2018 - 30/1/2023 với 1855 quan sát
BTC: Tỷ suất sinh lợi thị trường Bitcoin (%)
ETH: Tỷ suất sinh lợi thị trường Ethereum (%)
BCH: Tỷ suất sinh lợi thị trường BitcoinCash (%)
EOS: Tỷ suất sinh lợi thị trường EOS (%)
XRO: Tỷ suất sinh lợi thị trường XRP(%)
library(readxl)
library(DT)
crypto <- read_excel("D:\\data nckh\\datareturn.xlsx")
datatable(crypto)Sau đây tôi sẽ viết một function dùng để tính toán thống kê mô tả các đặc trưng đo lường, hệ số độ nhọn (Skewness) cho biết độ lệch của phân phối, hệ số độ lồi(Kurtosis) cho biết phân phối có đuôi phình to hay không, các giá trị p_value của các kiểm định như là kiểm định phân phối chuân (JarqueBera), kiểm định tính dừng (Augmented Dickey–Fuller), Kiểm định tự tương quan (Ljung_Box) với độ trễ là 20 ngày.
fun1 <- function(data,var) {
library(tidyverse)
library(tseries)
data %>% summarize(
n=n(),
mean = mean({{var}},na.rm = TRUE),
median = median({{var}},na.rm = TRUE),
quartiles_25 = quantile({{var}}, c(0.25)),
quartiles_75 = quantile({{var}}, c(0.75)),
sd = sd({{var}}),
variance = var({{var}}),
max_value = max({{var}}),
min_value = min({{var}}),
skewness = moments::skewness({{var}}),
kurtosis = moments::kurtosis({{var}}),
jarque_bera_p_value = jarque.bera.test({{var}})$p.value,
adf_p_value=adf.test({{var}})$p.value,
ljung_p_value=Box.test({{var}},lag=20,type=c("Ljung-Box"))$p.value
)
}Tiến hành sử dụng function đã được lập từ trước tính toán thống kê mô tả với chuỗi tỷ suất sinh lợi của Bitcoin như sau:
btc <- fun1(crypto,crypto$BTC)## Warning: There was 1 warning in `summarize()`.
## ℹ In argument: `adf_p_value = adf.test(crypto$BTC)$p.value`.
## Caused by warning in `adf.test()`:
## ! p-value smaller than printed p-valuedatatable(btc,caption=c("Thống kê mô tả và kiểm định chuỗi lợi suất Bitcoin"))Ta thấy bộ dữ liệu tỉ suất sinh lợi Bitcoin có 1855 quan sát, với giá trị trung bình là 0.02%. Trung vị là 0.09% cho biết 50% quan sát về tỉ suất sinh lợi của bitcoin dưới 0.09%, 50% quan sát trên 0.09%. Tứ phân vị thứ nhất đạt -1.63% cho biết 25% tỉ suất sinh lợi của Bitcoin đạt dưới -1.63% và 75% tỉ suất sinh lợi của Bitcoin đạt trên -1.63%. Độ lệch chuẩn cho biến độ biến động của tỷ suất sinh lợi Bitcoin là 3.94%. Giá trị tỷ suất sinh lợi lớn nhất đạt 17.86%, nhỏ nhất là -48.09%. Hệ số Skewness là -1.13 < 1 cho thấy phân phối tỷ suất sinh lợi lệch về bên trái so với phân phối chuẩn. Hệ số kurtosis đạt 16.89 >3 cho thấy phân phối tỷ suất sinh lợi Bitcoin có đuôi phình to hơn phân phối chuẩn. p_value của thống kê Jarque_bera bằng 0 và nhỏ hơn mức ý nghĩa, do đó tỷ suất sinh lợi BTC không có phân phối chuẩn. P_value của kiểm định ADF_TEST là 0.01% nhỏ hơn mức ý nghĩa, do đó bác bỏ H0, nghĩa là chuỗi dữ liệu bitcoin dừng tại mức ý nghĩa 1%, điều này cho thấy dữ liệu phù hợp để đưa vào các mô hình. Cuối cùng, p_value của kiểm định Ljung_box là 0.3 lớn hơn mức ý nghĩa, chấp nhận H0, nghĩa là chuỗi dữ liệu Bitcoin không có hiện tượng tự tương quan chuỗi ở độ trễ 20 ngày.
Thực hiện tương tự cho chuỗi tỷ suất sinh lợi của thị trường Ethereum
eth <- fun1(crypto,crypto$ETH)## Warning: There was 1 warning in `summarize()`.
## ℹ In argument: `adf_p_value = adf.test(crypto$ETH)$p.value`.
## Caused by warning in `adf.test()`:
## ! p-value smaller than printed p-valuedatatable(eth,caption=c("Thống kê mô tả và kiểm định chuỗi lợi suất Ethereum"))Ta thấy bộ dữ liệu tỉ suất sinh lợi Ethereum có 1855 quan sát, với giá trị trung bình là 0.03%. Trung vị là 0.11% cho biết 50% quan sát về tỉ suất sinh lợi của Ethereum dưới 0.11%, 50% quan sát trên 0.11%. Tứ phân vị thứ nhất đạt -2.22% cho biết 25% tỉ suất sinh lợi của Ethereum đạt dưới -2.22% và 75% tỉ suất sinh lợi của Bitcoin đạt trên -2.22%. Độ lệch chuẩn cho biến độ biến động của tỷ suất sinh lợi Ethereum là 5.16%. Giá trị tỷ suất sinh lợi lớn nhất đạt 23.07%, nhỏ nhất là -58.96%. Hệ số Skewness là -1.08 < 1 cho thấy phân phối tỷ suất sinh lợi lệch về bên trái so với phân phối chuẩn. Hệ số kurtosis đạt 14.55 >3 cho thấy phân phối tỷ suất sinh lợi Ethereum có đuôi phình to hơn phân phối chuẩn. p_value của thống kê Jarque_bera bằng 0 và nhỏ hơn mức ý nghĩa, do đó tỷ suất sinh lợi ETH không có phân phối chuẩn. P_value của kiểm định ADF_TEST là 0.01% nhỏ hơn mức ý nghĩa, do đó bác bỏ H0, nghĩa là chuỗi dữ liệu etherum dừng tại mức ý nghĩa 1%, điều này cho thấy dữ liệu phù hợp để đưa vào các mô hình. Cuối cùng, p_value của kiểm định Ljung_box là 0.002 nhỏ hơn mức ý nghĩa, bác bỏ H0, nghĩa là chuỗi dữ liệu Ethereum có hiện tượng tự tương quan chuỗi trong vòng 20 ngày (độ trễ 20 ngày.)
1.3 Tìm hiểu các Indicator
1.3.1 GDP (Gross Domestic Product)
GDP (Gross Domestic Product) - Tổng sản phẩm quốc nội: Tổng giá trị của tất cả hàng hóa và dịch vụ được sản xuất trong một quốc gia trong một năm. Ta thực hiện lấy dữ liệu này trên worldbank như sau:
library(WDI)
gdp_data <- WDI(country = "all", indicator = "NY.GDP.MKTP.CD", start = 1990, end = 2020)
gdp_data <- na.omit(gdp_data)Các biến có trong gdp_data là:
iso2c: Mã quốc gia 2 chữ cái (ISO 3166-1 alpha-2 code).
iso3c: Mã quốc gia 3 chữ cái (ISO 3166-1 alpha-3 code).
country”: Tên quốc gia.
year”: Năm.
NY.GDP.MKTP.CD”: Tổng sản phẩm quốc nội (GDP) theo giá trị thị trường (USD).
Lấy dữ liệu GDP của Brazil như sau:
Brazil <- gdp_data[gdp_data$country=="Brazil",]Vẽ biểu đồ đường thể hiện xu hướng sự thay đổi GDP của Brazil qua từng năm:
library(ggplot2)
ggplot(data=Brazil,aes(x=year,y=NY.GDP.MKTP.CD))+geom_line(color="blue")
Nhìn vào biểu đồ ta thấy GDP của Braxzil có xu hướng tăng lên từ nă 2000 đến 2010 và có xu hướng giảm xuống từ năm 2011 đến năm 2020.
Tính toán chỉ tiêu tốc độ phát triển qua từng năm của GDP Brazil
log_returns <- Brazil$NY.GDP.MKTP.CD %>%
lag() %>% # Lấy giá trị trước đó
log() - log(Brazil$NY.GDP.MKTP.CD) # Tính toán log return
log_returns <- as.data.frame(log_returns)
rate <- na.omit(log_returns)
Year <- Brazil$year[c(1:30)]
rbrazil <- data.frame(rate,Year)Ta có dataframe rbrazil có biến rate là tốc độ phát triển GDP của Brazil theo từng năm. Ta vẽ biểu đồ đường như sau:
ggplot(data=rbrazil,aes(x=Year,y=log_returns))+geom_line(color="blue")
Biểu đồ tốc độ phát triển cho chúng ta thấy dường như data tốc độ phát triển của GDP Brazil có tính dừng.
2 Tuần 5 : TỔNG HỢP
Trước hết ta thực hiện xử lí dữ liệu như sau:
library('palmerpenguins') #gọi package 'palmerpenguins' xuất hiện
h <- penguins # gán penguins vào h
h <- na.omit(h) #loại bỏ các giá trị NA
names(h) <- c("s","i","bl","bd","f","bm","sx","y") #đặt tên cho các biến có trong dữ liệu
library(DT)
datatable(h)library(utf8)
library(tidyverse)2.1 Trực quan hóa đặc trưng đo lường giữa các biến.
2.1.1 Trung bình
Đầu tiên ta nhóm bộ dữ liệu chim cánh cụt theo loài và tính trung bình cho biến bill_leng_mm(độ dài mỏ). Gán kết quả vào tb.
tb <- h %>% group_by(s) %>% summarise(m=mean(bl))
tb## # A tibble: 3 × 2
## s m
## <fct> <dbl>
## 1 Adelie 38.8
## 2 Chinstrap 48.8
## 3 Gentoo 47.6Kết quả loài Adelie có độ dài mỏ trung bình là 38.82mm. Loài Chinstrap có độ dài mỏ trung bình là 48.83mm. Loài Gentoo có độ dài mỏ trung bình là 47.56mm.
Tiếp theo ta chuyển tb thành một dataframe để thực hiện vẽ đồ thị cho dataframe tb.
tb <- as.data.frame(tb)
library(ggplot2)
ggplot(data=tb,aes(x=s,y=m))+geom_col(fill="purple")
Ta thấy loài Adelie có trung bình độ dài mỏ thấp nhất, loài Gentoo có trung bình độ dài mỏ đứng thứ 2, cuối cùng loài Chinstrap có trung bình độ dài mỏ lớn nhất.
2.1.2 Độ lệch chuẩn
Tiếp theo ta nhóm bộ dữ liệu chim cánh cụt theo loài và tính độ lệch chuẩn độ dài cánh của từng loài.
sd <- h %>% group_by(s) %>% summarise(dolech=sd(f))
sd## # A tibble: 3 × 2
## s dolech
## <fct> <dbl>
## 1 Adelie 6.52
## 2 Chinstrap 7.13
## 3 Gentoo 6.59Kết quả cho thấy cánh của loài Adelie có độ lệch chuẩn là 6.52. Độ dài cánh của loài Chinstrap có độ lệch chuẩn là 7.13 và cuối cùng độ lệch chuẩn cánh của loài Gentoo là 6.58. Qua đó cho thấy độ dài cánh loài Chinstaps có độ lệch chuẩn lớn nhất tương đương mức độ phân tán của độ dài cánh quanh trung bình là lớn nhất.
sd <- as.data.frame(sd)
ggplot(data=sd,aes(x=s,y=dolech))+geom_col(fill="blue")+
geom_text(aes(label=dolech),color = 'yellow', vjust = 2, size = 5)+
labs(x="loài",y="Độ lệch chuẩn cánh")
Ta thấy rằng loài Chinstrap có độ lệch chuẩn của cánh lớn nhất, tương đương sự phân tán chiều dài cánh quan giá trị trung bình là lớn nhất. Sụ phân tán độ dài cánh của loài Adelie và Gentoo quanh giá trị trung bình khá giống nhau.
2.1.3 Phương sai
Tiếp theo ta thực hiện tính phương sai của biến bill_length_mm(độ dài mỏ)
varbl <- h %>% group_by(i) %>% summarise(lech2=var(bl))
varbl## # A tibble: 3 × 2
## i lech2
## <fct> <dbl>
## 1 Biscoe 23.3
## 2 Dream 35.4
## 3 Torgersen 9.17Tiếp theo ta biến varbl thành dataframe và vẽ đồ thị
varbl <- as.data.frame(varbl)
ggplot(data=varbl,aes(x=i,y=lech2))+ geom_col(fill="blue")+
geom_text(aes(label=lech2),color="white",vjust=2,size=4)+
labs(x="Nơi sinh sống",y="Phương sai độ dài mỏ")
Ta thấy độ biến động độ dài mỏ của cá thể cánh cụt sống ở đảo Dream là lớn nhất với phương sai xấp xỉ 35.37. Độ biến động độ dài mỏ của cá thể sống ở Torgersen là nhỏ nhất với phương sai xấp xỉ 9.17.
2.2 Phân tổ không đều, lọc và vẽ biểu đồ
summary(h$f)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 172 190 197 201 213 231daicanh <- cut(h$f,breaks=c(0,210,231))
table(daicanh)## daicanh
## (0,210] (210,231]
## 237 96Chúng ta thấy rằng có 237 cá thể cánh cụt có chiều dài từ (0-210mm], có 96 cá thể có chiều dài cánh từ (210-231mm]. Chúng ta trực quan hóa điều này lên biểu đồ phân tán như sau:
h %>% ggplot(aes(x=bl,y=f))+
geom_point(color="red",size=3)+
geom_point(data = h %>% filter(f > 210 ), size = 3, color = 'blue') +
facet_wrap(~s)+
xlab('Độ dài mỏ')+
ylab('Độ dài cánh')
Ta thấy rằng 237 cá thể cánh cụt có chiều dài cánh từ (0-210mm] được chấm bằng màu đỏ và 96 cá thể có chiều dài cánh từ (210-237mm) được chấm bằng màu xanh nước biển trên đồ thị phân tán. Phân biểu đồ theo từng loài ta có thể thấy rằng hầu như cá thể cánh cụt có chiều dài cánh lớn hơn 210mm đều là loài Gentoo.
Trực quan hóa điều này lên đồ thị phân tán như sau:
h %>% ggplot(aes(x=bl,y=f,shape=i))+
geom_point(color="red",size=3)+
geom_point(data = h %>% filter(f > 210 ), size = 3, color = 'blue') +
xlab('Độ dài mỏ')+
ylab('Độ dài cánh')
Tiếp theo tôi lọc ra các đặc điểm các cá thể chim cánh cụt được nghiên cứu vào năm 2007
thu2007 <- h[h$y==2007,]
str(thu2007)## tibble [103 × 8] (S3: tbl_df/tbl/data.frame)
## $ s : Factor w/ 3 levels "Adelie","Chinstrap",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ i : Factor w/ 3 levels "Biscoe","Dream",..: 3 3 3 3 3 3 3 3 3 3 ...
## $ bl: num [1:103] 39.1 39.5 40.3 36.7 39.3 38.9 39.2 41.1 38.6 34.6 ...
## $ bd: num [1:103] 18.7 17.4 18 19.3 20.6 17.8 19.6 17.6 21.2 21.1 ...
## $ f : int [1:103] 181 186 195 193 190 181 195 182 191 198 ...
## $ bm: int [1:103] 3750 3800 3250 3450 3650 3625 4675 3200 3800 4400 ...
## $ sx: Factor w/ 2 levels "female","male": 2 1 1 1 2 1 2 1 2 2 ...
## $ y : int [1:103] 2007 2007 2007 2007 2007 2007 2007 2007 2007 2007 ...
## - attr(*, "na.action")= 'omit' Named int [1:11] 4 9 10 11 12 48 179 219 257 269 ...
## ..- attr(*, "names")= chr [1:11] "4" "9" "10" "11" ...Kết quả thu được tệp dữ liệu thu2007 với 103 hàng cà 8 cột, tức là có 103 cá thể được nghiên cứu vào năm 2007. Trực quan hóa điều này lên đồ thị phân tán như sau:
h %>% ggplot(aes(x=bl,y=f,shape=i))+
geom_point(color="red",size=3)+
geom_point(data = h %>% filter( y==2007 ), size = 3, color = 'blue') +
facet_wrap(~s)+
xlab('Độ dài mỏ')+
ylab('Độ dài cánh')
Ta vẽ đồ thị phân tán giữa biến chiều dài mỏ và độ dài cánh. Tiếp theo ta vẽ thêm 1 lớp đồ thị mới lọc ra những cá thể chim cánh cụt được nghiên cứu vào năm 2007 và tô màu xanh. Kết quả ta thấy những chấm xanh trên đồ thị là những cá thể được nghiên cứu vào năm 2007
Tiếp theo ta lọc ra những cá thể chim cánh cụt sống ở đảo Biscoe như sau:
Biscoe <- filter(h,i=="Biscoe")
datatable(Biscoe)TRực quan hóa điều này trên đồ thị như sau:
h %>% ggplot(aes(x=f,y=bm))+
geom_point(color="red",size=3)+
geom_point(data = h %>% filter( i=="Biscoe" ), size = 3, color = 'yellow') + facet_wrap(~s)+
xlab('Độ dài cánh')+
ylab('Khối lượng cơ thể')
Biểu đồ trên là biểu đồ phân tán giữa 2 biến độ dài cánh và khối lượng cơ thể chim cánh cụt. Sau khi vẽ xong, ta đã thêm một layer mới lọc ra những chú chim cánh cụt sống ở đảo Biscoe. Kết quả những chấm màu vàng thể hiện điều đó. Và ta thấy rằng hầu như loài Gentoo đều sống ở đảo Biscoe. Một số cá thể loài Adelie cũng sống ở đảo này. Ngoài ra loài Chinstrap hầu như không sống ở đảo Biscoe.
2.3 Đồ thị histogram, đồ thị phân tán và hệ số tương quan
Trong phần này tôi thực hiện tính ma trận tương quan giữa các biến. Ngoài ra tôi còn vẽ đồ thị thể hiện phân phối giữa các biến và đồ thị phân tán. Các điều này được thể hiện trên duy nhất 1 hình ảnh.
2.3.1 Hai biến độ dài cánh (f) và khối lượng cơ thể (bm)
names(h)## [1] "s" "i" "bl" "bd" "f" "bm" "sx" "y"attach(h)
fvsbl <- cbind(f,bm)
library(psych)##
## Attaching package: 'psych'## The following object is masked from 'package:fBasics':
##
## tr## The following objects are masked from 'package:ggplot2':
##
## %+%, alphapairs.panels(fvsbl)
Ta thấy hệ số tương quan giữa 2 biến độ dài cánh và khối lượng cơ thể rất cao, bằng 0.87. Qua đó, có thể kết luận 2 biến này tương quan với nhau khá mạnh. Đường chéo thể hiện phân phối của 2 biến này. Bên trái đường chéo đó chính là đồ thị phân tán và đường hồi quy.
2.3.2 Ba biến độ dài cánh (f), độ dài mỏ (bl) và độ sâu mỏ (bd)
Tiếp theo ta vẽ ma trận tương quan giữa 3 biến độ dài cánh, độ dài mỏ và độ sâu mỏ như sau:
fblbd <- cbind(f,bl,bd)
pairs.panels(fblbd)
Ta thấy hệ số tương quan giữa 2 biến độ dài cánh (f) và độ dài mỏ khá cao, bằng 0.65. Qua đó, có thể kết luận 2 biến này tương quan thuận với nhau khá mạnh. Hệ số tương quan giữa 2 biến bl và bd khá thấp, bằng -0.23 cho thấy 2 biến độ dài mỏ và độ sâu mỏ tương quan nghịch với nhau và tương quan khá yếu. Đường chéo thể hiện phân phối của 3 biến này. Bên trái đường chéo đó chính là đồ thị phân tán và đường hồi quy giữa các biến.
2.3.3 Tất cả các biến trong dataset penguins
Tiếp theo ta vẽ tiếp ma trận tương quan giữa tất cả các biến định lượng trong bộ dữ liệu penguins, như sau:
tatca = cbind(f,bl,bd,bm)
pairs.panels(tatca)
Bên phải đường chéo chính là hệ số tương quan. Dựa vào ma trận tương quan ta có thể thấy các biến có tương quan nghich với nhau bao gồm: tương quan khá yếu giữa độ dài mỏ(bl) và độ sâu mỏ(bd) với hệ số tương quan là -0.23, độ sâu mỏ(bd) và độ dài cánh(f) với hệ số tương quan -0.58,độ sâu mỏ (bd) và khối lương cơ thể(bm) với hệ số tương quan là -0.47. Tiếp theo, các biến có tương quan thuận với nhau bao gồm: độ dài mỏ (bl) và độ dài cánh(f) với hệ số tương quan là 0.65, độ dài mỏ (bl) và khối lượng cơ thể(bm) với hệ số tương quan là 0.59, độ dài cánh và(f) khối lượng cơ thể(bm) với hệ số tương quan là 0.87. Trong đó 2 biến f và bm có tương quan mạnh nhât, hai biến bl và bd có tương quan yếu nhất.
2.4 Kiểm định sự tương quan
Ta có bài toán kiểm định như sau:
Bài toán kiểm định
- H0: Không có tương quan (R=0)
- H1: Có sự tương quan (R khác 0)
2.4.1 Hai biến độ dài cánh và khối lượng cơ thể
Thực hiện kiểm định 2 biến độ dài cánh và khối lượng cơ thể:
cor.test(h$f,h$bm)##
## Pearson's product-moment correlation
##
## data: h$f and h$bm
## t = 32.562, df = 331, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.8447622 0.8963550
## sample estimates:
## cor
## 0.8729789Ta thấy p_value < 2e10^-16 < 5%, bác bỏ H0. Vậy 2
biến body_mass_g và flipper_length_mm có tương
quan tuyến tính với nhau. 2,2e*10^-16 <
5%, bác bỏ H0. Vậy 2 biếnbody_mass_gvàflipper_length_mm`
có tương quan tuyến tính với nhau.
Vẽ đồ thị phân tán giữa 2 biến này như sau:
library(tidyverse)
h %>% ggplot(aes(x=f,y=bm))+geom_point(color="purple",size=2)+
geom_smooth(formula=y~x,method='lm',color='yellow')+
xlab('Độ dài mỏ')+
ylab('Khối lượng cơ thể')
2.4.2 Hai biến độ dài cánh và độ dài mỏ
Thực hiện kiểm định giữa 2 biến độ dài cánh (f) và độ dài mỏ (bl) như sau:
cor.test(h$f,h$bl)##
## Pearson's product-moment correlation
##
## data: h$f and h$bl
## t = 15.691, df = 331, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.5868092 0.7106869
## sample estimates:
## cor
## 0.6530956Ta thấy p_value < 2.2e-16 < 5%, bác bỏ H0. Vậy 2
biến bill_length_mmmm và flipper_length_mm có
tương quan tuyến tính với nhau.
Thực hiện vẽ đồ thị phân tán và đường hồi quy như sau:
library(tidyverse)
h %>% ggplot(aes(x=f,y=bl))+geom_point(color="orange",size=2)+
geom_smooth(formula=y~x,method='lm',color='blue')+
xlab('Độ dài cánh')+
ylab('Độ dài mỏ')
2.4.3 Hai biến độ dài mỏ và khối lượng cơ thể
Thực hiện kiểm định giữa 2 biến độ dài mỏ và khối lượng cơ thể như sau:
cor.test(h$bl,h$bm)##
## Pearson's product-moment correlation
##
## data: h$bl and h$bm
## t = 13.276, df = 331, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.5145745 0.6554058
## sample estimates:
## cor
## 0.5894511Ta thấy p_value < 2.2e-16 < 5%, bác bỏ H0. Vậy 2
biến bill_length_mmmm và body_mass_g có tương
quan tuyến tính với nhau.
Vẽ đồ thị phân tán và đường hồi quy giữa 2 biến này như sau:
h %>% ggplot(aes(x=bl,y=bm))+geom_point(color="red",size=2)+
geom_smooth(formula=y~x,method='lm',color='yellow')+
xlab('Độ dài mỏ')+
ylab('Khối lượng cơ thể')
2.5 Tạo dataframe mới
Tiếp theo, giả định rằng dữ liệu penguins được 2 người thu thập với tổng cộng 333 cột(không tính các cột missing value) . Tôi sẽ tạo ra một dataframe có giá trị là 1 và 2 lần lượt tương ứng dữ liệu do người thứ nhất và người thứ hai thu thập.
Dữ liệu từ cột 1 đến cột 146 do người thứ nhất thu thập, được mã hóa bởi số 1
Dữ liệu từ cột 147 đến cột 333 do người thứ hai thu thập, được mã hóa bởi số 2
Quy trình thực hiện như sau:
Đầu tiên gọi biến
svà gán vào da (mỗi một phần tử trong biếnstương ứng với 1 cá thể chim cánh cụt)Tiếp theo lấy 146 phần tử đầu tiên (do người thứ nhất thu thập) gán vào biến
t1, phần từ tiếp theo từ 147 đến 333 (do người thứ 2 thu thập) gán vào biếnt2.tiếp theo tạo một data frame
au1gồm 2 biếnloai(chứa các phần tử của biến s từ 1-150) vàauthor(chứa số 1)Tương tự tạo data frame au2 gồm 2 biến
loai(chứa các phần tử của biếnspeciestừ 151-333) và author (chứa số 2 tương ứng do người thứ 2 thực hiện)Gộp 2 dataframe
au1au2lại bằng lệnh rbind() và gán vàoau. Kết quả tạo ra dataframe mới làauthỏa mãn yêu cầu.
da <- h$s
t1 <- da[1:146]
t2 <- da[147:333]
au1 <- data.frame(loai=t1,author=1)
au2 <- data.frame(loai=t2,author=2)
au <- rbind(au1,au2)
library(DT)
datatable(au)Tiếp theo ta vẽ đồ thị cột thể hiện điều này như sau:
ggplot(data=au,aes(x=loai,y=after_stat(count)))+geom_bar(fill="brown")+
geom_text(aes(label = scales::percent(after_stat(count/sum(count)))), stat = 'count', color = 'red', vjust = - .5) +
facet_grid(.~author) 
Ta thấy theo giả thuyết người thứ nhất chỉ nghiên cứu loài Adelie chiểm 43.8%, người thứ 2 nghiên cứu loài Chinstrap (20.4%) và Gentoo(35.7%)
3 Tuần 4
3.1 Scatter plot
Trong phần này tôi thực hiện vẽ biểu đồ trực quan hóa đồ thị phân tán cho dataset penguins: Trước khi vẽ, ta thực hiện tập hợp các lệnh xử lí dữ liệu như sau:
library('palmerpenguins') #gọi package 'palmerpenguins' xuất hiện
h <- penguins # gán penguins vào h
h <- na.omit(h) #loại bỏ các giá trị NA
names(h) <- c("s","i","bl","bd","f","bm","sx","y") #đặt tên cho các biến có trong dữ liệu
library(DT)
datatable(h)library(utf8)
library(tidyverse)Tiếp theo, ta nghi ngờ rằng 2 biến bill_length_mm (Độ dài mỏ chim) và flipper_length_mm (độ dài cánh) có tương quan với nhau, ta thực hiện vẽ scatter plot như sau:
library(ggplot2)
h%>%ggplot(aes(x=bl,y=f))+geom_point(color='red',size=3)+
xlab("Độ dài mỏ")+
ylab("Độ dài cánh")
Nhìn qua đồ thị ta có thể nghi ngờ rằng 2 biến này có tương quan thuận
với nhau. Ta làm rõ hơn điều này bằng cách xây dựng mô hình hồi quy 2
biến như sau:
model1 <- lm(f~bl,data=h)
summary(model1)##
## Call:
## lm(formula = f ~ bl, data = h)
##
## Residuals:
## Min 1Q Median 3Q Max
## -43.413 -7.837 0.652 8.360 21.321
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 127.3304 4.7291 26.93 <2e-16 ***
## bl 1.6738 0.1067 15.69 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 10.63 on 331 degrees of freedom
## Multiple R-squared: 0.4265, Adjusted R-squared: 0.4248
## F-statistic: 246.2 on 1 and 331 DF, p-value: < 2.2e-16Ta thấy hệ số tương quan bình phương R-squared = 0.4265, do đó 2 biến bill_length_mm và flipper_length_mm có tương quan thuận với nhau. Mô hình hồ quy f = 127.3304 + 1.6738bl. Vẽ thêm đường hồi quy cho 2 biến này, ta được kết quả như sau:
h%>%ggplot(aes(x=bl,y=f))+geom_point(color='red',size=3)+
geom_smooth(formula=y~x,method='lm',color='blue')+
xlab('Độ dài mỏ')+
ylab('Độ dài cánh')
Đây chính là đường hồi quy của hàm f = 127.3304 + 1.6738bl. Ý nghĩa hệ số góc: Nếu độ dài mỏ chim cánh cụt tăng lên 1mm thì độ dài cánh chim cánh cụt tăng 1.6738mm
Ta vẽ đồ thị phân tán giữa độ dài mỏ chim độ dài cánh và tô màu theo từng loài, như sau:
h%>%ggplot(aes(x=bl,y=f))+geom_point(aes(color=s,na.rm=T))+
geom_smooth(aes(color=s),formula=y~x,method='lm',na.rm = T)+
facet_wrap(~s)+
labs(x="Độ dài mỏ",y="Độ dài cánh")## Warning in geom_point(aes(color = s, na.rm = T)): Ignoring unknown aesthetics:
## na.rm
Giờ đây ta có thể thấy rằng: loài Adelie có độ dài mỏ chim dưới 45mm và độ dài cánh hầu hết dưới 200mm, do đó độ dài cánh và độ dài mỏ của loài Adelie nhỏ nhất trong 3 loài. Độ dài mỏ của của 2 loài Chinstrap và Gentoo tương đương nhau từ 40mm-60mm nhưng độ dài cánh của chúng là khác nhau. Cánh của loài Gentoo từ 210 -230mm, trong khi đó cánh của loài Chinstap từ 185mm-210mm. Do đó cùng số đo về độ dài mỏ nhưng cánh của loài Gentoo thì lớn hơn cánh của loài Chinstrap.
h%>%ggplot(aes(x=bl,y=f))+geom_point(aes(color=i,na.rm=T))+
geom_smooth(aes(color=i),formula=y~x,method='lm',na.rm = T)+
facet_wrap(~i)+
xlab('Độ dài mỏ')+
ylab('Độ dài cánh')## Warning in geom_point(aes(color = i, na.rm = T)): Ignoring unknown aesthetics:
## na.rm
Qua biểu đồ ta có thể thấy độ dốc của đường hồi quy của cá thể loài thuộc đảo Biscoe dốc nhất do đó cá thể chim cánh cụt ở đây có độ dài mỏ càng tăng thì sự tăng trưởng về độ dài cánh càng mạnh và hơn hẳn so với 2 đảo còn lại. Ngoài ra, cá thể cánh cụt có mỏ dài hơn 40mm và cánh dài hơn 210mm thì thường sinh sống ở đảo Biscoe. Cá thể cánh cụt có mỏ dài từ 30mm-55mm nhưng và cánh dài từ 170mm-210mm thì thường sống ở đảo Dream. Cuối cùng, đảo Torgersen là nơi cư trú của cá thể chim cánh cụt mỏ từ 0-45mm và cánh từ 170mm-200mm.
Tiếp theo, ta chỉ quan tâm tới những cá thể cánh cụt có độ dài cánh trên 210mm và tô màu xanh nước biến cho những cá thể này. Câu lệnh như sau:
h %>% ggplot(aes(x=bl,y=f,shape=i))+
geom_point(color="red",size=3)+
geom_point(data = h %>% filter(f > 210 ), size = 3, color = 'blue') +
xlab('Độ dài mỏ')+
ylab('Độ dài cánh')
Đầu tiên tôi vẽ biểu đồ phân tán thể hiện mối liên hệ giữa độ dài mỏ và độ dài cánh theo nơi sinh sống và tô màu đỏ cho toàn bộ quan sát. Sau đó tôi thêm một layer mới vẫn là vẽ đồ thị đó nhưng lần này ta tô màu cho những quan sát có độ dài cánh lớn hơn 210mm là màu xanh nước biển. Kết quả ta có thể thấy những chấm tròn màu xanh nước biển chính là những cá thể có độ dài cánh lớn hơn 210mm. Và ta có thể thấy toàn bộ cá thể có cánh trên 210mm đều sống ở đảo Biscoe.
Tiếp theo ta chỉ quan tâm tới các cá thể chim cánh cụt giống đực. Ta tiếp tục vẽ đồ thị phân tán giữa 2 biến độ dài mỏ và chiều dài cánh và chia theo từng loài
h %>% ggplot(aes(x=bl,y=f,shape=s))+
geom_point(color='purple',size=3)+
geom_point(data = h %>% filter(sx=="male"), size = 3, color = 'green') +
xlab('Độ dài mỏ')+
ylab('Độ dài cánh')
Đầu tiên ta thực hiện vẽ đồ thị phân tán và tô màu cho toàn bộ cá thể là màu tím. Sau đó tôi thêm 1 layer mới vẫn là vẽ đồ thị đó nhưng lần này tôi tô màu cho cá thể đực là màu xanh lá. Ta thấy rằng những chấm màu xanh lá thể hiện những cá thể cánh cụt giống đực và những chấm màu tím thể hiện cá thể cánh cụt giống cái. Nhìn vào biểu đồ ta thấy cá thể đực phân bồ ở tất cả mọi nơi.
3.2 Biểu đồ cột
Ta vẽ biểu đồ cột cho biến nơi sinh sống của loài chim cánh cụt như sau:
h |> ggplot(aes(x=i))+geom_bar(fill="yellow")
Ta thấy rằng có hơn 150 cá thể chim cánh cụt sống ở Biscoe, khoảng 120 cá thể sống ở Dream và gần 50 cá thế sống ở ToTorgersen
Ta không muốn dùng số đếm, ta dùng tỉ lệ phần trăm như sau:
library(scales)##
## Attaching package: 'scales'## The following objects are masked from 'package:psych':
##
## alpha, rescale## The following object is masked from 'package:purrr':
##
## discard## The following object is masked from 'package:readr':
##
## col_factorh |> group_by(i) |>
summarise(n = n()) |>
mutate(pr = percent(n/sum(n), accuracy = ,01)) |>
ggplot(aes(x = i, y = pr)) +
geom_col(fill = 'blue') +
theme_classic() +
labs(x = 'Noi sinh song', y = 'Ty le %')
Dựa vào biểu đồ ta có thể thấy 49% cá thể chim cánh cụt sống ở Biscoe, 37% sinh sống ở đảo Dream và 14% sinh sống ở đảo Torgersen.
Tiếp theo, tôi tiến hành vẽ biểu đồ sự phân bố các cá thể chim cánh cụt ở từng hòn đảo, và phân chia biểu đồ theo biến Species (loài) để phân tích.
h |> ggplot(aes(x = i, y = after_stat(count))) +
geom_bar(fill = 'purple') +
geom_text(aes(label = scales::percent(after_stat(count/sum(count)))), stat = 'count', color = 'red', vjust = - .5) +
facet_grid(.~s) +
labs(x = 'Noi sinh song', y = 'So luong')
Dựa vào biểu đố trên ta có thể thấy loài Adelie phân bố hầu hết ở cả 3 hòn đảo: trong đó cá thể sống ở Biscoe chiếm 13.21%, cá thể sống ở Dream chiếm 16.52% và ở Torgersen chiếm 14.11%. Loài chinstrap chỉ phân bố ở đảo Dream chiểm 20.42% và loài Gentoo chỉ sổng ở đảo Biscoe chiếm 35.74%.
Tiếp theo tôi thực hiện vẽ biểu đồ thể hiện sự phân bố cá thể chim cánh cụt theo giới tính, chia biểu đồ theo từng loài. Câu lệnh như sau:
table(h$s,h$sx)##
## female male
## Adelie 73 73
## Chinstrap 34 34
## Gentoo 58 61h |> ggplot(aes(x = sx, y = after_stat(count))) +
geom_bar(fill = 'green') +
geom_text(aes(label = scales::percent(after_stat(count/sum(count)))), stat = 'count', color = 'red', vjust = - .5) +
facet_grid(.~s) +
# theme_classic() +
labs(x = 'Giới tính', y = 'Số lượng')
Ta có thể thấy, loài Adelie các nhà nghiên cứu lấy mẫu giống đực và cái với tỷ lệ bằng nhau, giống đực chiếm 21.92% và giống cái chiếm 21.92%. Loài Chinstrap tương tự, giống đực chiếm 10.21% và giống cái chiếm 10.21%. Riêng loài Gentoo, giống cái ít hơn giống đực (17.42%<18.32%)
Tiếp theo, ta thêm một biến mới có tên là size. Biến này được tạo ra với mục đích là để phân tổ không đều cho biến bill_length_mm (độ dài mỏ). Khoảng từ (32-39.5] được xem là mỏ ngắn, khoảng từ (39.5,48.6] được xem là vừa và khoảng từ (48.6,59.6] được xem là mỏ dài
h$size <- cut(h$bl,breaks = c(32,39.5,48.6,59.6),lables=c("ngắn","vừa","dài"))
table(h$size)##
## (32,39.5] (39.5,48.6] (48.6,59.6]
## 86 164 83Ta thấy có 86 cá thể mỏ ngắn, 164 cá thể mỏ vừa và 83 cá thể mỏ dài. Trực quan hóa điều này trên đồ thị và phân chia theo loài, ta được kết quả như sau:
h %>% ggplot(aes(x=size,y=after_stat(count)))+geom_bar(fill="blue")+
geom_text(aes(label = scales::percent(after_stat(count/sum(count)))), stat = 'count', color = 'red', vjust = - .5) +
facet_wrap(~s)+
labs(x="Độ dài mỏ",y="số lượng")
Ta thấy ở loài Adelie, cá thể có độ dài mỏ ngắn chiếm 25.8%, cá thể mỏ vừa chiếm 18%. Ở loài Chinstrap, không có cá thể nào mỏ ngắn, cá thể mỏ vừa chiếm 8.7% và cá thể mỏ lớn chiếm 11.7%. Ở loài Gentoo cũng không có cá thể nào mỏ ngắn, có 22.5% cá thể mỏ vừa và 13.2% cá thể mỏ dài.
Tiếp theo, ta thực hiện phân tích biến năm thu thập dữ liệu
h |> count(y) |>
mutate(pC = percent(n/sum(n),accuracy = 0.01)) |>
ggplot(aes(x = y, y = n)) +
geom_col(fill = 'blue') +
geom_text(aes(label = pC),color = 'yellow', vjust = 2, size = 5) +
ylab("Số cá thể") +
xlab("Năm thu thập số liệu")
Ta có thể thấy năm 2007 các nhà khoa học thu thập được 30.93% cá thể, năm 2008 thực hiện nghiên cứu 33.93% cá thể và năm 2009 nghiên cứu nhiều nhất với 35.14% trên tổng số. Chia theo từng khu vực sinh sống ta được kết quả như sau:
table(h$y,h$i)##
## Biscoe Dream Torgersen
## 2007 43 45 15
## 2008 63 34 16
## 2009 57 44 16Ta thấy năm 2007 nghiên cứu 43 cá thể ở đảo Biscoe, 45 cá thể ở đảo Dream và 15 cá thể ở đảo Torgersen. Năm 2008 nghiên cứu 63 cá thể ở đảo Biscoe, 34 cá thể ở đảo Dream và 15 cá thể ở đảo Torgersen. Năm 2009 nghiên cứu 57 cá thể ở đảo Biscoe, 44 cá thể ở đảo Dream và 16 cá thể ở đảo Torgersen. Vẽ điều này lên biểu đồ với tỷ lệ phần trăm như sau:
h %>% ggplot(aes(x=y,y=after_stat(count)))+geom_bar(fill="blue")+
geom_text(aes(label = scales::percent(after_stat(count/sum(count)))), stat = 'count', color = 'red', vjust = - .5) +
facet_wrap(~i)
labs(x="Năm thu thập số liệu",y="số lượng")## $x
## [1] "Năm thu thập số liệu"
##
## $y
## [1] "số lượng"
##
## attr(,"class")
## [1] "labels"Ta thấy rằng các cá thể được nghiên cứu nhiều nhất ở đảo Biscoe, cá thể được nghiên cứu ít nhất ở đảo Torgersen qua từng năm.
h %>% count(y) %>%
mutate(Pc=percent(n/sum(n),accuracy=1))## # A tibble: 3 × 3
## y n Pc
## <int> <int> <chr>
## 1 2007 103 31%
## 2 2008 113 34%
## 3 2009 117 35%3.3 Ma trận tương quan
3.3.1 Đồ thị histogram, đồ thị phân tán và hệ số tương quan.
Trong phần này tôi thực hiện tính ma trận tương quan giữa các biến.
Ngoài ra tôi còn vẽ đồ thị thể hiện phân phối giữa các biến và đồ thị
phân tán. Các điều này được thể hiện trên duy nhất 1 hình ảnh. Trước
hết, muốn vẽ ma trận tương quan giữa các biến có trong bộ dữ liệu
penguins, ta thực hiện library package
PerformanceAnalytics. Sau đó ta thực hiện vẽ ma trận tương
quan, các câu lệnh như sau:
library("PerformanceAnalytics")## Loading required package: xts## Loading required package: zoo##
## Attaching package: 'zoo'## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numeric##
## ######################### Warning from 'xts' package ##########################
## # #
## # The dplyr lag() function breaks how base R's lag() function is supposed to #
## # work, which breaks lag(my_xts). Calls to lag(my_xts) that you type or #
## # source() into this session won't work correctly. #
## # #
## # Use stats::lag() to make sure you're not using dplyr::lag(), or you can add #
## # conflictRules('dplyr', exclude = 'lag') to your .Rprofile to stop #
## # dplyr from breaking base R's lag() function. #
## # #
## # Code in packages is not affected. It's protected by R's namespace mechanism #
## # Set `options(xts.warn_dplyr_breaks_lag = FALSE)` to suppress this warning. #
## # #
## #################################################################################
## Attaching package: 'xts'## The following objects are masked from 'package:dplyr':
##
## first, last##
## Attaching package: 'PerformanceAnalytics'## The following objects are masked from 'package:fBasics':
##
## kurtosis, skewness## The following object is masked from 'package:graphics':
##
## legendtuongquan <- h[,c(3:6)]
chart.Correlation(tuongquan, histogram=TRUE,pch=19)## Warning in par(usr): argument 1 does not name a graphical parameter
## Warning in par(usr): argument 1 does not name a graphical parameter
## Warning in par(usr): argument 1 does not name a graphical parameter
## Warning in par(usr): argument 1 does not name a graphical parameter
## Warning in par(usr): argument 1 does not name a graphical parameter
## Warning in par(usr): argument 1 does not name a graphical parameter
Đầu tiên ta có thể thấy bên trái đường chéo là đồ thị phân tán giữa các biến. Đường chéo chính thể hiện đồ thị phân phối histogram của các biếnbiến
Bên phải đường chéo chính là hệ số tương quan. Dựa vào ma trận tương quan ta có thể thấy các biến có tương quan nghich với nhau bao gồm: tương quan khá yếu giữa độ dài mỏ(bl) và độ sâu mỏ(bd) với hệ số tương quan là -0.23, độ sâu mỏ(bd) và độ dài cánh(f) với hệ số tương quan -0.58,độ sâu mỏ (bd) và khối lương cơ thể(bm) với hệ số tương quan là -0.47. Tiếp theo, các biến có tương quan thuận với nhau bao gồm: độ dài mỏ (bl) và độ dài cánh(f) với hệ số tương quan là 0.65, độ dài mỏ (bl) và khối lượng cơ thể(bm) với hệ số tương quan là 0.59, độ dài cánh và(f) khối lượng cơ thể(bm) với hệ số tương quan là 0.87. Trong đó 2 biến f và bm có tương quan mạnh nhât, hai biến bl và bd có tương quan yếu nhất.
3.3.2 Đồ thị corplot
library("PerformanceAnalytics")
prs <- h[c(3:6)]
res <- cor(prs)
library(corrplot)## corrplot 0.92 loadedcorrplot(res, type = "upper", order = "hclust",
tl.col = "black", tl.srt = 45)
Ta thấy rằng có 3 ô màu cam thể hiện hệ số tương quan âm. Tức các cặpcặp biến độ dài mỏ (bl) và độ sâu mỏ (bd), độ dài cánh (f) và độ sâu mỏ (bd), độ sâu mỏ (bd) và khối lượng cơ thể (bm) đều tương quan nghịch với nhau. Ngoài ra, giữa 2 biến f và bd có màu cam đậm thể hiện 2 biến này tương quan nghịch mạnh hơn các biến còn lại. Những ô tròn màu xanh còn lại thể hiện sự tương quan thuận lẫn nhau. Giữa 2 biến f và bm có màu xanh đậm thể hiện rằng 2 biến này tương quan thuận chặt chẽ với nhau.
3.4 Vẽ bản đồ các quốc gia trên thế giới
3.4.1 Dataset “World”
Dataset world là bộ dữ liệu về địa điểm tất cả mọi nơi trên trái đất, mỗi một nơi sẽ có kinh độ và vỹ độ riêng biệt. Bao gồm các biến như sau:
long: kinh độ (longitude) của các điểm trên đường biên quốc gia
lat: vĩ độ (latitude) của các điểm trên đường biên quốc gia
group: chỉ số của các group (nhóm) được sử dụng để vẽ các đường biên quốc gia. Có thể có nhiều group cho cùng một quốc gia nếu có các hòn đảo hoặc vùng lãnh thổ của quốc gia.
order: chỉ số sắp xếp các điểm trên mỗi group theo thứ tự vẽ
region: tên khu vực địa lý của quốc gia (Asia, Europe, …)
subregion: tên nhóm phân vùng địa lý của quốc gia (Eastern Asia, Southern Europe, …)
Bây giờ ta thực hiện lệnh gọi dataset world.map như sau:
world.map <- map_data("world")3.4.2 Trực quan hóa hình dạng của các quốc gia
Bây giờ, ta thực hiện trực quan hóa hình dáng của tất cá các nước trên thế giới như sau:
ggplot(data=world.map, aes(x=long, y=lat, group=group, fill=factor(group))) + geom_polygon(col="white") + theme(legend.position="none")
Câu lệnh trên tôi sử dụng nhằm vẽ bản đồ của tất cả các nước trên thế giới trên một mặt phẳng nghiêng.
Tiếp theo tôi thực hiện vẽ biểu đồ của các quốc gia trên thế giới như sau:
countries <- c("Russia","Ukraine") #gọi quốc gia cần vẽ
map <- map_data("world",region=countries)
ggplot(data=map, aes(x=long, y=lat, group=group, fill=factor(group))) + geom_polygon(col="white")+ theme(legend.position="none")
Ta thấy rằng Nga là quốc gia được tô màu xanh, Ukraine là quốc gia được tô màu đỏ. Diện tích của Nga lớn hơn rất nhiều so với Ukraine. Nhìn vào vị trí 2 quốc gia ta có thể hiểu lý do vì sao xung đột giữa Nga - Ukraine diễn ra trong bối cảnh Ukraine cho phép Mỹ đặt căn cứ quân sự tại quốc gia của mình.
Tiếp theo ta vẽ bản đồ của 3 nước Đông Dương như sau:
countries <- c("Vietnam","Laos","Cambodia") #gọi quốc gia cần vẽ
map <- map_data("world",region=countries)
ggplot(data=map, aes(x=long, y=lat, group=group, fill=factor(group))) + geom_polygon(col="white")+ theme(legend.position="none")
Quốc gia được tô màu hồng là Việt Nam, tô màu xanh là Laos và màu vàng là Campuchia. Nhìn hình vẽ ta có thể đoán được rằng diện tích của Việt Nam lớn hơn so với Laos và Campuchia. Diện tích của Laos lớn hơn diện tích của Campuchia.Vị trí của cả 3 nước sát nhau.
Vẽ bản đồ của các nước Đông Nam Á như sau:
countries <- c("Vietnam","Laos","Cambodia","Thailand","Myanmar","Indonesia","Malaysia","Philippines","Singapore","Brunei") #gọi quốc gia cần vẽ
map <- map_data("world",region=countries)
ggplot(data=map, aes(x=long, y=lat, group=group, fill=factor(group))) + geom_polygon(col="white")+ theme(legend.position="none")
Việc vẽ bản đồ các nước theo các cách trên đây cho ra hình dáng chưa bao gồm các tỉnh, thành phố… cho nên ta thực hiện vẽ biểu đồ thể hiện rõ ranh giới giữa các tỉnh của Việt Nam. Trước khi vẽ, ta cần gọi các package hỗ trợ như sau:
library(ggplot2)
library(raster)## Loading required package: sp## The legacy packages maptools, rgdal, and rgeos, underpinning this package
## will retire shortly. Please refer to R-spatial evolution reports on
## https://r-spatial.org/r/2023/05/15/evolution4.html for details.
## This package is now running under evolution status 0##
## Attaching package: 'raster'## The following object is masked from 'package:dplyr':
##
## selectlibrary(sp)
library(tidyverse)
library(viridis)## Loading required package: viridisLite##
## Attaching package: 'viridis'## The following object is masked from 'package:scales':
##
## viridis_pallibrary(rgdal)## Please note that rgdal will be retired during October 2023,
## plan transition to sf/stars/terra functions using GDAL and PROJ
## at your earliest convenience.
## See https://r-spatial.org/r/2023/05/15/evolution4.html and https://github.com/r-spatial/evolution
## rgdal: version: 1.6-7, (SVN revision 1203)
## Geospatial Data Abstraction Library extensions to R successfully loaded
## Loaded GDAL runtime: GDAL 3.5.2, released 2022/09/02
## Path to GDAL shared files: C:/Users/ADMIN/AppData/Local/R/win-library/4.2/rgdal/gdal
## GDAL binary built with GEOS: TRUE
## Loaded PROJ runtime: Rel. 8.2.1, January 1st, 2022, [PJ_VERSION: 821]
## Path to PROJ shared files: C:/Users/ADMIN/AppData/Local/R/win-library/4.2/rgdal/proj
## PROJ CDN enabled: FALSE
## Linking to sp version:1.6-1
## To mute warnings of possible GDAL/OSR exportToProj4() degradation,
## use options("rgdal_show_exportToProj4_warnings"="none") before loading sp or rgdal.##
## Attaching package: 'rgdal'## The following object is masked from 'package:fBasics':
##
## getDescriptionTa thực hiện vẽ biểu đồ Việt Nam bao gồm ranh giới các tỉnh như sau:
vnm = getData("GADM", country="Vietnam", level=1)## Warning in getData("GADM", country = "Vietnam", level = 1): getData will be removed in a future version of raster
## . Please use the geodata package insteadvn = fortify(vnm, regions="VARNAME_1")## Regions defined for each Polygonsggplot(data=vn,aes(x=long,y=lat,group=group))+geom_polygon(aes(fill=id))+theme(legend.position="none")
4 Tuần 3
4.1 Xoay trục dữ liệu sau đó so sánh 2 bộ dữ liệu
4.1.1 Chọn dataset có cấu trúc dữ liệu theo hàng và xoay trục bằng pivot_longer()
0 Bảng dữ liệu về an toàn hàng không (airline safety) là một tập dữ liệu liên quan đến các thông tin về các chuyến bay trong một số hãng hàng không trên thế giới được tích hợp trong package “fivethirtyeight”. Mỗi hàng là 1 quan sát tương ứng với mỗi hãng hàng không, các cột bao gồm các biến mô tả đặc điểm của hãng hàng không đó, nên dữ liệu này được tổ chức theo hàng. Bảng dữ liệu này thường được sử dụng để phân tích các chỉ tiêu về an toàn của hãng hàng không, đánh giá mức độ rủi ro trong các chuyến bay và đối chiếu với các tiêu chuẩn an toàn cơ bản của ngành hàng không.Các biến có trong dữ liệu bao gồm:
airline: Tên của hãng hàng không
incl_reg_subsidiaries: hãng hàng không có công ty con hay không?
avail_seat_km_per_week: Số ghế được phục vụ bởi hãng hàng không trong một tuần
incidents_00_14: Số vụ tai nạn được báo cáo từ năm 2000 đến 2014
fatal_accidents_00_14: Số vụ tai nạn được báo cáo từ năm 2000 đến 2014 với ít nhất một nạn nhân chết
fatalities_00_14: Tổng số lượng hành khách thiệt mạng trong các vụ tai nạn được báo cáo từ năm 2000 đến 2014
incidents_85_99: Số vụ tai nạn được báo cáo từ năm 1985 đến 1999
fatal_accidents_85_99: Số vụ tai nạn được báo cáo từ năm 1985 đến 1999 với ít nhất một nạn nhân chết
fatalities_85_99”: Tổng số lượng hành khách thiệt mạng trong các vụ tai nạn/khủng bố được báo cáo từ năm 1985 đến 1999
library(fivethirtyeight)## Some larger datasets need to be installed separately, like senators and
## house_district_forecast. To install these, we recommend you install the
## fivethirtyeightdata package by running:
## install.packages('fivethirtyeightdata', repos =
## 'https://fivethirtyeightdata.github.io/drat/', type = 'source')library(tidyverse)
data("airline_safety")
str(airline_safety)## tibble [56 × 9] (S3: tbl_df/tbl/data.frame)
## $ airline : chr [1:56] "Aer Lingus" "Aeroflot" "Aerolineas Argentinas" "Aeromexico" ...
## $ incl_reg_subsidiaries : logi [1:56] FALSE TRUE FALSE TRUE FALSE FALSE ...
## $ avail_seat_km_per_week: num [1:56] 3.21e+08 1.20e+09 3.86e+08 5.97e+08 1.87e+09 ...
## $ incidents_85_99 : int [1:56] 2 76 6 3 2 14 2 3 5 7 ...
## $ fatal_accidents_85_99 : int [1:56] 0 14 0 1 0 4 1 0 0 2 ...
## $ fatalities_85_99 : int [1:56] 0 128 0 64 0 79 329 0 0 50 ...
## $ incidents_00_14 : int [1:56] 0 6 1 5 2 6 4 5 5 4 ...
## $ fatal_accidents_00_14 : int [1:56] 0 1 0 0 0 2 1 1 1 0 ...
## $ fatalities_00_14 : int [1:56] 0 88 0 0 0 337 158 7 88 0 ...names(airline_safety)=c("a","irs","askp","i85","fa85","f85","i00","fa00","f00")
asl <- airline_safety %>%
pivot_longer(cols = c(i85,i00,fa85,fa00),
names_to = "ad",
values_to = "number")Đầu tiên tôi gọi package fivethirtyeight nhằm trích xuất
bộ dữ liệu airline_safety có trong package này. Sau đó đặt
tên cho các biến bằng lệnh names(). Sau đó tôi thực hiện
lệnh pivot_longer() cho các biến liên quan đến “số vụ tai
nạn” bao gồm i85 là số vụ tai nạn được báo cáo
từ năm 1985 đến 1999,i00 là Số vụ
tai nạn được báo cáo từ năm 2000 đến 2014, “fa85” là
số vụ tai nạn được báo cáo từ năm 1985 đến 1999 với ít nhất
một nạn nhân chết ,0fa00 là số vụ
tai nạn được báo cáo từ năm 2000 đến 2014 với ít nhất một nạn nhân
chết với mục đích chuyển hết các biến này từ hàng về dạng
cột. Lệnh names_to="ad" nhằm đặt tên cho cột mới tạo chứa
tên của 4 biến này là “ad” và lệnh values_to="number" nhằm
đặt tên cho cột chứa giá trị của các biến này là “number”. Toàn bộ quá
trình tôi thực hiện gán vào biến asl.
Kết quả thu được tệp dữ liệu asl với 224 hàng và 7 cột,
trên thực tế bộ dữ liệu này dài hơn dữ liệu gốc gấp 4
lần vì mỗi hàng trong bảng dữ liệu ban đầu đại diện cho 4
hàng trong asl và ít cột hơn. Mỗi hàng cho 1 quan sát số vụ
tai nạn theo 4 thời điểm khác nhau. Hai cột mới tạo ra là “ad” có 4 biểu
hiện tương ứng với 4 thời điểm và “number” chứa giá trị của các
biến.
table(asl$ad)##
## fa00 fa85 i00 i85
## 56 56 56 56Biến ad có 4 biểu hiện thỏa mãn yêu cầu.
as <- airline_safety
sum <- as$f00+as$f85
aa <- data.frame(as,sum)
ff <- aa %>%
pivot_longer(cols = c(f85,f00,sum),
names_to = "dn",
values_to = "count")Tiếp theo tôi thêm biến sum là tổng của 2 biến f0 và
f85 nhằm tính tổng số người tử vong của cả 2 giai đoạn. Mục
đích là để so sánh dữ liệu cũ và dữ liệu sau khi xoay trục liệu rằng dữ
liệu nào sẽ cho kết quả trực quan tốt hơn trong phần sau 1.1.2. Dữ liệu
mới là aa gồm 10 biến và 56 quan sát.
Tương tự tôi tiếp thực hiện lệnh pivot_longer cho các
biến liên quan đến “số người tử vong” nhằm chuyển các biến này từ hàng
về cột. Các biến bao gồm: f85 là Tổng số lượng
hành khách thiệt mạng trong các vụ tai nạn/khủng bố được báo cáo từ năm
1985 đến 1999, f00 là Tổng số
lượng hành khách thiệt mạng trong các vụ tai nạn được báo cáo từ năm
2000 đến 2014 và biến sum tổng số
người tử vong của cả hai giai đoạn này . Lệnh
names_to="dn" nhằm đặt tên cho cột mới chứa 3 biến này là
“dn” và lệnh value_to="count" nhằm đặt tên cho cột chứa giá
trị của 2 biến này là “count”. Toàn bộ quá trình được gán vào biến
ff.
Kết quả thu được tệp dữ liệu ff gồm 9 biến và 168 quan
sát. Trên thực tế bộ dữ liệu này dài hơn bộ dữ liệu gốc 3 lần và ít biến
hơn (9<10). Trong đó 2 biến mới là dn gồm 3 biểu hiện là
“f00” và f85” và “sum”, biến count chứa các giá trị tương
ứng của 3 biểu hiện này.
table(ff$dn)##
## f00 f85 sum
## 56 56 56Biến dn cps 3 biểu hiện thỏa mãn yêu cầu.
4.1.2 So sánh dữ liệu gốc so với dữ liệu đã xoay trục
Sau khi thực hiện xoay trục dữ liệu thành cột dài hơn, tôi có 2 bộ dữ
liệu mới là asl và
ff. . Tiếp theo, tôi thực hiện việc so sánh
xem liệu rằng với dữ liệu gốc airline_safety và dữ liệu mới
ff thì dữ liệu nào sẽ dễ phân tích hơn thông qua việc trực
quan hóa nó
Tôi thực hiện trực quan hóa biểu đồ cột bao gồm 4 biến về số vụ tai nạn bao gồm đối với dữ liệu gốc như sau:
library(ggplot2)
ggplot(data=aa) +
geom_col(aes(x =irs, y = sum), width = 1)
Biểu đồ trên thể hiện tổng số trường hợp tử vong trong cả 2 giai đoạn
theo biến irs (hãng hàng không có công ty con
hay không). Rõ ràng với bộ dữ liệu ban đầu, chúng ta chỉ
có thể nhận xét rằng những hãng hàng không lớn có công ty con hoạt động
thì số ca tử vong do tai nạn máy bay sẽ ít hơn những hãng hàng không
không có công ty con. Nhìn vào dữ liệu ff được xoay trục
dưới đây, chúng ta sẽ thấy điều khác biệt:
ff %>%
filter(dn != "sum") %>%
ggplot() +
geom_col(
aes(x = irs, y = count, fill = dn),
width = 1
)
Biểu đồ của tệp dữ liệu ff đã trở nên khác biệt và đẹp
mắt hơn dữ liệu gốc. Nhìn vào kết quả trực quan này, ta không chỉ nhận
ra rằng những hãng hàng không lớn có công ty con hoạt động thì số ca tử
vong do tai nạn máy bay sẽ ít hơn những hãng hàng không không có công ty
con mà chúng ta còn so sánh được số ca tử vong do tai nạn máy bay qua
từng giai đoạn. Biểu đồ cột chồng trên cho ta thấy số ca tử vong trong
giai đoạn 1985 đến 1999 xấp xỉ gấp đôi số người tử vong trong giai đoạn
2000-2014 trong cả 2 biểu hiện của biến irs.
4.1.3 Phân tích
sự phụ thuộc giữa 2 biến i85 và
askp
plot(i85 ~ askp, data = airline_safety,col="green")
Câu lệnh trên dùng để kiểm tra tính tuyến tính bằng biểu đồ phân tán
(scatter plot) giữa hai biến avail_seat_km_per_week và
incidents_85_99 để xem xét mối tương quan giữa chúng. Nhìn
hình ta chưa đoán được rằng giữa số ghế ngồi của hãng và số vụ tai nạn
có tương quan với nhau hay không. Do đó ta thực hiện thêm bài toán kiểm
định sự tương quan giữa 2 biến như sau:
Bài toán kiểm định
- H0: Không có tương quan (R=0)
- H1: Có sự tương quan (R khác 0)
Thực hiện như sau:
cor.test(airline_safety$i85,airline_safety$askp)##
## Pearson's product-moment correlation
##
## data: airline_safety$i85 and airline_safety$askp
## t = 2.1395, df = 54, p-value = 0.03693
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.01795744 0.50530365
## sample estimates:
## cor
## 0.2795383Ta thấy p_value = 0.03693 < 5% bác bỏ giả thiết
H0, nghĩa là 2 biến này có tương quan với nhau.
4.1.4 Dùng hàm pivot_wider chuyển dữ liệu về lại trạng thái ban đầu
head(asl)## # A tibble: 6 × 7
## a irs askp f85 f00 ad number
## <chr> <lgl> <dbl> <int> <int> <chr> <int>
## 1 Aer Lingus FALSE 320906734 0 0 i85 2
## 2 Aer Lingus FALSE 320906734 0 0 i00 0
## 3 Aer Lingus FALSE 320906734 0 0 fa85 0
## 4 Aer Lingus FALSE 320906734 0 0 fa00 0
## 5 Aeroflot TRUE 1197672318 128 88 i85 76
## 6 Aeroflot TRUE 1197672318 128 88 i00 6library(tidyverse)
hung <- asl %>% pivot_wider(names_from=ad,
values_from=number)
head(hung)## # A tibble: 6 × 9
## a irs askp f85 f00 i85 i00 fa85 fa00
## <chr> <lgl> <dbl> <int> <int> <int> <int> <int> <int>
## 1 Aer Lingus FALSE 320906734 0 0 2 0 0 0
## 2 Aeroflot TRUE 1197672318 128 88 76 6 14 1
## 3 Aerolineas Argentinas FALSE 385803648 0 0 6 1 0 0
## 4 Aeromexico TRUE 596871813 64 0 3 5 1 0
## 5 Air Canada FALSE 1865253802 0 0 2 2 0 0
## 6 Air France FALSE 3004002661 79 337 14 6 4 2head(airline_safety)## # A tibble: 6 × 9
## a irs askp i85 fa85 f85 i00 fa00 f00
## <chr> <lgl> <dbl> <int> <int> <int> <int> <int> <int>
## 1 Aer Lingus FALSE 320906734 2 0 0 0 0 0
## 2 Aeroflot TRUE 1197672318 76 14 128 6 1 88
## 3 Aerolineas Argentinas FALSE 385803648 6 0 0 1 0 0
## 4 Aeromexico TRUE 596871813 3 1 64 5 0 0
## 5 Air Canada FALSE 1865253802 2 0 0 2 0 0
## 6 Air France FALSE 3004002661 14 4 79 6 2 337Kết quả tạo ra được datteset hung giống với dataset ban
đầu airline_safety với cấu trúc dữ liệu theo hàng.
4.2 Phân tích dữ liệu với dataset penguins
4.2.1 Tạo bảng thống kê chi tiết
Gọi lại dataset penguins (các biến đã được giải thích ở
tuần trước)
library(palmerpenguins)
data(penguins)
penguins <- na.omit(penguins)
names(penguins) <-c("s","i","bl","bd","f","bm","sx","y")Thực hiện thống kê mô tả về độ dài cánh cho chuỗi dữ liệu
penguins theo từng loài một cách chi tiết và đổi tên các
đặc trưng đo lường như sau:
library(tidyverse)
df1 <- penguins %>% rename(Loai=s) %>%
group_by(Loai)%>%
summarise(TB=mean(f),PS=var(f),LC=sd(f),NN=mean(f),LN=max(f))
df1## # A tibble: 3 × 6
## Loai TB PS LC NN LN
## <fct> <dbl> <dbl> <dbl> <dbl> <int>
## 1 Adelie 190. 42.5 6.52 190. 210
## 2 Chinstrap 196. 50.9 7.13 196. 212
## 3 Gentoo 217. 43.4 6.59 217. 231Tới đây ta thu được bảng dữ liệu về đặc trưng đo lường nhưng còn xấu. Tôi thực hiện làm đẹp bảng dữ liệu hơn và có thể thực hiện tìm kiếm nếu số hàng của bảng giá trị lớn, như sau:
library(DT)
datatable(df1,rownames=FALSE,colnames=c("Loài", "Trung bình", "Phương sai", "Độ lệch chuẩn", "Nhỏ nhất","Lớn nhất"),
caption="Bảng 1: Kết quả khảo sát các loài chim cánh cụt")4.2.2 Hồi quy 2 biến
Tôi nghi ngờ rằng 2 biến độ dài cánh
flipper_length_mm và khối lượng cơ thể
body_mass_g có sự tương quan với nhau. Do đó tôi thực hiện
vẽ đồ thị phân tán, như sau:
plot(bm~f,data=penguins,col="blue")
Nhìn qua đồ thị ta có thể nghi ngờ rằng 2 biến này có tương quan thuận với nhau. Ta làm rõ hơn điều này bằng cách xây dựng mô hình hồi quy 2 biến như sau:
model1 <- lm(bm~f,data=penguins)
summary(model1)##
## Call:
## lm(formula = bm ~ f, data = penguins)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1057.33 -259.79 -12.24 242.97 1293.89
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -5872.09 310.29 -18.93 <2e-16 ***
## f 50.15 1.54 32.56 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 393.3 on 331 degrees of freedom
## Multiple R-squared: 0.7621, Adjusted R-squared: 0.7614
## F-statistic: 1060 on 1 and 331 DF, p-value: < 2.2e-16Ta thấy rằng hệ số tương quan R-square=0.76210, cho thấy 2 biến này có tương quan khá mạnh. Kết quả hồi quy như sau bm= - 5872 + 50.15f. Ý nghĩa hệ số góc: Nếu độ dài cánh tăng lên 1mm thì khối lượng cơ thể chim cánh cụt tăng 50g.
4.2.3 Kiểm định sự tương quan
Sau khi có kết quả hồi quy ở trên, ta thấy hàm hồi quy có ý nghĩa, để chắc chắn hơn ta kiểm định sự tương quan giữa 2 biến như sau
Bài toán kiểm định
- H0: Không có tương quan (R=0)
- H1: Có sự tương quan (R khác 0)
Thực hiện như sau:
cor.test(penguins$f,penguins$bm)##
## Pearson's product-moment correlation
##
## data: penguins$f and penguins$bm
## t = 32.562, df = 331, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.8447622 0.8963550
## sample estimates:
## cor
## 0.8729789Ta thấy p_value < 2.2e-16 < 5%, bác bỏ H0. Vậy 2
biến body_mass_g và flipper_length_mm có tương
quan tuyến tính với nhau.
4.2.4 Kiểm định phân phối chuẩn
Để kiểm định một chuỗi dữ liệu có phải là phân phối chuẩn hay không, chúng ta có rất nhiều phương pháp. Trong phần này, tôi sử dụng kiểm định Jarque-Bera để kiểm tra điều đó
Ta có bài toán kiểm định
- H0: Chuỗi dữ liệu có phân phối chuẩn
- H1: Chuỗi dữ liệu không có phân phối chuẩn
Thực hiện kiểm định cho biến độ dài cánh chim cánh cụt:
library(tseries)
jarque.bera.test(penguins$f)##
## Jarque Bera Test
##
## data: penguins$f
## X-squared = 20.05, df = 2, p-value = 4.427e-05Dựa vào bảng trên ta thấy kết quả kiểm định là 20.05 và
p_value= 4.427*10^-5 < 5% bác bỏ H0, vậy biến
flipper_length_mm không có phân phối chuẩn
Thực hiện kiểm định cho biến khối lượng cơ thể chim cánh cụt:
jarque.bera.test(penguins$bm)##
## Jarque Bera Test
##
## data: penguins$bm
## X-squared = 19.874, df = 2, p-value = 4.835e-05p_value=4.835*10^-5 < 5%, bác bỏ H0. Do đó biến
body_mass_g không có phân phối chuẩn
4.2.5 Vẽ biểu đồ
dung <- table(penguins$s)
dung##
## Adelie Chinstrap Gentoo
## 146 68 119Bảng tần số trên thể hiện số chim cánh cụt theo từng loài, có 146 cá thể thuộc loài Adelie, 68 cá thể thuôc loài Chinstrap và 119 cá thể thuộc loài GentoGentoo. Ta thực hiện vẽ biểu đồ cột cho biến này như sau:
library(fBasics)
mausac <- qualiPalette(12,"Set3")
barplot(dung,col=mausac,main="phân bố chim cánh cụt theo loài")
Kết quả cho thấy loài Adelie chiếm số lượng lớn nhất, loàn Chinstrap
chiếm số lượng nhỏ nhất trong bộ dữ liệu.
dang <- table(penguins$sx,penguins$s)
dang##
## Adelie Chinstrap Gentoo
## female 73 34 58
## male 73 34 61Tên đây là bảng tần số về tỉ lệ giới tính theo các vùng Tiếp theo ta thực hiện vẽ biểu đồ trực quan điểu này. Thực hiện lấy dải màu thuộc “Set3” bằng lệnh qualiPalette().
mausac <- qualiPalette(12,"Set3")
barplot(dang,col=mausac,main="Tỉ lệ giới tính theo từng loài",
beside=TRUE,horiz=TRUE)
Qua biểu đồ ta thấy bộ số liệu nghiên cứu có tỉ lệ cá thể thể đực và cái khá giống nhau ở cả 3 loài.
5 Tuần 2
5.1 Giải thích dữ liệu:
Dữ liệu penguins được cung cấp bởi Tiến sĩ Kristen Gorman và Trạm Palmer, Nam Cực. Dữ liệu thu thập các đặc điểm của các loài chim cánh cụt sống ở Nam Cực. Dữ liệu này được tích hợp trong package “palmerpenguins” bao gồm 8 biến:
species: các loại chim cánh cụt
island : hòn đảo chim cánh cụt sinh sống
bill_length_mm: chiều dài mỏ chim cánh cụt
bill_depth_mm: chiều rộng mỏ chim cánh cụt
body_mass_g: khối lượng cơ thể chim cánh cụt
sex: giống (gồm đực và cái)
flipper_length_mm: độ dài cánh của chim cánh cụt
year: năm thu thập số liệu
5.2 Lọc và phân tổ không đều dữ liệu
Trước hết, gọi package chứa dữ liệu penguins xuất hiện bằng lệnh library()
library(palmerpenguins)5.2.1 Lọc dữ liệu
names(penguins) <-c("s","i","bl","bd","f","bm","sx","y")
penguins <- na.omit(penguins) Đầu tiên, tôi thực hiện đặt tên cho các biến có trong dữ liệu bằng lệnh names(), kết quả thu được biến s (species), biến i (island), biến bl (bill_length_mm), biến bd (bill_depth_mm), biến bm (body_mass_g), biến sx (sex), biến y (year) . Sau đó tôi loại bỏ các giá trị missing value cho tệp dữ liệu penguins bằng lệnh na.omit(penguins), kết quả loại bỏ đi được những cột và hàng chứa giá trị NA.
Gentoo <- penguins[penguins$s=="Gentoo",]
str(Gentoo)## tibble [119 × 8] (S3: tbl_df/tbl/data.frame)
## $ s : Factor w/ 3 levels "Adelie","Chinstrap",..: 3 3 3 3 3 3 3 3 3 3 ...
## $ i : Factor w/ 3 levels "Biscoe","Dream",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ bl: num [1:119] 46.1 50 48.7 50 47.6 46.5 45.4 46.7 43.3 46.8 ...
## $ bd: num [1:119] 13.2 16.3 14.1 15.2 14.5 13.5 14.6 15.3 13.4 15.4 ...
## $ f : int [1:119] 211 230 210 218 215 210 211 219 209 215 ...
## $ bm: int [1:119] 4500 5700 4450 5700 5400 4550 4800 5200 4400 5150 ...
## $ sx: Factor w/ 2 levels "female","male": 1 2 1 2 2 1 1 2 1 2 ...
## $ y : int [1:119] 2007 2007 2007 2007 2007 2007 2007 2007 2007 2007 ...
## - attr(*, "na.action")= 'omit' Named int [1:11] 4 9 10 11 12 48 179 219 257 269 ...
## ..- attr(*, "names")= chr [1:11] "4" "9" "10" "11" ...Tiếp theo tôi thực hiện việc lọc ra một bộ dữ liệu mới về đặc điểm của chim cánh cụt thuộc loài Gento. Kết quả thu được tệp dữ liệu mới (Gentoo) gồm 199 cá thể Gentoo và 8 biến.
notGentoo <- penguins[!penguins$s=="Gentoo",]
str(notGentoo)## tibble [214 × 8] (S3: tbl_df/tbl/data.frame)
## $ s : Factor w/ 3 levels "Adelie","Chinstrap",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ i : Factor w/ 3 levels "Biscoe","Dream",..: 3 3 3 3 3 3 3 3 3 3 ...
## $ bl: num [1:214] 39.1 39.5 40.3 36.7 39.3 38.9 39.2 41.1 38.6 34.6 ...
## $ bd: num [1:214] 18.7 17.4 18 19.3 20.6 17.8 19.6 17.6 21.2 21.1 ...
## $ f : int [1:214] 181 186 195 193 190 181 195 182 191 198 ...
## $ bm: int [1:214] 3750 3800 3250 3450 3650 3625 4675 3200 3800 4400 ...
## $ sx: Factor w/ 2 levels "female","male": 2 1 1 1 2 1 2 1 2 2 ...
## $ y : int [1:214] 2007 2007 2007 2007 2007 2007 2007 2007 2007 2007 ...
## - attr(*, "na.action")= 'omit' Named int [1:11] 4 9 10 11 12 48 179 219 257 269 ...
## ..- attr(*, "names")= chr [1:11] "4" "9" "10" "11" ...Tiếp theo tôi thực hiện việc lọc ra bộ dữ liệu mới về đặc điểm của chim cánh cụt không thuộc loài Gentoo, nghĩa là thuộc các loài còn lại. Kết quả thu được bộ dữ liệu mới (notGentoo) gồm 213 cá thể không thuộc loài Gentoo và 8 đặc điểm. Câu lệnh được thực hiện như trên
Chinstrap <- penguins[penguins$s=="Chinstrap",]
str(Chinstrap)## tibble [68 × 8] (S3: tbl_df/tbl/data.frame)
## $ s : Factor w/ 3 levels "Adelie","Chinstrap",..: 2 2 2 2 2 2 2 2 2 2 ...
## $ i : Factor w/ 3 levels "Biscoe","Dream",..: 2 2 2 2 2 2 2 2 2 2 ...
## $ bl: num [1:68] 46.5 50 51.3 45.4 52.7 45.2 46.1 51.3 46 51.3 ...
## $ bd: num [1:68] 17.9 19.5 19.2 18.7 19.8 17.8 18.2 18.2 18.9 19.9 ...
## $ f : int [1:68] 192 196 193 188 197 198 178 197 195 198 ...
## $ bm: int [1:68] 3500 3900 3650 3525 3725 3950 3250 3750 4150 3700 ...
## $ sx: Factor w/ 2 levels "female","male": 1 2 2 1 2 1 1 2 1 2 ...
## $ y : int [1:68] 2007 2007 2007 2007 2007 2007 2007 2007 2007 2007 ...
## - attr(*, "na.action")= 'omit' Named int [1:11] 4 9 10 11 12 48 179 219 257 269 ...
## ..- attr(*, "names")= chr [1:11] "4" "9" "10" "11" ...Tiếp theo tôi lọc ra các đặc điểm các loài chim thuộc loài Chinstrap. Với chuỗi các câu lệnh dưới đây, dữ liệu các loài Chinstrap được tạo ra với 68 hàng và 8 cột, tức là có 68 cá thể chim thuộc loài Chinstrap cùng với 8 đặc điểm. Câu lệnh được thực hiện như trên
thu2007 <- penguins[penguins$y==2007,]
str(thu2007)## tibble [103 × 8] (S3: tbl_df/tbl/data.frame)
## $ s : Factor w/ 3 levels "Adelie","Chinstrap",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ i : Factor w/ 3 levels "Biscoe","Dream",..: 3 3 3 3 3 3 3 3 3 3 ...
## $ bl: num [1:103] 39.1 39.5 40.3 36.7 39.3 38.9 39.2 41.1 38.6 34.6 ...
## $ bd: num [1:103] 18.7 17.4 18 19.3 20.6 17.8 19.6 17.6 21.2 21.1 ...
## $ f : int [1:103] 181 186 195 193 190 181 195 182 191 198 ...
## $ bm: int [1:103] 3750 3800 3250 3450 3650 3625 4675 3200 3800 4400 ...
## $ sx: Factor w/ 2 levels "female","male": 2 1 1 1 2 1 2 1 2 2 ...
## $ y : int [1:103] 2007 2007 2007 2007 2007 2007 2007 2007 2007 2007 ...
## - attr(*, "na.action")= 'omit' Named int [1:11] 4 9 10 11 12 48 179 219 257 269 ...
## ..- attr(*, "names")= chr [1:11] "4" "9" "10" "11" ...Tiếp theo tôi lọc ra các đặc điểm các cá thể chim cánh cụt được nghiên cứu vào năm 2007. Kết quả thu được tệp dữ liệu thu2007 với 103 hàng cà 8 cột, tức là có 103 cá thể được nghiên cứu vào năm 2007.
duc2008 <- penguins[penguins$sx=="male"&penguins$y==2008,]
str(duc2008)## tibble [57 × 8] (S3: tbl_df/tbl/data.frame)
## $ s : Factor w/ 3 levels "Adelie","Chinstrap",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ i : Factor w/ 3 levels "Biscoe","Dream",..: 1 1 1 1 1 1 1 1 1 3 ...
## $ bl: num [1:57] 40.1 42 41.4 40.6 37.6 41.3 41.1 41.6 41.1 41.8 ...
## $ bd: num [1:57] 18.9 19.5 18.6 18.8 19.1 21.1 18.2 18 19.1 19.4 ...
## $ f : int [1:57] 188 200 191 193 194 195 192 192 188 198 ...
## $ bm: int [1:57] 4300 4050 3700 3800 3750 4400 4050 3950 4100 4450 ...
## $ sx: Factor w/ 2 levels "female","male": 2 2 2 2 2 2 2 2 2 2 ...
## $ y : int [1:57] 2008 2008 2008 2008 2008 2008 2008 2008 2008 2008 ...
## - attr(*, "na.action")= 'omit' Named int [1:11] 4 9 10 11 12 48 179 219 257 269 ...
## ..- attr(*, "names")= chr [1:11] "4" "9" "10" "11" ...Tiếp theo tôi lọc ra các loài chim cánh cụt giống đực và được nghiên cứu vào năm 2008. Kết quả thu được dữ liệu duc2008 với 57 hàng và 8 cột thể hiện 57 cá thể chim cánh cụt đực và 8 đặc điểm của chúng ứng với 8 biến
female190to <- penguins[penguins$f<190&penguins$sx=="female"&penguins$i=="Torgersen",]
str(female190to)## tibble [13 × 8] (S3: tbl_df/tbl/data.frame)
## $ s : Factor w/ 3 levels "Adelie","Chinstrap",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ i : Factor w/ 3 levels "Biscoe","Dream",..: 3 3 3 3 3 3 3 3 3 3 ...
## $ bl: num [1:13] 39.5 38.9 41.1 36.6 34.4 36.2 34.6 36.7 38.6 35.7 ...
## $ bd: num [1:13] 17.4 17.8 17.6 17.8 18.4 16.1 17.2 18.8 17 17 ...
## $ f : int [1:13] 186 181 182 185 184 187 189 187 188 189 ...
## $ bm: int [1:13] 3800 3625 3200 3700 3325 3550 3200 3800 2900 3350 ...
## $ sx: Factor w/ 2 levels "female","male": 1 1 1 1 1 1 1 1 1 1 ...
## $ y : int [1:13] 2007 2007 2007 2007 2007 2008 2008 2008 2009 2009 ...
## - attr(*, "na.action")= 'omit' Named int [1:11] 4 9 10 11 12 48 179 219 257 269 ...
## ..- attr(*, "names")= chr [1:11] "4" "9" "10" "11" ...Tiếp theo tôi lọc ra các cá thể chim cánh cụt có chiều dài cánh dưới 190mm, giống cái, và cư trú ở đảo Torgesen. Kết quả thu được dữ liệu female190to gồm 13 cá thể thỏa mãn yêu cầu và 8 cột thể hiện 8 đặc điểm của chúng.
5.2.2 Phân tổ không đều và lập bảng tần số, tần suất cho các tổ
penguins$size <- cut(penguins$bl,breaks = c(32,39.5,48.6,59.6),lables=c("ngắn","vừa","dài"))
x <- table(penguins$size)
x##
## (32,39.5] (39.5,48.6] (48.6,59.6]
## 86 164 83Trước tiên tôi thêm một biến mới có tên là size. Biến này chứa dữ liệu liên tục về chiều dài mỏ chim cánh cụt cho biết cá thể nào cỏ mỏ “dài”, cá thể nào có mỏ “ngắn” và cá thể nào “vừa” và gán dữ liệu vào biến x
Khoảng từ (32,39.5] là được xem là “ngắn”.
Khoảng từ (39.5,48.6] là “vừa”.
Khoảng từ (48.6,59.6] là “dài”.
Tiếp theo tôi thực hiện lập bảng tần số cho biến mới này. Kết quả có 86 con cá thể mỏ ngắn, 164 cá thể mỏ vừa và 83 cá thể mỏ dài.
prop.table(x)##
## (32,39.5] (39.5,48.6] (48.6,59.6]
## 0.2582583 0.4924925 0.2492492Lập bảng tần suất cho biến x. Kết quả: 25.83% chim cánh cụt mỏ ngắn, 49.25% chim cánh cụt mỏ vừa và 24.92% chim cánh cụt mỏ dài.
summary(penguins$bd)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 13.10 15.60 17.30 17.16 18.70 21.50penguins$bdsize <- cut(penguins$bd,breaks = c(13.09,15.6,18.7,21.5),lables=c("nhỏ","vừa","lớn"))
y <- table(penguins$bdsize)
y##
## (13.1,15.6] (15.6,18.7] (18.7,21.5]
## 85 169 79Ngoài ra tôi còn thêm một biến mới với tên là bdsize. Biến chứa dữ liệu liên tục về chiều rộng mỏ chim cánh cụt cho biết cá thể nào mỏ “nhỏ”, cá thể nào mỏ “vừa” và cá thể nào mỏ “lớn” sau đó lập bảng tần số cho các khoảng và gán vào biến y
Khoảng từ (13.09,15.6] được xem là mỏ “nhỏ”
Khoảng từ (15.6,18.7] được xem là mỏ “vừa”
Khoảng từ (18.7,21.5] được xem là mỏ “lớn”
Kết quả thu được 85 cá thể mỏ “nhỏ”, 169 cá thể mỏ “vừa” và 79 cá thể mỏ “lớn”.
prop.table(y)##
## (13.1,15.6] (15.6,18.7] (18.7,21.5]
## 0.2552553 0.5075075 0.2372372Lập bảng tần suất cho biến y vừa tìm. Kết quả: 25.53% cá thể mỏ nhỏ, 50.75% cá thể mỏ vừa và 23.72% cá thể mỏ lớn
5.3 Thống kê mô tả các đặc trưng đo lường
5.3.1 Thống kê mô tả tổng hợp các đặc trưng với hàm summary()
summary(penguins$f)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 172 190 197 201 213 231Tôi thực hiện thống kê mô tả các đặc trưng đo lường của các biến. Đầu tiên với lệnh summary(penguins$f), tôi thu được các đặc trưng của biến độ dài cánh của chim cánh cụt như sau:
Độ dài cánh nhỏ nhất là 172mm
Tứ phân vị thứ nhất là 190 cho biết có 25% cá thể chim cánh cụt có cánh dài dưới 190mm và 75% trên 190mm.
Trung vị là 197 cho biết có 50% cá thể chim cánh cụt có cánh dài trên 197mm và 50% dưới 197mm.
Trung bình là 201 cho biết độ dài cánh trung bình là 201mm.
Tứ phân vị thứ 3 là 213 cho biết có 75% cá thể chim cánh cụt có cánh dài dưới 213mm và 25% trên 213mm.
5.3.2 Thống kê mô tả riêng lẻ các đặc trưng đo lường
Bây giờ tôi tiếp tực thực hiện việc thống kê mô tả các biến nhưng chỉ tập trung vào 1 đặc trưng đo lường. Câu lênh mean(penguins$f) cho ra giá trị trung bình của biến độ dài cánh là 200.967mm.
mean(penguins$f)## [1] 200.967Câu lệnh sd(penguins$bl) cho biết chiều dài mỏ chim cánh cụt có độ lệch chuẩn là 5.468668mm tức là chênh lệch giữa độ dài mỏ chim được quan sát so với giá trị trung bình độ dài mỏ chim là 5.468668mm.
sd(penguins$bl)## [1] 5.468668Câu lệnh var(penguins$bm) chi biết độ lệnh chuẩn bình phương của khối lượng cơ thể chim cánh cụt được quan sát so với khối lượng trung bình của chim cánh cụt là 648372.5mm
var(penguins$bm)## [1] 648372.5Câu lệnh quantile(penguins$bm,0.6) cho biết có 60% cá thể chim cánh cụt có khối lượng dưới 4300g và 40% trên 4300g
quantile(penguins$bm,0.6)## 60%
## 43005.3.3 Thống kê mô tả với hàm aggregate và toán tử đường ống “%>%”
Tiếp theo tôi tiến hành tính trung bình của độ dài cánh của chim cánh cụt theo từng loài. Kết quả độ dài cánh trung bình loài Adelie nhỏ nhất xấp xỉ 190.103mm, đứng thứ 2 là loài Chinstrap với 195.8235mm và lớn nhất là loài Gentoo với 217.2353mm. Câu lệnh thực hiện như sau:
aggregate(penguins$f,list(penguins$s),FUN="mean")## Group.1 x
## 1 Adelie 190.1027
## 2 Chinstrap 195.8235
## 3 Gentoo 217.2353Tiếp theo tôi tiến hành tính độ lệch chuẩn của độ dài cánh của chim cánh cụt theo từng loài. Kết quả thu được:
Độ dài cánh loài Adelie chênh lệch xấp xỉ 6.522mm so với độ dài cánh trung bình.
Độ dài cánh loài Chinstrap chênh lệch xấp xỉ 7.132mm so với trung bình.
Độ dài cánh loài Gentoo chênh lệch xấp xỉ 6.585mm so với trung bình. Câu lệnh như sau:
aggregate(penguins$f,list(penguins$s),FUN="sd")## Group.1 x
## 1 Adelie 6.521825
## 2 Chinstrap 7.131894
## 3 Gentoo 6.585431Tiếp theo, tôi thực hiện lập các toán tử đường ống để tính toán nhanh hơn các đặc trưng đo lường. Trước hết ta phải gọi gói tidyverse xuất hiện để thực hiện việc tính toán này. Cụ thể tôi tính trung bình của độ dài mỏ chim thoe từng loài. Kết quả như sau:
Độ dài mỏ loài Adelie có trung bình là 38.82mm
Độ dài mỏ loài Chinstrap trung bình là 48.83mm
Độ dài mỏ loài Gentoo trung bình là 47.5mm
library(tidyverse)
penguins %>% group_by(s) %>% summarise(m=mean(bl))## # A tibble: 3 × 2
## s m
## <fct> <dbl>
## 1 Adelie 38.8
## 2 Chinstrap 48.8
## 3 Gentoo 47.6Tiếp theo, giả định rằng dữ liệu penguins được 2 người thu thập với tổng cộng 333 cột(không tính các cột missing value). Độ dài cánh . Tôi sẽ tạo ra một dataframe có giá trị là 1 và 2 tương ứng dữ liệu với người thứ nhất và người thứ hai thu thập.
Dữ liệu từ cột 1 đến cột 150 do người thứ nhất thu thập, được mã hóa bởi số 1
Dữ liệu từ cột 151 đến cột 333 do người thứ hai thu thập, được mã hóa bởi số 2
Quy trình thực hiện như sau:
Đầu tiên gọi biến s và gán vào stt (mỗi một phần tử trong biến s tương ứng với 1 cá thể chim cánh cụt)
Tiếp theo lấy 150 phần tử đầu tiên (do người thứ nhất thu thập) gán vào biến t1, 183 phần từ tiếp theo (do người thứ 2 thu thập) gán vào biến t2.
tiếp theo tạo một data frame tg1 gồm 2 biến chim (các phần tử của biến s từ 1-150) và nguoithuchien (1)
Tương tự tạo data frame tg2 gồm 2 biến chim(các phần tử của biến s từ 151-333) và nguoithuchien (2)
Gộp 2 dataframe tg1 tg2 lại bằng lệnh rbind() và gán vào tg. Kết quả tạo ra dataframe tg thỏa mãn yêu cầu.
stt <- penguins$s
t1 <- stt[1:150]
t2 <- stt[151:333]
tg1 <- data.frame(chim=t1,nguoithuchien=1)
tg2 <- data.frame(chim=t2,nguoithuchien=2)
tg <- rbind(tg1,tg2)
str(tg)## 'data.frame': 333 obs. of 2 variables:
## $ chim : Factor w/ 3 levels "Adelie","Chinstrap",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ nguoithuchien: num 1 1 1 1 1 1 1 1 1 1 ...5.4 Trực quan hóa dữ liệu
Trước hết để vẽ đồ thị ta cần gọi gói ggplot2 xuất hiện. Sau đó tôi thực hiện vẽ biểu đồ cột bằng lệnh geom_bar() nhằm đánh giá số lượng cá thể chim cánh cụt sinh sống ở mỗi hòn đảo khác nhau. Kết quả thu được đảo Biscoe là nơi cư trú nhiều nhất nhất của chim cánh cụt, đảo Torgersen có số lượng chim cánh cụt cư trú ít nhất.
library(ggplot2)
ggplot(data=penguins,aes(x=i,fill=i))+geom_bar()
6 Tuần 1
6.1 Giải thích dữ liệu:
Dữ liệu penguins được cung cấp bởi Tiến sĩ Kristen Gorman và Trạm Palmer, Nam Cực. Dữ liệu thu thập các đặc điểm của các loài chim cánh cụt sống ở Nam Cực. Dữ liệu này được tích hợp trong package “palmerpenguins” bao gồm 8 biến:
species: các loại chim cánh cụt
island : hòn đảo chim cánh cụt sinh sống
bill_length_mm: chiều dài mỏ chim cánh cụt
bill_depth_mm: chiều rộng mỏ chim cánh cụt
body_mass_g: khối lượng cơ thể chim cánh cụt
sex: giống (gồm đực và cái)
flipper_length_mm: độ dài cánh của chim cánh cụt
year: năm thu thập số liệu
6.2 Thao tác trên dataset:
6.2.1 Xử lý dữ liệu
Đầu tiên, tôi gọi package palmerpenguins xuất hiện bằng hàm library() và chọn bộ dữ liệu “penguins” được tích hợp trong package này.
Tiếp theo, tôi làm sạch dữ liệu bằng câu lệnh ‘na.omit()’ vầ gán bộ dữ liệu vào biến h, kết quả loại bỏ đi những giá trị missing value.
Ngoài ra, tôi thực hiện đặt tên cho các biến có trong bộ dữ liệu penguins bằng câu lệnh names().
Tôi thực hiện đọc 6 dòng đầu tiên và 6 dòng cuối cùng của dữ liệu lần lượt với lệnh “head()” và “tail()”. Các câu lệnh được viết như sau:
library('palmerpenguins') #gọi package 'palmerpenguins' xuất hiện
str(penguins) # cấu trúc của dataset penguins## tibble [333 × 10] (S3: tbl_df/tbl/data.frame)
## $ s : Factor w/ 3 levels "Adelie","Chinstrap",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ i : Factor w/ 3 levels "Biscoe","Dream",..: 3 3 3 3 3 3 3 3 3 3 ...
## $ bl : num [1:333] 39.1 39.5 40.3 36.7 39.3 38.9 39.2 41.1 38.6 34.6 ...
## $ bd : num [1:333] 18.7 17.4 18 19.3 20.6 17.8 19.6 17.6 21.2 21.1 ...
## $ f : int [1:333] 181 186 195 193 190 181 195 182 191 198 ...
## $ bm : int [1:333] 3750 3800 3250 3450 3650 3625 4675 3200 3800 4400 ...
## $ sx : Factor w/ 2 levels "female","male": 2 1 1 1 2 1 2 1 2 2 ...
## $ y : int [1:333] 2007 2007 2007 2007 2007 2007 2007 2007 2007 2007 ...
## $ size : Factor w/ 3 levels "(32,39.5]","(39.5,48.6]",..: 1 1 2 1 1 1 1 2 1 1 ...
## $ bdsize: Factor w/ 3 levels "(13.1,15.6]",..: 2 2 2 3 3 2 3 2 3 3 ...
## - attr(*, "na.action")= 'omit' Named int [1:11] 4 9 10 11 12 48 179 219 257 269 ...
## ..- attr(*, "names")= chr [1:11] "4" "9" "10" "11" ...h <- penguins # gán penguins vào h
h <- na.omit(h) #loại bỏ các giá trị NA
str(h)## tibble [333 × 10] (S3: tbl_df/tbl/data.frame)
## $ s : Factor w/ 3 levels "Adelie","Chinstrap",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ i : Factor w/ 3 levels "Biscoe","Dream",..: 3 3 3 3 3 3 3 3 3 3 ...
## $ bl : num [1:333] 39.1 39.5 40.3 36.7 39.3 38.9 39.2 41.1 38.6 34.6 ...
## $ bd : num [1:333] 18.7 17.4 18 19.3 20.6 17.8 19.6 17.6 21.2 21.1 ...
## $ f : int [1:333] 181 186 195 193 190 181 195 182 191 198 ...
## $ bm : int [1:333] 3750 3800 3250 3450 3650 3625 4675 3200 3800 4400 ...
## $ sx : Factor w/ 2 levels "female","male": 2 1 1 1 2 1 2 1 2 2 ...
## $ y : int [1:333] 2007 2007 2007 2007 2007 2007 2007 2007 2007 2007 ...
## $ size : Factor w/ 3 levels "(32,39.5]","(39.5,48.6]",..: 1 1 2 1 1 1 1 2 1 1 ...
## $ bdsize: Factor w/ 3 levels "(13.1,15.6]",..: 2 2 2 3 3 2 3 2 3 3 ...
## - attr(*, "na.action")= 'omit' Named int [1:11] 4 9 10 11 12 48 179 219 257 269 ...
## ..- attr(*, "names")= chr [1:11] "4" "9" "10" "11" ...names(h) <- c("s","i","bl","bd","f","bm","sx","y") #đặt tên cho các biến có trong dữ liệu ## Warning: The `value` argument of `names<-` must have the same length as `x` as of tibble
## 3.0.0.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.## Warning: The `value` argument of `names<-` can't be empty as of tibble 3.0.0.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.head(h) #đọc 6 dòng đầu tiên của dữ liệu ## # A tibble: 6 × 10
## s i bl bd f bm sx y `` ``
## <fct> <fct> <dbl> <dbl> <int> <int> <fct> <int> <fct> <fct>
## 1 Adelie Torgersen 39.1 18.7 181 3750 male 2007 (32,39.5] (15.6,18.7]
## 2 Adelie Torgersen 39.5 17.4 186 3800 female 2007 (32,39.5] (15.6,18.7]
## 3 Adelie Torgersen 40.3 18 195 3250 female 2007 (39.5,48.6] (15.6,18.7]
## 4 Adelie Torgersen 36.7 19.3 193 3450 female 2007 (32,39.5] (18.7,21.5]
## 5 Adelie Torgersen 39.3 20.6 190 3650 male 2007 (32,39.5] (18.7,21.5]
## 6 Adelie Torgersen 38.9 17.8 181 3625 female 2007 (32,39.5] (15.6,18.7]tail(h) #đọc 6 dòng cuối cùng của dữ liệu## # A tibble: 6 × 10
## s i bl bd f bm sx y `` ``
## <fct> <fct> <dbl> <dbl> <int> <int> <fct> <int> <fct> <fct>
## 1 Chinstrap Dream 45.7 17 195 3650 female 2009 (39.5,48.6] (15.6,18.7]
## 2 Chinstrap Dream 55.8 19.8 207 4000 male 2009 (48.6,59.6] (18.7,21.5]
## 3 Chinstrap Dream 43.5 18.1 202 3400 female 2009 (39.5,48.6] (15.6,18.7]
## 4 Chinstrap Dream 49.6 18.2 193 3775 male 2009 (48.6,59.6] (15.6,18.7]
## 5 Chinstrap Dream 50.8 19 210 4100 male 2009 (48.6,59.6] (18.7,21.5]
## 6 Chinstrap Dream 50.2 18.7 198 3775 female 2009 (48.6,59.6] (15.6,18.7]6.2.2 Lọc dữ liệu
Trong phần này, tôi thực hiện chuỗi các lệnh lọc dữ liệu theo các trường hợp riêng biệt. Các câu lệnh được đánh số thứ tự bên dưới bằng cứ pháp #number nhằm thuận tiện cho việc giải thích. Câu lệnh #1 tôi thực hiện lọc ra bộ dữ liệu bao gồm đặc tính của những chú chim cánh cụt thuộc loài Adelie, lọc ra được 146 quan sát.
Adelie <-h[h$s=="Adelie",] #1
head(Adelie)## # A tibble: 6 × 10
## s i bl bd f bm sx y `` ``
## <fct> <fct> <dbl> <dbl> <int> <int> <fct> <int> <fct> <fct>
## 1 Adelie Torgersen 39.1 18.7 181 3750 male 2007 (32,39.5] (15.6,18.7]
## 2 Adelie Torgersen 39.5 17.4 186 3800 female 2007 (32,39.5] (15.6,18.7]
## 3 Adelie Torgersen 40.3 18 195 3250 female 2007 (39.5,48.6] (15.6,18.7]
## 4 Adelie Torgersen 36.7 19.3 193 3450 female 2007 (32,39.5] (18.7,21.5]
## 5 Adelie Torgersen 39.3 20.6 190 3650 male 2007 (32,39.5] (18.7,21.5]
## 6 Adelie Torgersen 38.9 17.8 181 3625 female 2007 (32,39.5] (15.6,18.7]Câu lệnh #2 nhằm lọc ra bộ dữ liệu mới bao gồm đặc điểm của những chú chim cánh cụt Adelie có độ dài cánh là 180mm, kết quả thu được 4 quan sát.
Adelie180 <-h[h$s=="Adelie"&h$f==180,] #2
head(Adelie180)## # A tibble: 4 × 10
## s i bl bd f bm sx y `` ``
## <fct> <fct> <dbl> <dbl> <int> <int> <fct> <int> <fct> <fct>
## 1 Adelie Biscoe 37.7 18.7 180 3600 male 2007 (32,39.5] (15.6,18.7]
## 2 Adelie Biscoe 38.8 17.2 180 3800 male 2007 (32,39.5] (15.6,18.7]
## 3 Adelie Biscoe 40.5 18.9 180 3950 male 2007 (39.5,48.6] (18.7,21.5]
## 4 Adelie Dream 42.2 18.5 180 3550 female 2007 (39.5,48.6] (15.6,18.7]Câu lệnh #3 tôi gán dữ liệu biến b_l (chiều dài mỏ chim cánh cụt) vào vol. Tiếp theo tôi lọc ra những con chim cánh cụt có chiều dài mỏ nhỏ hơn 50mm và gán vào biến vol_50 ở câu lệnh #4, kết quả thu được biến vol_50 bao gồm 276 giá trị.
vol <- h$bl #3
vol50 <- vol[vol<50] #4
str(vol50)## num [1:276] 39.1 39.5 40.3 36.7 39.3 38.9 39.2 41.1 38.6 34.6 ...Câu lệnh #4 tôi thực hiện lọc biến vol (chiều dài mỏ chim cánh cụt). Lọc ra những con chim cánh cụt có chiều dài mỏ từ 20-50mm và gán trực tiếp vào biến vol20_50. Kết quả lọc ra được 276 giá trị
vol20_50 <- vol[vol>20&vol<50] #lọc ra những con chim cánh cụt có chiều dài mỏ chim từ 20-50m #5
str(vol20_50)## num [1:276] 39.1 39.5 40.3 36.7 39.3 38.9 39.2 41.1 38.6 34.6 ...Tiếp theo tôi thực hiện lọc ra những con chim cánh cụt có chiều dài mỏ chim lớn hơn 30mm và chiều dài cánh nhỏ hơn 180mm và gán trực tiếp vào vol2bien ở câu lệnh #6
vol2bien <- h[h$bl>30&h$f<180,] #6
str(vol2bien)## tibble [7 × 10] (S3: tbl_df/tbl/data.frame)
## $ s : Factor w/ 3 levels "Adelie","Chinstrap",..: 1 1 1 1 1 1 2
## $ i : Factor w/ 3 levels "Biscoe","Dream",..: 1 1 2 2 2 3 2
## $ bl: num [1:7] 37.8 37.9 39.5 37.2 33.1 40.2 46.1
## $ bd: num [1:7] 18.3 18.6 16.7 18.1 16.1 17 18.2
## $ f : int [1:7] 174 172 178 178 178 176 178
## $ bm: int [1:7] 3400 3150 3250 3900 2900 3450 3250
## $ sx: Factor w/ 2 levels "female","male": 1 1 1 2 1 1 1
## $ y : int [1:7] 2007 2007 2007 2007 2008 2009 2007
## $ NA: Factor w/ 3 levels "(32,39.5]","(39.5,48.6]",..: 1 1 1 1 1 2 2
## $ NA: Factor w/ 3 levels "(13.1,15.6]",..: 2 2 2 2 2 2 2
## - attr(*, "na.action")= 'omit' Named int [1:11] 4 9 10 11 12 48 179 219 257 269 ...
## ..- attr(*, "names")= chr [1:11] "4" "9" "10" "11" ...Chuyển đổi giá trị của biến f thành log và gán vào biền log(f)
h$logf <- log(h$f)
head(h$logf)## [1] 5.198497 5.225747 5.273000 5.262690 5.247024 5.1984976.2.3 Lập bảng tần số
Trong phần này, tôi thực hiện lập bảng tần số cho một số biến có trong tập dữ liệu.
Câu lệnh #7 tôi phân tổ biến chiều dài mỏ chim cánh cụt (biến bl) thành 5 phần và gán vào biến vol5.
Sau đó lập bảng tần số cho 5 tổ ở câu lệnh #8. kết quả thu được
Tổ 1 (32.1,37.6] với 50 con chim.
Tổ 2 (37.6,43.1] với 100 con,.
Tổ 3 (43.1,48.6] gồm 100 con.
Tổ 4(48.6,54.1] với 76 con.
Tổ 5 (54.1,59.6] gồm 7 con.
vol5 <- cut(vol,5) #chia biến vol thành 5 phần #7
table(vol5) #lập bảng tần số 5 phần đó #8## vol5
## (32.1,37.6] (37.6,43.1] (43.1,48.6] (48.6,54.1] (54.1,59.6]
## 50 100 100 76 7Câu lệnh #9 tôi thực hiện phân tổ cho biến vol20_50 (đã được lọc ở trên) thành 4 tổ. Biến vol20_50 là tập hợp những con chim cánh cụt có chiều dài mỏ từ 20-50mm.
Câu lệnh #10 tôi lập bảng tần số cho các tổ. Kết quả thu được:
Tổ 1 có giá trị từ (32.1,36.5] với 33 cá thểthể,
Tổ 2 từ (36.5,41] với 83 cá thể,
Tổ 3 từ (41,45.5] với 67 cá thể,
Tổ 4 từ (45.5,49.9] với 93 cá thể.
d <- cut(vol20_50,4) #chia biến vol20_50 thành 4 phần #9
table(d) #lập bảng tần số 4 phần đó #10## d
## (32.1,36.5] (36.5,41] (41,45.5] (45.5,49.9]
## 33 83 67 93Câu lệnh #11 tôi phân tổ cho biến flipper_leng_mm (f) .Biến này có ý nghĩa là độ dài cánh của chim cánh cụt. Tôi phân thành 3 tổ và lập bảng tần số theo từng loài (biến species). Kết quả thu được
tổ 1 từ (172,192] vói 108 cá thể, trong đó loài Adelie chiếm 90 con, loài Chinstrap gồm 18 con, loại Gentoo 0 con.
Tổ 2 từ (192,211] gồm 131 cá thể, trong đó có 56 con thuộc loài Adelie, 49 con thuộc loài Chinstrap và 26 con thuộc loài Gentoo.
Tổ 3 từ (211,231] gồm 94 cá thể, trong đó không có con nào thuộc loài Adelie, 1 con thuộc loài Chinstrap và 93 con thuộc loài Gentoo.
table(cut(h$f,3),h$s) #11##
## Adelie Chinstrap Gentoo
## (172,192] 90 18 0
## (192,211] 56 49 26
## (211,231] 0 1 93Câu lệnh #12 tôi thực hiện phân tổ cho biến Bill_depth_mm (bd). Biến này có ý nghĩa là chiều rộng mỏ chim cánh cụt. Tôi phân biến này thành 3 tổ và lập bảng tần số theo năm thu thập số liệu (biến year). Kết quả thu được
tổ 1 có giá trị từ (13.1,15.9] gồm 99 phần tử, trong đó năm 2007 thu thập được 30 con thuộc khoảng này, năm 2008 thu thập được 37 con, năm 2009 thu thập được 32 con.
Tổ 2 từ (15.9,18.7] bao gồm 155 giá trị, trong đó ghi nhận được 43 con vào năm 2007, 50 con vào năm 2008 và 62 con vào năm 2009.
Tổ thứ 3 từ (18.7,21.5] các nhà nghiên cứu thu thập được số liệu của 79 con nằm trong khoảng này, cụ thể năm 2007 ghi nhận 30 con, năm 2008 nghiên cứu 26 con, năm 2009 nghiên cứu 23 con.
table(cut(h$bd,3),h$y) #12##
## 2007 2008 2009
## (13.1,15.9] 30 37 32
## (15.9,18.7] 43 50 62
## (18.7,21.5] 30 26 236.3 Trực quan hóa dữ liệu:
Trong phần này, tôi thực hiện trực quan hóa dư liệu giữa 2 biến: độ dài cánh (flipper_length_mm) và khối lượng cơ thể chim cánh cụt (body_mass_g) nhằm đánh giá mối liên hệ giữa 2 biến này, sử dụng package ’ggplot2” và vẽ đồ thị phân tán bằng hàm geom_point(). Hơn nữa, tôi còn chia đồ thị phân tán theo từng loài (species) khác nhau. Các kết quả thu được như sau:
Thứ nhất, dựa vào đồ thị phân tán ta có thể thấy cánh càng cài thì khối lượng cơ thể càng lớn.
Thứ hai, nhờ việc thực hiện phân chia đồ thị theo mỗi loài ta có thể thấy loài Gentoo là loài lớn nhất trong 3 loài.
library('ggplot2')
library('palmerpenguins')
d <- penguins
ggplot(data=penguins,aes(x=f,y=bm))+geom_point(aes(color=s,shape=s))+facet_wrap(~s)