Estimando modelos autorregressivos de vetores variáveis ​​(VAR) no tempo

cupom com desconto - o melhor site de cupom de desconto cupomcomdesconto.com.br


[Esteartigofoipublicadopelaprimeiravezem[Thisarticlewasfirstpublishedon Jonas Haslbeck – r, e gentilmente contribuiu para os R-blogueiros]. (Você pode relatar um problema sobre o conteúdo desta página aqui)


Deseja compartilhar seu conteúdo com R-blogueiros? clique aqui se você tiver um blog ou aqui se não tiver.

Modelos para assuntos individuais estão se tornando cada vez mais populares na pesquisa psicológica. Um dos motivos é que é difícil fazer inferências a partir de dados entre pessoas para processos internos. Outra é que as séries temporais obtidas de indivíduos estão se tornando cada vez mais disponíveis devido à onipresença dos dispositivos móveis. O objetivo central da chamada modelagem idiográfica é explorar a dinâmica intrapessoal subjacente aos fenômenos psicológicos. Com esse objetivo em mente, muitos pesquisadores se propuseram a analisar as dependências multivariadas em séries temporais dentro da pessoa. O modelo mais simples e mais popular para essas dependências é o modelo VAR (Regressão Automática de Vetor) de primeira ordem, no qual cada variável no momento atual é prevista por (uma função linear de) todas as variáveis ​​(incluindo ela mesma) no tempo anterior. ponto.

Uma premissa chave do modelo VAR padrão é que seu parâmetro não muda com o tempo. No entanto, muitas vezes se interessa exatamente por essas mudanças ao longo do tempo. Por exemplo, alguém pode estar interessado em relacionar alterações nos parâmetros com outras variáveis, como alterações no ambiente de uma pessoa. Este poderia ser um novo emprego, as estações do ano ou o impacto de uma pandemia global. Em projetos menos exploratórios, pode-se examinar o impacto que determinadas intervenções (por exemplo, medicação ou terapia) têm nas interações entre os sintomas.

Nesta publicação, dou uma breve visão geral de como estimar um modelo VAR com variação de tempo com a abordagem de suavização de kernel, que discutimos neste recente tutorial. Este método baseia-se no pressuposto de que os parâmetros podem mudar suavemente ao longo do tempo, o que significa que os parâmetros não podem “pular” de um valor para outro. Depois, concentro-me em como estimar e analisar esse tipo de modelo VAR de variação temporal com o pacote R mgm.

Estimando modelos com variação de tempo via suavização de kernel

A idéia central da abordagem de suavização do kernel é a seguinte: Escolhemos pontos de tempo igualmente espaçados ao longo de toda a série temporal e, em seguida, estimamos modelos “locais” em cada um desses pontos de tempo. Todos os modelos locais, juntos, constituem o modelo de variação temporal. Com modelos “locais”, queremos dizer que esses modelos são amplamente baseados em pontos de tempo próximos ao ponto em questão. Isso é obtido ponderando as observações adequadamente durante a estimativa dos parâmetros. Esta ideia é ilustrada para um conjunto de dados de brinquedos na figura a seguir:

Aqui ilustramos apenas a estimativa do modelo local em $ t = 3 $. Vemos os 10 pontos temporais desta série cronológica no painel esquerdo. A coluna $ w_ {t_e = 3} $ em vermelho indica um possível conjunto de pesos que poderíamos usar para estimar o modelo local em $ t = 3 $: os dados em momentos próximos a $ t = 3 $ obtêm o peso mais alto, e pontos de tempo mais distantes ganham um peso cada vez menor. A função que define esses pesos é mostrada no painel direito. A coluna azul no painel esquerdo e a função azul correspondente à direita indicam outra possível ponderação. Usando essa ponderação, combinamos menos observações com o tempo. Isso nos permite detectar mais “variação no tempo” nos parâmetros, porque suavizamos ao longo de menos pontos no tempo. Por outro lado, no entanto, usamos menos dados, o que torna nossas estimativas menos confiáveis. Portanto, é importante escolher uma função de ponderação que atinja um bom equilíbrio entre a sensibilidade à “variação no tempo” e estimativas estáveis. No método apresentado aqui, usamos uma função de ponderação gaussiana (também chamada de núcleo), definido pelo seu desvio padrão (ou largura de banda). Voltaremos a como selecionar um bom parâmetro de largura de banda abaixo.

Leia Também  Animações na época do coronavírus

Nesta postagem do blog, eu me concentro em como estimar modelos que variam no tempo com o pacote R mgm. Para uma explicação mais detalhada do método, consulte nosso recente tutorial.

Carregando e inspecionando os dados

Para ilustrar a estimativa de modelos VAR variáveis ​​no tempo, eu uso uma série temporal do ESM de 12 variáveis ​​relacionadas ao humor que são medidas até 10 vezes por dia por 238 dias consecutivos (para detalhes sobre esse conjunto de dados, consulte Kossakowski et al. (2017)). As perguntas são “Sinto-me relaxado”, “Sinto-me abatido”, “Sinto-me irritado”, “Sinto-me satisfeito”, “Sinto-me sozinho”, “Sinto-me ansioso”, “Sinto-me entusiasmado”, “Sinto-me suspeito” , “Sinto-me alegre”, “sinto-me culpado”, “sinto-me indeciso” e “sinto-me forte”. Cada pergunta é respondida em uma escala Likert de 7 pontos, variando de “não” a “muito”.

O conjunto de dados é carregado com o mgm-pacote. Primeiro, subconjuntamos as 12 variáveis ​​de humor:

library(mgm)
mood_data <- as.matrix(symptom_data$data[, 1:12])
mood_labels <- symptom_data$colnames[1:12]
colnames(mood_data) <- mood_labels
time_data <- symptom_data$data_time

Vemos que o conjunto de dados possui 1476 observações:

dim(mood_data)
## [1] 1476   12
head(mood_data)
##      Relaxed Down Irritated Satisfied Lonely Anxious Enthusiastic Suspicious Cheerful
## [1,]       5   -1         1         5     -1      -1            4          1        5
## [2,]       4    0         3         3      0       0            3          1        4
## [3,]       4    0         2         3      0       0            4          1        4
## [4,]       4    0         1         4      0       0            4          1        4
## [5,]       4    0         2         4      0       0            4          1        4
## [6,]       5    0         1         4      0       0            3          1        3
##      Guilty Doubt Strong
## [1,]     -1     1      5
## [2,]      0     1      4
## [3,]      0     2      4
## [4,]      1     1      4
## [5,]      1     2      3
## [6,]      1     2      3

time_data contém informações temporais sobre cada medição. Usaremos o dia em que a medição ocorreu (dayno), o prompt de medição (beepno) e o horário geral (time_norm)

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

Selecionando a largura de banda ideal

Uma maneira de selecionar um bom parâmetro de largura de banda é ajustar modelos que variam no tempo com diferentes parâmetros candidatos de largura de banda em um conjunto de dados de treinamento e avaliar seu erro de previsão em um conjunto de dados de teste. A função bwSelect() implementa esse esquema de seleção de largura de banda. Aqui não mostramos a especificação dessa função, pois ela possui os mesmos argumentos de entrada que a função tvmvar() que descrevemos em um momento abaixo, além de uma sequência candidata de valores de largura de banda e algumas especificações de como dividir os dados em dados de treinamento e teste. Além disso, a seleção da largura de banda orientada a dados pode levar uma quantidade considerável de tempo para ser executada, o que não permitiria a execução do código durante a leitura da postagem do blog. Portanto, neste tutorial, apenas corrigimos a largura de banda para o valor retornado por bwSelect(). No entanto, você pode encontrar o código para executar a seleção de largura de banda com bwSelect() no conjunto de dados atual aqui).

Leia Também  Instale o R sem suporte para longas dobras (noLD) no Ubuntu
bandwidth <- .34

Estimando modelos VAR com variação temporal

Agora podemos especificar a estimativa do modelo VAR variável no tempo. Fornecemos os dados como entrada e especificamos o tipo de variáveis ​​e quantas categorias eles têm com o type e level argumentos. No nosso exemplo de conjuntos de dados, todas as variáveis ​​são contínuas e, portanto, definimos type = rep("g", 12) para Gaussiano contínuo e defina o número de categorias como 1 por convenção. Optamos por selecionar os parâmetros de regularização com validação cruzada com lambdaSel = "CV"e especificamos que o modelo VAR deve incluir um único atraso com lags = 1. Os argumentos beepvar e dayvar forneça o dia e o número de notificações em um determinado dia para cada medição, o que é necessário para especificar a matriz de design do VAR. Além disso, fornecemos os carimbos de hora de todas as medições com timepoints = time_data$time_norm para dar conta de medições ausentes. Observe, no entanto, que ainda assumimos um tamanho de atraso constante igual a 1. Os carimbos de hora são usados ​​apenas para garantir que a ponderação realmente dê a esses pontos de tempo o peso mais alto mais próximo do ponto de estimativa atual (para obter detalhes, consulte a Seção 2.5 neste documento). ) Até agora, a especificação é a mesma que para o mvar() função que se encaixa em modelos VAR mistos estacionários.

Para o modelo de variação temporal, precisamos especificar dois argumentos adicionais. Primeiro, com estpoints = seq(0, 1, length = 20) especificamos que gostaríamos de estimar 20 modelos locais ao longo de toda a série temporal (que é normalizada para [0,1]) O número de pontos de estimativa pode ser escolhido arbitrariamente grande, mas em algum momento a adição de mais pontos de estimativa não compensa os custos computacionais adicionais, porque os modelos locais subsequentes são essencialmente idênticos. Por fim, especificamos a largura de banda com o bandwidth argumento.

# Estimate Model on Full Dataset
set.seed(1)
tvvar_obj <- tvmvar(data = mood_data,
                    type = rep("g", 12),
                    level = rep(1, 12), 
                    lambdaSel = "CV",
                    lags = 1,
                    beepvar = time_data$beepno,
                    dayvar = time_data$dayno,
                    timepoints = time_data$time_norm, 
                    estpoints = seq(, 1, length = 20), 
                    bandwidth = bandwidth,
                    pbar = FALSE)
## Note that the sign of parameter estimates is stored separately; see ?tvmvar

Podemos colar o objeto de saída no console

# Check on how much data was used
tvvar_obj
## mgm fit-object 
## 
## Model class:  Time-varying mixed Vector Autoregressive (tv-mVAR) model 
## Lags:  1 
## Rows included in VAR design matrix:  876 / 1476 ( 59.35 %) 
## Nodes:  12 
## Estimation points:  20

que fornece um resumo do modelo e também mostra quantas linhas estavam na matriz de design do VAR (876) em comparação com o número de pontos no tempo no conjunto de dados (1476). O número anterior é menor, porque um modelo VAR (1) só pode ser estimado se, por um determinado ponto no tempo, também estiver disponível o atraso no ponto 1 no tempo anterior. Esse não é o caso da primeira medição em um determinado dia ou se houver respostas ausentes durante o dia.

Computando erros de previsão variáveis ​​no tempo

Da mesma forma que os modelos VAR estacionários, podemos calcular erros de previsão. Isso pode ser feito com o predict() , que pega o objeto de modelo, os dados e duas variáveis ​​indicando o número do dia e o número da notificação. O fornecimento de dados e variáveis ​​de notificação independentemente do objeto de modelo permite calcular erros de previsão para novas amostras.

Leia Também  Plataformas de mineração Zigmabit oferecem o ROI mais rápido do mercado

O argumento errorCon = c("R2", "RMSE") especifica que a proporção da variação explicada ($ R ^ 2 $) e o erro médio quadrático da raiz (RMSE) devem ser retornados como erros de previsão. O argumento final tvMethod especifica como os erros de previsão com variação temporal devem ser calculados. A opção tvMethod = "closestModel" faz previsões para um ponto no tempo usando o modelo local mais próximo a ele. A opção escolhida aqui, tvMethod = "weighted", fornece uma média ponderada das previsões de todos os modelos locais, ponderada usando a função de ponderação centralizada no local do ponto de tempo em questão. Normalmente, ambos os métodos dão resultados muito semelhantes.

cupom com desconto - o melhor site de cupom de desconto cupomcomdesconto.com.br
pred_obj <- predict(object = tvvar_obj, 
                    data = mood_data, 
                    beepvar = time_data$beepno,
                    dayvar = time_data$dayno,
                    errorCon = c("R2", "RMSE"),
                    tvMethod = "weighted")

A saída principal são os dois objetos a seguir:
pred_obj$tverrors é uma lista que inclui os erros de estimativa para cada ponto de estimativa / modelo local; pred_obj$errors contém o erro médio entre os pontos de estimativa.

Visualizando partes do modelo

O modelo VAR (1) de variação temporal consiste em parâmetros $ (p + p ^ 2) times E $, em que $ p $ é o número de variáveis ​​e $ E $ é o número de pontos de estimativa. Visualizar todos os parâmetros de uma vez é, portanto, desafiador. Em vez disso, pode-se escolher os parâmetros que são de maior interesse para a questão de pesquisa em questão. Aqui, escolhemos duas visualizações diferentes. Primeiro, usamos o qgraph() função do pacote R qgraph para inspecionar os parâmetros de interação VAR nos pontos de estimativa 1, 10 e 20:

library(qgraph)

par(mfrow=c(1,3))
for(tp in c(1,10,20)) qgraph(t(tvvar_obj$wadj[, , 1, tp]), 
                             layout = "circle",
                             edge.color = t(tvvar_obj$edgecolor[, , 1, tp]), 
                             labels = mood_labels, 
                             mar = rep(5, 4), 
                             vsize=14, esize=15, asize=13,
                             maximum = .5, 
                             pie = pred_obj$tverrors[[tp]][, 3],
                             title = paste0("Estimation point = ", tp))

parcela do pedaço sem nome-pedaço-8

Vemos que alguns parâmetros nos modelos VAR variam consideravelmente ao longo do tempo. Por exemplo, o efeito de autocorrelação do Relaxed parece estar diminuindo ao longo do tempo, o efeito positivo de Strong on Satisfied aparece apenas no ponto de estimativa 20, e também o efeito negativo de Satisfied on Guilty só aparece no ponto de estimativa 20.

Podemos ampliar esses parâmetros individuais, plotando-os em função do tempo:

# Obtain parameter estimates with sign
par_ests <- tvvar_obj$wadj[, , 1, ]
par_ests[tvvar_obj$edgecolor[, , 1, ]=="red"] <- par_ests[tvvar_obj$edgecolor[, , 1, ]=="red"] * -1

# Select three parameters to plot
m_par_display <- matrix(c(1, 1, 
                          4, 12, 
                          10, 4), ncol = 2, byrow = T)
# Plotting
plot.new()
par(mar = c(4, 4, , 1))
plot.window(xlim=c(1, 20), ylim=c(-.25, .55))
axis(1, c(1, 5, 10, 15, 20), labels = T)
axis(2, c(-.25, , .25, .5), las = 2)
abline(h = , col = "grey", lty = 2)
title(xlab = "Estimation points", cex.lab = 1.2)
title(ylab = "Parameter estimate", cex.lab = 1.2)

for(i in 1:nrow(m_par_display)) {
  par_row <- m_par_display[i, ]
  P1_pointest <- par_ests[par_row[1], par_row[2], ]
  lines(1:20, P1_pointest, lwd = 2, lty = i) 
}

legend_labels <- c(expression("Relaxed"["t-1"]  %->%  "Relaxed"["t"]),
                   expression("Strong"["t-1"]  %->%  "Satisfied"["t"]),
                   expression("Satisfied"["t-1"]  %->%  "Guilty"["t"]))
legend(1, .49, 
       legend_labels,
       lwd = 2, bty = "n", cex = 1, horiz = T, lty = 1:3)

parcela do pedaço sem nome-pedaço-9

Vemos que o efeito de Relaxed sobre si mesmo no próximo ponto no tempo é relativamente forte no início da série temporal, mas depois diminui para zero e permanece zero em torno do ponto de estimativa 13. O efeito cruzado de Strong on Satisfied no o próximo ponto no tempo é igual a zero até o ponto de estimativa 9, mas parece aumentar monotonicamente. Finalmente, o efeito cruzado de Satisfied on Guilty também é igual a zero até o ponto de estimativa 13 e depois diminui monotonicamente.

Estabilidade das estimativas

Semelhante aos modelos estacionários, pode-se avaliar a estabilidade de parâmetros variáveis ​​no tempo usando distribuições de amostragem com rapidez. o mgm pacote permite fazer isso com o resample() função. Aqui você pode encontrar código sobre como fazer isso, por exemplo, neste tutorial.

Variando no tempo ou não?

Claramente, “variação no tempo” é um continuum que vai de estacionário a extremamente variável no tempo. No entanto, em alguns casos, pode ser necessário decidir se os parâmetros de um modelo VAR variam de maneira confiável no tempo. Para chegar a essa decisão, pode-se usar um teste de hipótese com a hipótese nula de que o modelo não varia no tempo. Aqui está uma maneira de realizar esse teste de hipótese: Começa-se ajustando um modelo estacionário de VAR aos dados e depois simula repetidamente os dados desse modelo estimado. Para cada um desses conjuntos de dados de séries temporais simulados, calcula-se um erro de previsão em pool para o modelo de variação temporal. A distribuição desses erros de previsão serve como distribuição de amostragem de erros de previsão sob a hipótese nula. Agora, é possível calcular o erro de estimativa combinada do modelo VAR de variação temporal nos dados empíricos e usá-lo como uma estatística de teste. Este teste é explicado em mais detalhes aqui, e o código para implementar esse teste para o conjunto de dados usado neste tutorial pode ser encontrado aqui.

Sumário

Nesta postagem do blog, mostrei como estimar um modelo de VAR com variação de tempo com uma abordagem de suavização de kernel, baseada na suposição de que todos os parâmetros são uma função suave do tempo. Além de estimar o modelo, discutimos a seleção de um parâmetro de largura de banda apropriado, como calcular erros de previsão (variáveis ​​no tempo) e como visualizar diferentes aspectos do modelo. Por fim, forneci indicadores para o código que mostra como avaliar a estabilidade das estimativas via bootstrapping e como executar um teste de hipótese que se pode usar para selecionar entre modelos VAR estacionários e variáveis ​​no tempo.


Gostaria de agradecer a Fabian Dablander por seus comentários sobre este post do blog.

var vglnk = {key: ‘949efb41171ac6ec1bf7f206d57e90b8’};

(função (d, t) {
var s = d.createElement
s.src = ‘//cdn.viglink.com/api/vglnk.js’;
var r = d.getElementsByTagName
} (documento, ‘script’));

Para Deixe um comentário para o autor, siga o link e comente no blog: Jonas Haslbeck – r.

R-bloggers.com oferece atualizações diárias por email sobre notícias e tutoriais do R sobre o aprendizado do R e muitos outros tópicos. Clique aqui se você deseja publicar ou encontrar um emprego em ciência da dados / R.


Deseja compartilhar seu conteúdo com R-blogueiros? clique aqui se você tiver um blog ou aqui se não tiver.



cupom com desconto - o melhor site de cupom de desconto cupomcomdesconto.com.br