We do not have a real dataset according to a given block graph which is not a tree. Therefore we first generate the observations. How is the data generated is described in Vignette “Additional functionalities”.
Create the block graph.
graph(c(1,3,1,2,2,3,
g<-3,4,4,5,5,3,
3,7,3,6,6,7), directed=FALSE)
set.vertex.attribute(g, "name", V(g), c("a", "b", "c", "d", "e", "f", "g"))
g<-plot(g)
Create the edge weights to be assigned to the edges of the graph.
# all deltas are squares already
c(0.2, 0.8, 0.6) # d_13^2, d_12^2, d_23^2
C1<- c(0.3, 0.5, 0.1) # d_34^2, d_45^2, d_35^2
C2<- c(0.4, 0.05, 0.25) # d_37^2, d_36^2, d_67^2 C3<-
Attach the edge weights to the edges.
HRMBG(g)
hrmbgobj<-#> From HRMnetwork: Edges have been assigned names
setParams(hrmbgobj, c(C1, C2, C3))
hrmbgobj<-#> From setParams.HRMtree: Names have been attributed to the vector 'value' in the order corresponding to the order of the edges: The fist element has the name of the first edge, the second element the name of the second edge, etc.
#> From setParams.HRMtree: The parameters have been attached to the edges according to their names
hrmbgobj#> $graph
#> IGRAPH 0d44d7f UN-- 7 9 --
#> + attr: name (v/c), name (e/c)
#> + edges from 0d44d7f (vertex names):
#> [1] a--c a--b b--c c--d d--e c--e c--g c--f f--g
#>
#> $depParams
#> e1 e2 e3 e4 e5 e6 e7 e8 e9
#> 0.20 0.80 0.60 0.30 0.50 0.10 0.40 0.05 0.25
#>
#> attr(,"class")
#> [1] "HRMnetwork" "HRMBG" "MME"
Create the matrix \(\Lambda\), whose entry \(\lambda_{ij}\) is the sum of the edge weights on the unique shortest path between node \(i\) and node \(j\).
HRLambda(hrmbgobj)
hrmlam<-
hrmlam#> a b c d e f g
#> a 0.00 0.80 0.20 0.50 0.30 0.25 0.60
#> b 0.80 0.00 0.60 0.90 0.70 0.65 1.00
#> c 0.20 0.60 0.00 0.30 0.10 0.05 0.40
#> d 0.50 0.90 0.30 0.00 0.50 0.35 0.70
#> e 0.30 0.70 0.10 0.50 0.00 0.15 0.50
#> f 0.25 0.65 0.05 0.35 0.15 0.00 0.25
#> g 0.60 1.00 0.40 0.70 0.50 0.25 0.00
Generate 1000 observations from Huesler-Reiss distribution with parameter matrix \(\Lambda\) and some independent random noise.
rHRM(hrmbgobj, hrmlam, 1000, noise = TRUE) X<-
Create the subsets for local estimation.
RootDepSet()
rdsobj<- setRootDepSet(rdsobj, subset = list(c("a", "b", "d", "e", "f", "g"),
rdsobj<-c("a", "b", "d", "e", "f", "g"),
c("a", "b", "d", "e", "f", "g"),
c("a", "b", "d", "e", "f", "g"),
c("a", "b", "d", "e", "f", "g"),
c("a", "b", "d", "e", "f", "g")),
c("a", "b", "d", "e", "f", "g"))
#> From setRootDepSet.RootDepSet: The order of the subset must correspond to the
#> order of its corresponding root
Estimate the model treating the third variable as latent: create first an object of class HRMBG
and then use on it the method estimate
.
HRMBG(g)
hrmbg<-#> From HRMnetwork: Edges have been assigned names
suppressMessages(estimate(hrmbg, X[,-3], rdsobj, k_ratio=0.2))
hrmbg<-$depParams
hrmbg#> e1 e2 e3 e4 e5 e6 e7
#> 0.13332429 0.38743378 0.22935665 0.19510273 0.33491670 0.08015118 0.22420946
#> e8 e9
#> 0.06361528 0.17779621
We have suppressed the messages. The messages are informative, if they don’t stop the process they are not errors.