Thanks so much! It is so nice that you have answered.
Yes, I did two interesting experiments. The first one is that I used an arbitrarily chosen simple normal nugget effect (N(0, sigma)) to generate one single 100*100 grid realization, without intercept and any covariates, just the nugget effect (N(0, sigma)). Then, I averaged those 100*100 grids to 10*10 coarser averaged areal value, use INLA() to function to fit those coarser 10*10 averaged areal values. Even for this one single realization, the precision for the gaussian observation is 1/100 times the sigma of the finer resolution (100*100) nugget variance that I arbitrarily chose in the first place. This is aligned with a perfect world without any spatial random fields or SPDE model or numerical representation involved.
However, if I added a Matérn covariance spatial random field in the simulation model, plus a simple intercept and a simple normal nugget effect like above, and fit the model on the composite averaged observations. The INLA() function cannot give accurate precision for the observation value for the averaged coarser areal observations (very different from 1/100 times the fine resolution nugget variance, which is expected). Thus, I am wondering if it is caused by the interaction with the numerical representation of the spatial random field itself, since the INLA() function can derive a good estimation on the pure composite nugget effects without any spatial random fields involved. Looking forward to your reply!
Thanks,
Bowen