[R-sig-Geo] convert "errorsarlm" to "lm"

Mon Oct 23 10:48:07 CEST 2017

On Mon, 23 Oct 2017, niv de malach wrote:

> I am not sure that I fully understand you... I Appologize for my ignorance

Even though piecewiseSEM uses a spatial example, I would argue that, 
unless we know with certainty that the multivariate observations are not 
in any way autocorrelated, we should assume that they are autocorrelated, 
and may also be coorrelated through missing covariates and/or 
inappropriate functional form. This will potentially affect the inferences 
we draw from correlations, as any unmodelled spatial autocorrelation will 
affect the count of "independent" observations, so that the effective df 
will be smaller than the df assumed if observations are independent. In 
path modelling (unlike simultaneous equation models, which feed back 
between equations), we are dealing with unidirectional causality, but 
inferences depend on assumptions about the number of independent 
observations. MRF fit a plausible spatially structured random effect (but 
by each equation separately - there are spatial econometrics approaches to 
SEM which model a multivariate error process).

As a speculation, you might be able to look at an ade4::multispati 
approach, packages ade4 or adespatial (using Moran eigenvectors instead of 
ML SAR - errorsarlm). See for example DOI: 10.1890/11-1183.1. I suspect 
that the union of sets of eigenvalues chosen for your lm() fits, added to 
the covariates, might account for the spatial patterning present in 
residuals (a PCNM-Moran eigenvector-spatial filtering approach).

>
> 1. Do you have some approach of transferring the object coefficients into
> "lm" assuming that they are fixed? (if so let me know how...)
>

If you read ?errorsalm, you see that tarX and tary are returned. They are 
(I - \lambda W) X and (I - \lambda W) y respectively. In ?bptest.sarlm, 
you see a sneaky and really untested example of their use - not 
recommended!!

> 2. I am not sure what do you mean by uncertainty (confidence intervals?
> significance tests? AIC\BIC scores?something totally different?). You
> suggested to use bayesian methods but I am not sure why. In this SEM type
> (local estimation approach) I can use the equation (before transformation)
> for inference (with all their "uncertainty" whatever that be). Putting them
> together as "SEM" is used mainly for testing whether there is a correlation
> between some variables and the residuals of equation describing other (non
> linked) variables . Using a special test this package test what is the
> probability that the model is "wrong" based on the "missing links".
> In general it would be better to get also an AIC\BIC for the total SEM
> model  (which is possible when regular lm objects are used) but my main
> goal in this analysis is parameter estimation rather than model selection
> (so the SEM is used only as a "sanity check" that
> there wasn't any important missing link)
>

Using the distribution of \lambda as something like N(lambda, se.lambda), 
sample many lambda values and re-generate tarX and tary. Re-fit your lm(), 
and examine the outcomes of piecewiseSEM for these possible values of 
\lambda (you have a joint situation with several equations, so treating 
each equation separately is not ideal - a PCNM approach would accommodate 
the shared spatial autocorrelation.

Assuming that issues of this kind only need technical fixes may lead to 
very erroneous inference. I'd guess that there is a viable paper in 
extending PCNM to path modelling, and that is a better way to go than an 
ad-hoc kludge.

Roger

> Thanks
>
> On Sun, Oct 22, 2017 at 10:23 PM, Roger Bivand <Roger.Bivand at nhh.no> wrote:
>
>> On Sun, 22 Oct 2017, niv de malach wrote:
>>
>> Regarding your question about MRF, I don't know what it is...
>>>
>>
>> see for example:
>>
>> http://www.fromthebottomoftheheap.net/2017/10/19/first-
>> steps-with-mrf-smooths/
>>
>> here using mgcv::gam() to fit a model with an MRF smooth. In effect,
>> errorsarlm() is fitting something similar,
>>
>>
>>> I know that the package piecewiseSEM accept the following object types:
>>>
>>> ‘lm‘, ‘rq‘, ‘glm‘, ‘glm.nb‘, ‘gls‘, ‘pgls‘, ‘merMod‘, ‘merModLmerTest‘,
>>> ‘lme‘, ‘glmmPQL‘, ‘glmmadmb‘, and ‘glmmTMB‘.
>>>
>>> so  "transformation" of errorsarlm to any of these object should solve the
>>> problem.
>>>
>>>
>> No! You have to pass through the uncertainty associated with fitting the
>> spatial regression coefficient. The same would apply with any other
>> MRF-based approach. It is possible to pretend that this coefficient is
>> fixed (I'll return to the example tomorrow). I don't think that any of the
>> listed object types come from fitting functions supporting MRF, though. You
>> would likely need to re-fit your piecewiseSEM over samples from the
>> coefficients, to check that the path outcomes were robust (just a
>> speculation). Any Bayesians who could help?
>>
>> Roger
>>
>>
>>> Below there is a code where a dataset is simulated and analyzed using 'lm'
>>> and 'errorsarlm'. Only 'lm' is accepted for
>>> SEM analysis.
>>>
>>> Thanks,
>>> Niv
>>>
>>>
>>> #################### creating a simulated dataset
>>> ##############################
>>>
>>> set.seed(1)
>>> sample=50 # sample size (number of points)
>>>
>>> # spatial location
>>> Lattitude=rnorm(sample, mean = 0, sd = 4)
>>> Longitude=rnorm(sample, mean = 0, sd = 4
>>>
>>> # creating exploratory variable
>>> <https://maps.google.com/?q=ry+variable&entry=gmail&source=g>s
>>> x1=rnorm(sample, mean = 30, sd = 2
>>> <https://maps.google.com/?q=30,+sd+%3D+2&entry=gmail&source=g>0)
>>> x2=rnorm(sample, mean = -20, sd = 10)
>>> x3=rnorm(sample, mean = 5, sd = 3)
>>> x4=rnorm(sample, mean = 25,
>>> <https://maps.google.com/?q=mean+%3D+25,&entry=gmail&source=g> sd = 5)
>>> x5=rnorm(sample, mean = -10, sd = 10)
>>> x6=rnorm(sample, mean = 50, sd = 30)
>>>
>>> # response variables
>>> y1=10*x1+5*x2+2*x3+30*Lattitude+rnorm(sample, mean = 0, sd = 20) # y1 is
>>> function of x1-x3 +location+random error
>>> y2=10*x4+5*x5+2*x6+10*Longitude+rnorm(sample, mean = 0, sd = 20) # y1 is
>>> function of x4-x5 +location+random error
>>>
>>> mydata=data.frame(x1,x2,x3,x4,x5,x6,y1,y2,Lattitude,Longitude)
>>>
>>> ################# building errorsarlm objects and lm objects
>>> ####################################
>>>
>>> library(spdep)
>>>
>>> # creating neighbors for errorsarlm function
>>> xy=SpatialPoints(as.matrix(mydata[,c(9,10)]))
>>> myneighbors=dnearneigh(as.matrix(mydata[,c(9,10)]), 0, 1000, row.names =
>>> NULL, longlat = TRUE)# 1000 kilometer distance
>>> myweight=nb2listw(myneighbors,style="W") #
>>>
>>> # two errorsarlm objects
>>> sar1=errorsarlm(y1~x1+x2+x3,data=mydata,myweight)
>>> sar2=errorsarlm(y2~x4+x5+x6,data=mydata,myweight)
>>>
>>> # two lm objects
>>> lm1=lm(y1~x1+x2+x3,data=mydata)
>>> lm2=lm(y2~x4+x5+x6,data=mydata)
>>>
>>> ############### trying to use sarlm and lm in "piecewiseSEM"
>>> ####################
>>>
>>> library(piecewiseSEM) # package for structural equation modeling
>>> SEM_LM=sem.fit(list(lm1,lm2),data=mydata) # the package accept lm objects
>>> (and random models such as lme)
>>> SEM_errorsarLM=sem.fit(list(sar1,sar2),data=mydata) # the package does
>>> not
>>> accept errprsarlm objects
>>>
>>> # Error in nobs.default(x) : no 'nobs' method is available #
>>>
>>>
>>>
>>>
>>> ᐧ
>>>
>>> On Fri, Oct 20, 2017 at 2:06 PM, Roger Bivand <Roger.Bivand at nhh.no>
>>> wrote:
>>>
>>> On Wed, 18 Oct 2017, niv de malach wrote:
>>>>
>>>> Hi all,
>>>>
>>>>> I need to convert an "errorsarlm" object into a simple "lm" object
>>>>> because
>>>>> the package I use ("piecewiseSEM") could not deal with "errorsarlm". Is
>>>>> there a way to transform the object?
>>>>>
>>>>>
>>>> Please give a reproducible example [1] showing what you are doing,
>>>> preferably with a built-in data set. In any case, an lm object made from
>>>> a
>>>> sarlm object will fix the spatial coefficient and loose its uncertainty.
>>>>
>>>> [1] https://stackoverflow.com/questions/5963269/how-to-make-a-
>>>> great-r-reproducible-example
>>>>
>>>> Can you use an model class accepted in piecewiseSEM that fits an MRF
>>>> random effect?
>>>>
>>>> Roger
>>>>
>>>>
>>>> Thanks,
>>>>> Niv
>>>>> ᐧ
>>>>>
>>>>>         [[alternative HTML version deleted]]
>>>>>
>>>>> _______________________________________________
>>>>> R-sig-Geo mailing list
>>>>> R-sig-Geo at r-project.org
>>>>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>>>>>
>>>>>
>>>> --
>>>> Roger Bivand
>>>> Department of Economics, Norwegian School of
>>>> <https://maps.google.com/?q=an+School+of&entry=gmail&source=g>Economics,
>>>> Helleveien 30, N-5045 Bergen, Norway.
>>>> voice: +47 55 95 93 55; e-mail: Roger.Bivand at nhh.no
>>>> Editor-in-Chief of The R Journal, https://journal.r-project.org/
>>>> index.html
>>>> http://orcid.org/0000-0003-2392-6140
>>>> https://scholar.google.no/citations?user=AWeghB0AAAAJ&hl=en
>>>>
>>>
>>>
>> --
>> Roger Bivand
>> Department of Economics, Norwegian School of Economics,
>> Helleveien 30, N-5045 Bergen, Norway.
>> voice: +47 55 95 93 55; e-mail: Roger.Bivand at nhh.no
>> Editor-in-Chief of The R Journal, https://journal.r-project.org/index.html
>> http://orcid.org/0000-0003-2392-6140
>> https://scholar.google.no/citations?user=AWeghB0AAAAJ&hl=en
>>
>
> ᐧ
>

-- 
Roger Bivand
Department of Economics, Norwegian School of Economics,
Helleveien 30, N-5045 Bergen, Norway.
voice: +47 55 95 93 55; e-mail: Roger.Bivand at nhh.no
Editor-in-Chief of The R Journal, https://journal.r-project.org/index.html
http://orcid.org/0000-0003-2392-6140
https://scholar.google.no/citations?user=AWeghB0AAAAJ&hl=en