Identifying lead/lags using multivariate regression analysis

Question

I have three time-series variables (x,y,z) measured in 3 replicates. x and z are the independent variables. y is the dependent variable. t is the time variable. All the three variables follow diel variation, they increase during the day and decrease during the night. An example with a simulated dataset is below.

   library(nlme)
    library(tidyverse)
    n <- 100 
    t <- seq(0,4*pi,,100)
    a <- 3
    b <- 2
    c.unif <- runif(n)
    amp <- 2

datalist = list()
for(i in 1:3){
y <- 3*sin(b*t)+rnorm(n)*2
x <- 2*sin(b*t+2.5)+rnorm(n)*2
z <- 4*sin(b*t-2.5)+rnorm(n)*2
data = as_tibble(cbind(y,x,z))%>%mutate(t = 1:100)%>% mutate(replicate = i)
datalist[[i]] <- data
}
df <- do.call(rbind,datalist)
ggplot(df)+
  geom_line(aes(t,x),color='red')+geom_line(aes(t,y),color='blue')+
  geom_line(aes(t,z),color = 'green')+facet_wrap(~replicate, nrow = 1)+theme_bw()

I can identify the lead/lag of y with respect to x and z individually. This can be done with ccf() function in r. For example

ccf(x,y)
ccf(z,y)

But I would like to do it in a multivariate regression approach. For example, nlme package and lme function indicates y and z are negatively affecting x

lme = lme(data = df, y~ x+ z , random=~1|replicate, correlation = corCAR1( form = ~ t| replicate))

It is impossible (in actual data) that x and z can negatively affect y.

I need the time-lead/lag and also I would like to get the standardized coefficient (t-value to compare the effect size), both from the same model.

Is there any multivariate model available that can give me the lead/lag and also give me regression coefficient?

Answer 1

We might be considering the " statistical significance of Cramer Rao estimation of a lower bound". In order to find Xbeta-Xinfinity, taking the expectation of Xbeta and an assumed mean neu; will yield a variable, neu^squared which can replace Xinfinity. Using the F test-likelihood ratio, the degrees of freedom is p2-p1 = n-p2. Put it this way, the estimates are n=(-2neu^squared/neu^squared+n), phi t = y/Xbeta and Xbeta= (y-betazero)/a. The point estimate is derived from y=aXbeta + b: , Xbeta. The time lead lag is phi t and the standardized coefficient is n. The regression generates the lower bound Xbeta, where t=beta. Spectral analysis of the linear distribution indicates a point estimate beta zero = 0.27 which is a significant peak of variability. Scaling Xbeta by Betazero would be an appropriate idea.