原文链接： http://tecdat.cn/?p=3364

原文出处：拓端数据部落公众号

加载R包和数据集

 加载包后，咱们将此数据集中蕴含的12个情绪变量进行子集化：

mood\_data <- as.matrix(symptom\_data$data\[, 1:12\]) # 子集变量mood\_labels <- symptom\_data$colnames\[1:12\] # 子集变量标签colnames(mood\_data) <- mood\_labelstime\_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_norm1 13/08/12 226 1 08:58 08:58:56 09:00:15 0.0000000002 14/08/12 227 5 14:32 14:32:09 14:33:25 0.0051648743 14/08/12 227 6 16:17 16:17:13 16:23:16 0.0054705744 14/08/12 227 8 18:04 18:04:10 18:06:29 0.0057820975 14/08/12 227 9 20:57 20:58:23 21:00:18 0.0062857746 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正则化”时，倡议这样做，因为否则参数惩办的强度取决于预测变量的方差。 因为穿插验证计划应用随机抽取来定义，因而咱们设置种子以确保重现性。

在查看后果之前，咱们查看了1476个工夫点中有多少用于估算，这在调用输入对象的摘要中显示

> tvvar_objmgm fit-objectModel class: Time-varying mixed Vector Autoregressive (tv-mVAR) modelLags: 1Rows included in VAR design matrix: 876 / 1475 ( 59.39 %)Nodes: 12Estimation 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$errorsVariable Error.RMSE Error.R21 Relaxed 0.939 0.1552 Down 0.825 0.2973 Irritated 0.942 0.1194 Satisfied 0.879 0.2015 Lonely 0.921 0.1826 Anxious 0.950 0.0867 Enthusiastic 0.922 0.1698 Suspicious 0.818 0.2479 Cheerful 0.889 0.20010 Guilty 0.928 0.17511 Doubt 0.871 0.26812 Strong 0.896 0.195

可视化时变VAR模型

可视化下面预计的一部分随工夫变动的VAR参数：

# 两个网络图# 获取均值图的布局Q <- qgraph(t(mean_wadj), DoNotPlot=TRUE)saveRDS(Q$layout, "Tutorials/files/layout_mgm.RDS")#在选定的固定工夫点绘制图形tpSelect <- c(2, 10, 18)# 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"\] <- 2for(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))  } )    # 绘图暗影  polygon(x = c(1:20, 20:1), y = c(CIs\[1,\], rev(CIs\[2,\])), col=alpha(colour = cols\[i\], alpha = .3), border=FALSE)    } #

图  显示了下面预计的时变VAR参数的一部分。 顶行显示预计点8,15和18的VAR参数的可视化。蓝色实线箭头示意正关系，红色虚线箭头示意负关系。 箭头的宽度与相应参数的绝对值成比例。

 如果您有任何疑难，请在上面发表评论。 

最受欢迎的见解

1.R语言多元Logistic逻辑回归 利用案例

2.面板平滑转移回归(PSTR)剖析案例实现剖析案例实现")

3.matlab中的偏最小二乘回归（PLSR）和主成分回归（PCR）

4.R语言泊松Poisson回归模型剖析案例

5.R语言回归中的Hosmer-Lemeshow拟合优度测验

6.r语言中对LASSO回归，Ridge岭回归和Elastic Net模型实现

7.在R语言中实现Logistic逻辑回归

8.python用线性回归预测股票价格

9.R语言如何在生存剖析与Cox回归中计算IDI，NRI指标