knitr::opts_chunk$set(collapse=T,cache=F,eval=T) # #knitr::opts_chunk$set(eval = FALSE) # options(width = 120, max.print = 100)
The package siland aims to estimate the effects of landscape and local variables on a measurement of interest. The method enable to estimate effect intensities but also the scale of the effect for landscape variables, without any a priori. The method is based on estimation of Spatial Influence Function (SIF).The SIF function describes how the influence of a pixel/cell of a landscape variable is spatially distributed. We assume that the influence is maximal at the pixel location and decreases with the distance. The decrease is parametrised by one parameter which corresponds to the mean distance for the SIF.
It is worth noting that the method simultaneously estimate intensities of effect and scales of effect using a numerical likelihood maximisation procedure. This means that the mean distances of the SIFs are not given by the user but estimated for each landscape variable.
siland package is designed to be a user-friendly tool: data format are commonly used (data.frame and shape files from GIS data), model expression is simple (classical formula as lm expression), various outputs are very easily computed: summary of estimates and tests, graphical maps of effects and graphical checking of likelihood maximisation. Based on linear model, it includes numerous extensions of linear model: interactions, random effects, GLM (binomial and poisson). It is also possible to build multiannual model. We detailed hereafter some examples of simple use of the siland package. (More details about the method can be found in https://www.biorxiv.org/content/10.1101/692566v1).
siland package requires two data objects:
st_read
function of the package sf
. For instance if your GIS data (in shp format) are in the FILE file, named NAME (i.e.FILE contains NAME.shp, NAME.dbf,...), you can easily import your data in landdata using the command: landdata=st_read(dsn = "FILE",layer="NAME")
Here, we use the objects dataSiland and landSiland, included in siland package. You can load data using the data
command.
library(siland) data(dataSiland) data(landSiland)
dataSiland is a data.frame with columns: "obs", variable of interest and "x1" variable, a continuous local variable. The "Id" column indicates the identification number of the plots where observations were made. landSiland is an object of class sf where landscape is characterised by two features named L1 and L2. Landscape is described by polygons. For each polygon, L1 is equal to 1 if the polygon is of type L1, 0 otherwise (so is it for L2).
str(dataSiland) str(landSiland)
You can look at spatial configuration of the loaded data using the following commands :
library(ggplot2) library(sf) L1pol=st_geometry(landSiland[landSiland$L1==1,])#extract an sf object with only polygons of type L1 L2pol=st_geometry(landSiland[landSiland$L2==1,])#extract an sf object with only polygons of type L2 ggplot(landSiland)+ geom_sf(colour="grey",fill="white")+ geom_sf(data=L1pol,fill="red")+ geom_sf(data=L2pol,fill="blue")+ geom_point(data=dataSiland, aes(X,Y),col="green")
siland()
Estimations for landscape and local variables are based on discretization with cells (or pixels) of the landscape. This discretization is called a raster. Each unit (cell/pixel) of landscape variable has an influence which is maximal at its location and decreases as the distance increases. The decrease of this influence is described by a spatial influence function (SIF). For modelling a given landscape variable influence, we consider every cell where the landscape variable is distributed, and computed its influence by summing spatial influence of all cells. The spatial influence of a cell is determined by its intensity (negative or positive) and its SIF, which is a density function described by its mean distance (in siland package). The siland package allows to estimate the intensity of effect and the scale of effects i.e. here the mean distance of SIF.
Estimating effects of the local variable x1 and landscape variables L1 and L2 on the variable of interest obs can be done with the command siland()
(the argument wd
is presented below):
res=siland(obs~x1+L1+L2,land=landSiland,data=dataSiland,wd=30)
res
Here the parameters related to the landscape variable L1 are intensity L1
(-11.644) and the mean distance of its SIF, SIF.L1
(111.189). The vector res$paramSIF
gives for each landscape variable the mean distance of estimated SIFs. The data.frame res$landcontri
gives for each landscape variable, the cumulative spatial influence received by each observation (i.e. the sum of spatial influence received by each observation from all cells of the landscape). AICis AIC the of the model resB1 (Model: obs ~ x1 + L1 + L2) and
AIC(no landscape)is the AIC with only local variable (Model0: obs ~ x1).
(No landscape effect) p-value` is related to the test of landscape global effect significativity, i.e. H0= 'Model0: obs ~ x1' vs H1='Model: obs ~ x1 + L1 + L2' (Likelihood ratio test).
Parameter estimates and tests can be obtained using the command summary
:
summary(res)
Note that p.values are given conditionally to estimated parameters for SIFs.
The command plotsiland.sif
gives a representation of the estimated SIF for landscape variables. The shape of the SIF could be "exponential", "gaussian" or "uniform". The shape can be defined with the argument sif in the function siland()
and is exponential by default. The function plotsiland.sif()
allows to displays the estimated SIFs.
plotsiland.sif(res)
The command plotsiland.land
gives map of effects of the landscape.
By default, it gives the global contribution of all landscape variables:
plotsiland.land(x=res,land=landSiland,data=dataSiland)
One can obtain a map for a specific landscape variable by specifying its number with the argument var
:
plotsiland.land(x=res,land=landSiland,data=dataSiland,var=2)
As with all numerical maximization procedures, optimization problems may arise. The function siland.lik()
allows to point out possible problems of optimization.
likres=siland.lik(res,land= landSiland,data=dataSiland,varnames=c("L1","L2"))
likres
On this graphic, the -Log-likelihood is represented in function of buffer radii. The estimation is made by maximazing the likelihood i.e. by minimizing the -Log-likelihood. The orange line indicates the minimal value obtained during the estimation.
The black line represents the -loglikelihood for L1 when the buffer radius for L1 is set to given values in seqd argument (seqd=seq(2,200,length=10)
) and the other parameters are fixed to their estimation. The black dotted line indicates the value of L1 buffer size estimated during the global estimation procedure. The red continuous and dotted red lines simirlarly indicates the -loglikelihood and the value of buffer size estimated for L2. When minization correctly occurs, the minimal values of the profiled -loglikelihoods are equal and equal the minimal value of the global -Log-likelihood. This means that the minimums of continuous black and red lines are on the orange line (and that dotted black and red lines intersect the continuous black and red lines at their minimum, respectively). If it is not the case, the minimizing procedure has failed and it is necessary to proceed with a new estimation with different initialisation values. This is possible using the argument init
in function siland
. For instance, for starting estimation from buffer sizes of 20 and 150 for L1 and L2 landscape variables, respectively, use the command siland(obs~x1+L1+L2,land=landSiland,data=dataSiland,init=c(20,150))
.
The argument wd
indicates the mesh size used to construct spatial unit (cells/pixels) of landscape variable. In fact, landscape variable described as polygon are discretized on grid during the procedure estimation (like raster). The choice of wd is a tradeoff between computing precision and computing time (and memory size). The smallest wd
is, the better are the precision but the longer the computing time is (and the larger the required memory size is). It is worth to note that estimated parameters can be very sensitive to this mesh size. To obtain a reliable estimation, we recommand to ensure, after the estimation procedure, that wd
size is at least three times smaller than the smallest estimated SIF. If not, it is recommended to proceed with a new estimation with a smaller wd
size.
The argument maxD
corresponds to a radius. For each observation, only pixels with a distance less than maxD are used for estimation. For a given observation, pixels with a distance greater than maxD have in fact a very low intensity and can be ignored. To obtain a reliable estimation, we recommand to ensure, after the estimation procedure, that maxD
size is at least three times greater than the greatest estimated SIF.
The argument border
indicates whether an observation receives the spatial influence of cells belonging to the same plot (i.e. the plot where the observation was made). If border=T
, its receives influence only from the border of the plot (so no influence from the cells of the plot). By default, border = F
(all cells are considered). border
is a vector that defines if border
is True or False for each landscape variable. If border=T
(resp. border=F
) then the argument border
is set to True (resp. False) for all landscape variables.
The argument sif
defines the family function of SIF, i.e. the form of decrease for the landscape influence. The family can be exponential
, gaussian
, or uniform
. If influence is uniform, it implies that there is no decrease of influence and influence is uniform around each pixel of the raster. Note that all landscape variables have the same form.
Random effect can be included using the syntax (1| )
. Note that only local effect are concerned.
res2=siland(obs~x1+L1+L2+(1|Id),land=landSiland,data=dataSiland,wd=30)
summary(res2)
To consider distributions of the variable of interest that differ from Gaussian, use the option family
available in the siland()
function. Family can be "gaussian", "poisson" or "binomial" and the associated link function are identity, log and logit respectively.
Interaction between local and landscape variables using the syntax :
or *
(as well as localxlocal interaction). Interaction term modify the intensity of a landscape variable according to the local variable value. But the scale of effect of the landscape variable remains constant.
#Model with main and interaction effect res3=siland(obs~x1*L1+L2,land=landSiland,data=dataSiland,wd=30)
res3
res3.1=siland(obs~x1*L1+L2,land=landSiland,data=dataSiland,wd=30)
res3.1
It is possible to deal with data observed during several years, i.e. data associated with several landscapes with the function siland
.
The two data objects required are then :
a list of data.frames of located observations, one for each year. An important point is that data.frames column names have to be exactly the same and ordered in the same way.
a list of sf object containing a landscape description of each year.
Let us suppose we have two years of observations associated with two landscapes. Since the goal is only to show how to deal with such datasets, we take the same datastets for observations and for landscapes for the two years.
landSilandY1=landSiland landSilandY2=landSiland #landSilandY is a list with the landscape for each year landSilandY=list(year1=landSilandY1,year2=landSilandY2) dataSilandY1=dataSiland dataSilandY2=dataSiland dataSilandY1$year=factor("2018") dataSilandY2$year=factor("2019") head(dataSilandY1) head(dataSilandY2) dataSilandY=list(year1=dataSilandY1,year2=dataSilandY2) #dataSilandY is a list with the observed data for each year resY=siland(obs~year+x1+L1+L2, land = landSilandY,data=dataSilandY,wd=30) resY
resY
summary(resY)
It is possible to deal with multisite data, using the same procedure considering different sites instead of different years.
As for an object ot type GLM, functions AIC(), residuals() and fitted() are available. In fact, conditionnaly to the estimated buffers or SIFs, the fitted model is a GLM or a LMM (Linear Mixed Model) or a GLMM (Generalized Linear Mixed Model). So after an estimation, it is possible to analyse more precisely the estimated model with the object result stored in the output. Object result is an an object from glm() or lmer() or glmer() functions.
summary(res$result) BIC(res$result) fitted(res$result)[1:10] residuals(res$result)[1:10] class(res$result) class(res$result)
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.