davidbolin / MetricGraph

https://davidbolin.github.io/MetricGraph/
11 stars 0 forks source link

Inconsistent parameter estimates with unordered PtE with graph_lme(...) #1

Closed fpavone closed 9 months ago

fpavone commented 9 months ago

The function graph_lme(...) seems to have a different behavior depending whether the observation locations (PtE) are ordered (according to the edge indexing) or shuffled. In particular, when data are shuffled the model seems to be not correctly fitted, resulting in inconsistent estimates of parameters tau and kappa, for both the isoCov and the WM1 models.

I think this might be an unexpected behavior since I haven't found any explicit recommendation to order the data locations in the documentation.

An example:

library(MetricGraph)

set.seed(35242)

## Complete Graph (10 total nodes)
tot_nodes <- 10
V_comp <- matrix(0, tot_nodes, 2)
for(i in 1:tot_nodes){
  V_comp[i,] <- c(cos(2*pi/tot_nodes*(i-1)), sin(2*pi/tot_nodes*(i-1)))
}
E_comp <- t(combn(tot_nodes,2))
comp_graph <- metric_graph$new(V = V_comp, E = E_comp)

## Simulate random locations
nobs_edge <- 10

PtE <- NULL
for (i in 1:nrow(comp_graph$E)) {
  PtE <- rbind(PtE, cbind(rep(i, nobs_edge), runif(nobs_edge)))
}

PtE_shuffle <- PtE[sample(nrow(PtE), size = nrow(PtE), replace = FALSE),]

tau <- 1
kappa <- 1
alpha <- 1

## Simulate data
u <- sample_spde(tau = tau,
                 kappa = kappa,
                 graph = comp_graph,
                 alpha = alpha,
                 sigma_e = 0,
                 PtE = PtE)
data <- data.frame(edge_number = PtE[,1],
                   distance_on_edge =  PtE[,2],
                   u = u)

u_shuffle <- sample_spde(tau = tau,
                         kappa = kappa,
                         graph = comp_graph,
                         alpha = alpha,
                         sigma_e = 0,
                         PtE = PtE_shuffle)
data_shuffle <- data.frame(edge_number = PtE_shuffle[,1],
                           distance_on_edge =  PtE_shuffle[,2],
                           u = u_shuffle)

## Fit with ordered PtE
comp_graph$clear_observations()
comp_graph$add_observations(data = data, 
                            normalized = TRUE)
res_exp <- graph_lme(u ~ -1, 
                     graph = comp_graph, 
                     model = list(type = "isoCov"))
res_WM <- graph_lme(u ~ -1, 
                    graph = comp_graph, 
                    model = 'WM1')

## Fit with shuffled PtE
comp_graph$clear_observations()
comp_graph$add_observations(data = data_shuffle, 
                            normalized = TRUE)
res_exp_shuffle <- graph_lme(u ~ -1, 
                             graph = comp_graph, 
                             model = list(type = "isoCov"))
res_WM_shuffle <- graph_lme(u ~ -1, 
                            graph = comp_graph, 
                            model = 'WM1')

## Comparison
print(res_exp$coeff$random_effects)
#      tau     kappa 
# 1.0140101 0.5322749 
print(res_exp_shuffle$coeff$random_effects)
#         tau       kappa 
# 0.04955096 14.61043295 

print(res_WM$coeff$random_effects)
#     tau    kappa 
# 0.963405 0.730961 
print(res_WM_shuffle$coeff$random_effects)
 #       tau       kappa 
 # 10.07456 11910.14242 
vpnsctl commented 9 months ago

Hello Federico,

Thanks for your message, you really found one bug. However, the bug was not in graph_lme(). It was in sample_spde(). I have now fixed it. Can you confirm that it is working now? Please, install the most recent version of the package with:

remotes::install_github("davidbolin/metricgraph", ref = "devel")

Also, please note that the parameters you estimate when using "isoCov" do not need to be "close" to the ones you chose, since you sampled them from the Whittle-Matern model. This should be true, however, when you estimate using the "WM1" model.

I will close the issue. If you still have problem, please send another message and we can reopen the issue.

Best, Alexandre