AncReg | R Documentation |
This function performs ancestor regression for linear structural equation models \insertCiteschultheiss2024ancestorregressionstructuralvectorAncReg and vector autoregressive models \insertCiteancestorAncReg with explicit error control for false discovery, at least asymptomatically.
AncReg(x, degree = 0, targets = colnames(x), f = function(x) x^3)
x |
A named numeric matrix containing the observational data. |
degree |
An integer specifying the order of the SVAR process to be considered. Default is 0 for no time series. |
targets |
A character vector specifying the variables whose ancestors should be estimated. Default is all variables. |
f |
A function specifying the non-linearity used in the ancestor regression. Default is a cubic function. |
An object of class "AncReg" containing:
z.val |
A numeric matrix of test statistics. |
p.val |
A numeric matrix of p-values. |
summary.AncReg
, instant_graph
, summary_graph
,
instant_p.val
, summary_p.val
##### simulated example for inear structural equation models
# random DAGS for simulation
set.seed(1234)
p <- 5 #number of nodes
DAG <- pcalg::randomDAG(p, prob = 0.5)
B <- matrix(0, p, p) # represent DAG as matrix
for (i in 2:p){
for(j in 1:(i-1)){
# store edge weights
B[i,j] <- max(0, DAG@edgeData@data[[paste(j,"|",i, sep="")]]$weight)
}
}
colnames(B) <- rownames(B) <- LETTERS[1:p]
# solution in terms of noise
Bprime <- MASS::ginv(diag(p) - B)
n <- 5000
N <- matrix(rexp(n * p), ncol = p)
X <- t(Bprime %*% t(N))
colnames(X) <- LETTERS[1:p]
# fit ancestor regression
fit <- AncReg(X)
# collect ancestral p-values and graph
res <- summary(fit)
res
#compare true and estimated ancestral graph
trueGraph <- igraph::graph_from_adjacency_matrix(recAncestor(B != 0))
ancGraph <- igraph::graph_from_adjacency_matrix(res$graph)
oldpar <- par(mfrow = c(1, 2))
plot(trueGraph, main = 'true ancestral graph', vertex.size = 30)
plot(ancGraph, main = 'Ancestor Regression', vertex.size = 30)
##### SVAR-example with geyser timeseries
geyser <- MASS::geyser
# shift waiting such that it is waiting after erruption
geyser2 <- data.frame(waiting = geyser$waiting[-1], duration = geyser$duration[-nrow(geyser)])
# fit ancestor regression with 6 lags considered
fit2 <- AncReg(as.matrix(geyser2), degree = 6)
res2 <- summary(fit2)
res2
# visualize instantaneous ancestry
instGraph <- igraph::graph_from_adjacency_matrix(res2$inst.graph)
plot(instGraph, edge.label = round(diag(res2$inst.p.val[1:2, 2:1]), 2),
main = 'instantaneous effects', vertex.size = 90)
# visualize summary of lagged ancestry
sumGraph <- igraph::graph_from_adjacency_matrix(res2$sum.graph)
plot(sumGraph, edge.label = round(diag(res2$sum.p.val[1:2, 2:1]), 2),
main = 'summary graph', vertex.size = 90)
par(oldpar)
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.