I would suggest you to model your problem as a graph. You could convert your nxm matrix (n actions over m time points) into a nxn similarity matrix for every pairwise linear correlation between the vector time-series representing your actions. When the time-series of two actions are highly correlated it means that these actions are highly synchronised, which I think is what you want.
Such similarity/correlation matrix can be easily constructed in R using:
Y = cor(X, method="pearson")
where X is your original matrix. At this point you have a similarity/correlation matrix which you can visualise as a heatmap in order to see how your actions correlate. You can further convert this correlation matrix Y into an adjacency/sparse matrix by setting a threshold value below which all slots in the matrix will take zero values:
Z = Y[Y>=0.8]
I used a random threshold value of 0.8 but you should set that value empirically by "spying" your matrix at different thresholds (see here how to do that in R). Z is now a binary adjacency/sparse matrix representing an unweighted undirected graph. A connection between two nodes in that graph implies the presence of a high co-expression between two actions. You could alternatively create a weighted graph by explicitly setting all values below your threshold to zero while keeping the values above your threshold to their original correlation values instead of converting to 1:
Z = Y
Z[Z<0.8] = 0
By now you have a graph where:
- Every node represents an action
- Every undirected edge (connection) between nodes i and j indicates that action i is synchronous/correlated/co-expressed with action j. This connection may or may not have a weight (that's up to you).
Once you have that setup you can basically perform any graph-related algorithm to gain insights on your data, e.g.:
manually isolate every node along with all its parent nodes (adjacent nodes pointing to that node). That would allow you to make up a set of rules based a node's number of parents and their corresponding weights.
Find all possible paths between every pair of nodes {i,j} and evaluate their lengths by adding up all weights composing them. That would give you a metric of strength of transition (or causality) between all actions (nodes) in your graph.
model your graph as a Bayesian network. That would require to normalise and convert your weights into transition probabilities between states (actions) and would offer you a probabilistic inferential graphical model of your data.
Use graph metrics (degree centrality, betweeness, etc.) along with graph clustering algorithms (e.g. community structure partition) in order to detect important actions and groups of actions.
For your generic network-graph theoretical needs as well as solution 1 and/or solution 2, I would recommend packages igraph, rgraphviz, or networkd3. For Bayesian networks I would suggest gRain.
I hope that helps.