1
$\begingroup$

I would like to test differences in mean values of a Climate_variable between 4 different climatic Scenarios that are projected for many (3000) surveyed Sites for which I know the coordinates (Lon and Lat). My database is structured as follows:

Climate_variable Scenario Site Lon Lat
0.85 A S1 19.37 41.85
1.13 B S1 19.37 41.85
1.55 C S1 19.37 41.85
2.12 D S1 19.37 41.85
0.74 A S2 21.42 40.36
1.08 B S2 21.42 40.36
1.44 C S2 21.42 40.36
2.01 D S2 21.42 40.36

Because I have 4 observations for each Site but I am not interested in this effect, I wanted to go for a Linear Mixed Model with Site as random effect. However, climatic variables are often highly spatially autocorrelated so I also wanted to add a spatial autocorrelation structure using the coordinates of the sites.

This is the code I run using the nlme package on R:

mod <- lme(Climate_variable ~ Scenario,
           correlation = corExp(form = ~Lon+Lat),
           random = list(~1|Site))

However, there is an error message saying cannot have zero distances in "corSpatial", I guess because I have 4 observation with the same coordinates for each site.

So I tried this instead:

mod <- lme(Climate_variable ~ Scenario,
           correlation = corExp(form = ~Lon+Lat | Scenario),
           random = list(~1|Site))

But I get another error message: incompatible formulas for groups in 'random' and 'correlation'

How could I account for both the spatial autocorrelation and the non-independance of the observations for different Scenarios on the same Site?

$\endgroup$
2
  • 1
    $\begingroup$ It looks like Site is completely determined by Lon and Lat. If so, then it doesn't make sense to include both Lat/Lon as well as Site in the model - this is including redundant information. I would remove Site and just retain Lon and Lat. $\endgroup$
    – mkt
    Jul 11, 2022 at 9:17
  • 1
    $\begingroup$ @mkt Indeed, I finaly opted for another model specification without redundancy! $\endgroup$ Feb 27, 2023 at 15:23

2 Answers 2

1
$\begingroup$

I have to tell you that i do it different. First i make the null model indicating fixed and random. If you put site as random, you are telling the program to control for the different between sites, that i think is not what you want. right? You want to control for the difference between the scenarios. Then scenario is random. example from the package.

null.model <- lme(fixed = thick ~ 1, data = spdata, random = ~ 1 | dummy)

exp.sp <- update(null.model, correlation = corExp(1, form = ~ east + north), method = "ML")

Your error."There is an error message saying cannot have zero distances in "corSpatial", I guess because I have 4 observation with the same coordinates for each site." Exactly is because you have the same coordinates and when the program try to calculate the geographic distance to check for autocorrelation give 0, that for the program sound absurd. I think this method is not the correct if you only have one site. This is for checking spatial autocorrelation. It is better you only check for correlation in your scenarios in a simple way.

$\endgroup$
1
  • $\begingroup$ No I am actually interrested in the Scenario effect and not in the Site effect, so I want to control for differences between sites. $\endgroup$ Feb 27, 2023 at 15:26
0
$\begingroup$

Even though this analysis might not be relevant (see comment from @mkt), this code works:

mod <- lme(Climate_variable ~ Scenario, correlation = corExp(form = ~(Lon+Lat)/Scenario), random = list(~1|Site))

$\endgroup$
2
  • $\begingroup$ I have the same question as yours. Your code could be a solution. But do you know what the slash before Scenario means in the wording of this formula, as opposed to the vertical bar in the former script? $\endgroup$ May 6, 2023 at 18:07
  • $\begingroup$ The main question remains: is this model really inadequate as @mkt suggests? One can have a random effect of the Site, which one wants to be able to estimate, and at the same time, in the same model, a correlation structure of the residuals (which will necessarily be different for each Scenario), more precise than a simple random effect. But I am not sure. If someone can confirm or deny, thank you in advance. $\endgroup$ May 6, 2023 at 18:09

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.