原文

http://tecdat.cn/?p=3364?tecdat.cn

加載R包和數據集

上述癥狀數據集包含在R-package 中，并在加載時自動可用。 加載包后，我們將此數據集中包含的12個心情變量進行子集化：

mood_data <- as.matrix(symptom_data$data[, 1:12]) # Subset variables mood_labels <- symptom_data$colnames[1:12] # Subset variable labels colnames(mood_data) <- mood_labels time_data <- symptom_data$data_time

對象mood_data是一個1476×12矩陣，測量了12個心情變量：

> dim(mood_data) [1] 1476 12 > head(mood_data[,1:7]) Relaxed Down Irritated Satisfied Lonely Anxious Enthusiastic [1,] 5 -1 1 5 -1 -1 4 [2,] 4 0 3 3 0 0 3 [3,] 4 0 2 3 0 0 4 [4,] 4 0 1 4 0 0 4 [5,] 4 0 2 4 0 0 4 [6,] 5 0 1 4 0 0 3

time_data包含有關每次測量的時間戳的信息。數據預處理需要此信息。

> head(time_data) date dayno beepno beeptime resptime_s resptime_e time_norm 1 13/08/12 226 1 08:58 08:58:56 09:00:15 0.000000000 2 14/08/12 227 5 14:32 14:32:09 14:33:25 0.005164874 3 14/08/12 227 6 16:17 16:17:13 16:23:16 0.005470574 4 14/08/12 227 8 18:04 18:04:10 18:06:29 0.005782097 5 14/08/12 227 9 20:57 20:58:23 21:00:18 0.006285774 6 14/08/12 227 10 21:54 21:54:15 21:56:05 0.006451726

該數據集中的一些變量是高度偏斜的，這可能導致不可靠的參數估計。 在這里，我們通過計算自舉置信區間（KS方法）和可信區間（GAM方法）來處理這個問題，以判斷估計的可靠性。 由于本教程的重點是估計時變VAR模型，因此我們不會詳細研究變量的偏度。 然而，在實踐中，應該在擬合（時變）VAR模型之前始終檢查邊際分布。

估計時變VAR模型

通過參數lags = 1，我們指定擬合滯后1 VAR模型，并通過lambdaSel =“CV”選擇具有交叉驗證的參數λ。 最后，使用參數scale = TRUE，我們指定在模型擬合之前，所有變量都應縮放為零和標準差1。 當使用“1正則化”時，建議這樣做，因為否則參數懲罰的強度取決于預測變量的方差。 由于交叉驗證方案使用隨機抽取來定義折疊，因此我們設置種子以確保重現性。

在查看結果之前，我們檢查了1476個時間點中有多少用于估算，這在調用控制臺中的輸出對象時打印的摘要中顯示

> tvvar_obj mgm fit-object Model class: Time-varying mixed Vector Autoregressive (tv-mVAR) model Lags: 1 Rows included in VAR design matrix: 876 / 1475 ( 59.39 %) Nodes: 12 Estimation points: 20

估計的VAR系數的絕對值存儲在對象tvvar_obj $ wadj中，該對象是維度p×p×滯后×estpoints的數組。

參數估計的可靠性

res_obj <- resample(object = tvvar_obj, data = mood_data, nB = 50, blocks = 10,seeds = 1:50, quantiles = c(.05, .95))

res_obj $ bootParameters包含每個參數的經驗采樣分布。

計算時變預測誤差

函數predict（）計算給定mgm模型對象的預測和預測誤差。

預測存儲在pred_obj $預測中，并且所有時變模型的預測誤差組合在pred_obj中：

> pred_obj$errors Variable Error.RMSE Error.R2 1 Relaxed 0.939 0.155 2 Down 0.825 0.297 3 Irritated 0.942 0.119 4 Satisfied 0.879 0.201 5 Lonely 0.921 0.182 6 Anxious 0.950 0.086 7 Enthusiastic 0.922 0.169 8 Suspicious 0.818 0.247 9 Cheerful 0.889 0.200 10 Guilty 0.928 0.175 11 Doubt 0.871 0.268 12 Strong 0.896 0.195

可視化時變VAR模型

可視化上面估計的一部分隨時間變化的VAR參數：

# Two Network Plots # Get layout of mean graph Q <- qgraph(t(mean_wadj), DoNotPlot=TRUE) saveRDS(Q$layout, "Tutorials/files/layout_mgm.RDS") # Plot graph at selected fixed time points tpSelect <- c(2, 10, 18) # Switch to colorblind scheme tvvar_obj$edgecolor[, , , ][tvvar_obj$edgecolor[, , , ] == "darkgreen"] <- c("darkblue") lty_array <- array(1, dim=c(12, 12, 1, 20)) lty_array[tvvar_obj$edgecolor[, , , ] != "darkblue"] <- 2 for(tp in tpSelect) { qgraph(t(tvvar_obj$wadj[, , 1, tp]), layout = Q$layout, edge.color = t(tvvar_obj$edgecolor[, , 1, tp]), labels = mood_labels, vsize = 13, esize = 10, asize = 10, mar = rep(5, 4), minimum = 0, maximum = .5, lty = t(lty_array[, , 1, tp]), pie = pred_obj$tverrors[[tp]][, 3]) }

?

CIs <- apply(res_obj$bootParameters[par_row[1], par_row[2], 1, , ], 1, function(x) { quantile(x, probs = c(.05, .95)) } ) # Plot shading polygon(x = c(1:20, 20:1), y = c(CIs[1,], rev(CIs[2,])), col=alpha(colour = cols[i], alpha = .3), border=FALSE) } # end for: i

?

圖 顯示了上面估計的時變VAR參數的一部分。 頂行顯示估計點8,15和18的VAR參數的可視化。藍色實線箭頭表示正關系，紅色虛線箭頭表示負關系。 箭頭的寬度與相應參數的絕對值成比例。

如果您有任何疑問，請在下面發表評論。

總結

以上是生活随笔為你收集整理的网页中加载obj模型比较慢_R语言估计时变VAR模型时间序列的实证研究分析案例...的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。