Panel data with R

1.Call package:

pacman::p_load(rio,
               here,
               janitor,
               tidyverse,
               dplyr,
               magrittr,
               ggplot2,
               purrr,
               lubridate,
               mice,
               plotly)

2.Nhập file:

data(beav2,package ="MASS")

head(beav2)

##   day time  temp activ
## 1 307  930 36.58     0
## 2 307  940 36.73     0
## 3 307  950 36.93     0
## 4 307 1000 37.15     0
## 5 307 1010 37.23     0
## 6 307 1020 37.24     0

3.Chuyển object beav2 sang dạng time series:

#Check class của beav2 có phải dạng ts chưa
class(beav2)

## [1] "data.frame"

#Thêm cột hours để tạo ra dãy thời gian liên tục:
beav2<-beav2 %>% 
  mutate(hours = 24*(day-307)+time%/%100+(time%%100)/60)

#Tạo object TS:
temp<-ts(beav2$temp,start = 9.5,frequency = 6) 

4.Phân tích mô tả:

4.1.Xác định trending:

plot(temp)
abline(reg<-lm(temp~hours,data=beav2),col="blue")

Dữ liệu trên có dấu hiệu của xu hướng - trend và không có season. ## 4.2.So sánh 2 nhóm active 0 và 1:

temp1<-subset(beav2,activ==0)
temp2<-subset(beav2,activ==1)
#Dùng t.test để so sánh trung bình 2 nhóm:
t.test(temp1$temp,temp2$temp,var.equal = F)

## 
##  Welch Two Sample t-test
## 
## data:  temp1$temp and temp2$temp
## t = -18.548, df = 80.852, p-value < 2.2e-16
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -0.8927106 -0.7197342
## sample estimates:
## mean of x mean of y 
##  37.09684  37.90306

Kết quả cho thấy trung bình giữa 2 nhóm chênh lệch nhau khoảng 0.8 độ với chỉ số p có ý nghĩa. Điều đó cho thấy việc vận động (activ) làm tăng nhiệt độ cho cơ thể động vật (Một điều dễ nhận thấy)

5.Phân biệt các dạng mô hình phổ biến cho Time series data:

5.1.Simple Smoothing Average:

-Trung bình động (Simple Smoothing Average) là 1 số trung bình của 1 tập hợp dữ liệu được tính toán tích lũy. Thuật toán đơn giản nhằm xác định xu hướng tăng/giảm của chuỗi, có thể đưa ra quyết định nhanh và dễ giải thích. -Gỉa định: + sai số (error) không có dấu hiệu autocorrelation. + Trung bình = 0 và phương sai bất biến.

5.2.Mô hình ARIMA hay ARMA:

-ARIMA còn được gọi là các mô hình tự hồi quy trung bình trượt tích hợp. Mô hình sẽ tự hồi quy theo cơ chế MA (Moving Average) hoặc AR (Autoregression). Nhưng hiện nay, với sự phát triển của các kênh truyền thông, thông tin luôn được cập nhập liên tục nên các chuỗi thường được xây dựng trên cơ chế MA nhiều hơn AR. -Khác với SMA, ARIMA xây dựng với giả định rằng các sai số error có tương quan với nhau. -Gỉa định: + Chuỗi thời gian phải là chuỗi dừng.

6.Implementing model:

6.1.Sử dụng pp Simple Smoothing Average:

#Tạo object SMA:
temp_sma <- HoltWinters(temp, gamma = FALSE)
#Dự đoán giá trị của temp cho 6 period tiếp theo (1 period = 10 minutes)
temp_predict <- forecast:::forecast.HoltWinters(temp_sma, h = 6)

## Registered S3 method overwritten by 'quantmod':
##   method            from
##   as.zoo.data.frame zoo

#Plot graph:
plot(temp_predict)

print(temp_predict)

##          Point Forecast    Lo 80    Hi 80    Lo 95    Hi 95
## 26.16667       38.09779 37.91765 38.27793 37.82229 38.37329
## 26.33333       38.12558 37.85812 38.39304 37.71653 38.53463
## 26.50000       38.15337 37.81004 38.49670 37.62829 38.67845
## 26.66667       38.18116 37.76632 38.59600 37.54672 38.81560
## 26.83333       38.20895 37.72440 38.69351 37.46789 38.95001
## 27.00000       38.23674 37.68302 38.79046 37.38990 39.08358

6.2.Kiểm tra giả định:

#Có dấu hiệu autocorrelation không?
acf(temp_predict$residuals,lag.max = 10, na.action = na.pass) 

#Dùng Box test để kiếm tra:
Box.test(temp_predict$residuals, lag = 20, type = "Ljung-Box")

## 
##  Box-Ljung test
## 
## data:  temp_predict$residuals
## X-squared = 16.248, df = 20, p-value = 0.7011

#Error có phân phối chuấn không?
hist(temp_predict$residuals, main = "Phân phối của residuals")

#Phương sai có bất biến không?
plot(temp_predict$residuals)