AncReg: Ancestor Regression

View source: R/AncReg.R

AncRegR Documentation

Ancestor Regression

Description

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.

Usage

AncReg(x, degree = 0, targets = colnames(x), f = function(x) x^3)

Arguments

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.

Value

An object of class "AncReg" containing:

z.val

A numeric matrix of test statistics.

p.val

A numeric matrix of p-values.

References

\insertAllCited

See Also

summary.AncReg, instant_graph, summary_graph, instant_p.val, summary_p.val

Examples

##### 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)

AncReg documentation built on Aug. 8, 2025, 7:48 p.m.