Code - Note 6

Estimation based on large average for models on trees

Stefka Asenova

2023-02-11

Load the data.

#load("~/gremes/data/SeineData.RData")
data("SeineData", package = "gremes")

head(Seine)
#>       Paris      Meaux      Melun   Nemours     Sens
#> 1  2.124338  -1.736640  9.0907845 0.5187825 14.26062
#> 2  9.646991   2.632994  8.7440743 0.6394548 15.47040
#> 3 19.172176  21.584171  9.1692941 2.1770178 15.26537
#> 4  5.357301 -14.759895 -0.4941262 1.5233623 13.47500
#> 5  6.044298 -20.772417  1.0182723 2.1690976 13.39804
#> 6  3.805912 -40.779572 -7.3982814 1.2666607 14.60144

Generate the graph and name the nodes. Note that assigning names to nodes is crucial. Also the names of the nodes should correspond to the names of the columns in the dataset.

seg<- graph(c(1,2,
              2,3,
              2,4,
              4,5,
              5,6,
              5,7), directed = FALSE)
name_stat<- c("Paris", "2", "Meaux", "Melun", "5", "Nemours", "Sens")
seg<- set.vertex.attribute(seg, "name", V(seg), name_stat) # 
plot(seg)

Create an object combining the graph and the data. We can extract from it the nodes with latent variables.

tobj<- Tree(seg, Seine) 
#> From validate.Network: Edges have been assigned names
#> From validate.Network: There are nodes with latent variables
#> From validate.Network: Edges have been assigned names
#> From validate.Network: There are nodes with latent variables
Uc<- getNoDataNodes(tobj)

The messages say that edges are also named, and that there are latent variables in the model.

Estimate according to the Method of Moments: first create an object of the appropriate class and then call the method estimate.

mme_ave<- MMEave(seg)
#> From HRMnetwork: Edges have been assigned names
estimate(mme_ave, Seine, k_ratio=0.2)$depParams
#> From validate.Network: There are nodes with latent variables
#> From validate.Network: There are nodes with latent variables
#> From setRootDepSet.RootDepSet: The order of the subset must correspond to the
#>             order of its corresponding root
#> 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
#>        e1        e2        e3        e4        e5        e6 
#> 0.3991855 1.7674395 0.7114457 0.8146106 1.7026741 0.9087931

Estimate according to the method of maximum likelihood: first create an object of the appropriate class and then call the method estimate.

mle_ave<- MLEave(seg)
#> From HRMnetwork: Edges have been assigned names
estimate(mle_ave, Seine, k_ratio=0.2)$depParams
#> From validate.Network: There are nodes with latent variables
#> From validate.Network: There are nodes with latent variables
#> From setRootDepSet.RootDepSet: The order of the subset must correspond to the
#>             order of its corresponding root
#>        e1        e2        e3        e4        e5        e6 
#> 0.3416488 1.6908232 0.9294585 1.0357446 1.6046801 0.9647892

The estimates are slightly different from each other. Note also that the estimates are a bit higher compared to MME and MLE in Estimation - Notes 1,2,3, but they all display similar pattern: edge weights two and five are above one and \(e_3\), \(e_4\) and \(e_6\) have similar values, around 0.6-0.9.

Although local estimation is possible using a collection of subsets around every node with observable variable, we do not pursue this method here.