Simple Ripley's CRS test for market point patters

classic Classic list List threaded Threaded
9 messages Options
Reply | Threaded
Open this post in threaded view
|

Simple Ripley's CRS test for market point patters

R-sig-geo mailing list
Dear R-Sig-Geo Members,     I"ve like to find any simple way for apply CRS test for market point patters, for this I try to create a script below:
#Packages require(spatstat)require(sp)

# Create some points that represents ant nests in UTMxp<-c(371278.588,371250.722,371272.618,371328.421,371349.974,371311.95,371296.265,371406.46,371411.551,371329.041,371338.081,371334.182,371333.756,371299.818,371254.374,371193.673,371172.836,371173.803,371153.73,371165.051,371140.417,371168.279,371166.367,371180.575,371132.664,371129.791,371232.919,371208.502,371289.462,371207.595,371219.008,371139.921,371133.215,371061.467,371053.69,371099.897,371108.782,371112.52,371114.241,371176.236,371159.185,371159.291,371158.552,370978.252,371120.03,371116.993)
yp<-c(8246507.94,8246493.176,8246465.974,8246464.413,8246403.465,8246386.098,8246432.144,8246394.827,8246366.201,8246337.626,8246311.125,8246300.039,8246299.594,8246298.072,8246379.351,8246431.998,8246423.913,8246423.476,8246431.658,8246418.226,8246400.161,8246396.891,8246394.225,8246400.391,8246370.244,8246367.019,8246311.075,8246255.174,8246255.085,8246226.514,8246215.847,8246337.316,8246330.197,8246311.197,8246304.183,8246239.282,8246239.887,8246241.678,8246240.361,8246167.364,8246171.581,8246171.803,8246169.807,8246293.57,8246183.194,8246189.926)
# Then I create the size of each nest - my covariate used as marked processarea<-c(117,30,4,341,15,160,35,280,108,168,63,143,2,48,182,42,88,56,27,156,288,45,49,234,72,270,91,40,304,56,35,4,56.7,9,4.6,105,133,135,23.92,190,12.9,15.2,192.78,104,255,24)

# Make a countour - only as exerciseW <- convexhull.xy(xp,yp)
#Create a ppp objectp_xy<-cbind(xp,yp)syn.ppp<-ppp(x=coordinates(p_xy)[,1],y=coordinates(p_xy)[,2],window=W, marks=area)syn.ppp <- as.ppp(syn.ppp)plot(syn.ppp, main=" ")
First, I've like to study CSR of market point process (my hypothesis is that different size have a same spatial distribution) when area >= 0, area < 25 and area >=25, area < 55, for this I make:
# Area 0-25env.formi1<-envelope(syn.ppp,nsim=99,fun=Kest, area >= 0, area < 25)plot(env.formi1,lwd=list(3,1,1,1), main="") 
# Area 25-55env.formi2<-envelope(syn.ppp,nsim=99,fun=Kest, area >=25, area < 55)plot(env.formi2,lwd=list(3,1,1,1), main="") 
My first problem is that I have the same plot in both conditions and I don't know why.
Second, if I try to estimate the market intensity observed pattern
est.int <- density(syn.ppp)est_xy <-  rmpoispp(est.int)plot(est_xy, main=" ")
My output is only my points position without marked area in my ppp object created.
My question is what is the problem with this Ripley's reduced second moment function approach? There are any way for study my point process when the area is a covariate of my point process?
Thanks in advanced
Alexandre


-- ======================================================================Alexandre dos SantosProteção Florestal IFMT - Instituto Federal de Educação, Ciência e Tecnologia de Mato GrossoCampus CáceresCaixa Postal 244Avenida dos Ramires, s/nBairro: Distrito Industrial Cáceres - MT                      CEP: 78.200-000Fone: (+55) 65 99686-6970 (VIVO) (+55) 65 3221-2674 (FIXO)e-mails:[hidden email]         [hidden email] Lattes: http://lattes.cnpq.br/1360403201088680 OrcID: orcid.org/0000-0001-8232-6722   -   ResearcherID: A-5790-2016Researchgate: www.researchgate.net/profile/Alexandre_Santos10                       LinkedIn: br.linkedin.com/in/alexandre-dos-santos-87961635Mendeley:www.mendeley.com/profiles/alexandre-dos-santos6/
======================================================================



        [[alternative HTML version deleted]]

_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Reply | Threaded
Open this post in threaded view
|

Re: Simple Ripley's CRS test for market point patters

Marcelino de la Cruz Rot
Dear Alexandre,

I think that the solution to your first problem would be something like
this:

#Area 0-25
env.formi1<-envelope(syn.ppp[syn.ppp$marks>=0 & syn.ppp$marks<25],
nsim=99,fun=Kest)
plot(env.formi1,lwd=list(3,1,1,1), main="")

Then, for your second problem, you should clarify what do you understand
by "market intensity observed pattern". In fact, rmpoispp()  simulates
random  multitype ppp's, i.e., it does not estimate parameters of
quantitatively marked point processes.

Maybe be you would consider interesting

Smooth.ppp(syn.ppp)

and

envelope(syn.ppp, nsim=99,fun=markcorr,
simulate=expression(rlabel(syn.ppp)))


Cheers,


Marcelino







El 22/07/2019 a las 23:36, Alexandre Santos via R-sig-Geo escribió:

> Dear R-Sig-Geo Members,     I"ve like to find any simple way for apply CRS test for market point patters, for this I try to create a script below:
> #Packages require(spatstat)require(sp)
>
> # Create some points that represents ant nests in UTMxp<-c(371278.588,371250.722,371272.618,371328.421,371349.974,371311.95,371296.265,371406.46,371411.551,371329.041,371338.081,371334.182,371333.756,371299.818,371254.374,371193.673,371172.836,371173.803,371153.73,371165.051,371140.417,371168.279,371166.367,371180.575,371132.664,371129.791,371232.919,371208.502,371289.462,371207.595,371219.008,371139.921,371133.215,371061.467,371053.69,371099.897,371108.782,371112.52,371114.241,371176.236,371159.185,371159.291,371158.552,370978.252,371120.03,371116.993)
> yp<-c(8246507.94,8246493.176,8246465.974,8246464.413,8246403.465,8246386.098,8246432.144,8246394.827,8246366.201,8246337.626,8246311.125,8246300.039,8246299.594,8246298.072,8246379.351,8246431.998,8246423.913,8246423.476,8246431.658,8246418.226,8246400.161,8246396.891,8246394.225,8246400.391,8246370.244,8246367.019,8246311.075,8246255.174,8246255.085,8246226.514,8246215.847,8246337.316,8246330.197,8246311.197,8246304.183,8246239.282,8246239.887,8246241.678,8246240.361,8246167.364,8246171.581,8246171.803,8246169.807,8246293.57,8246183.194,8246189.926)
> # Then I create the size of each nest - my covariate used as marked processarea<-c(117,30,4,341,15,160,35,280,108,168,63,143,2,48,182,42,88,56,27,156,288,45,49,234,72,270,91,40,304,56,35,4,56.7,9,4.6,105,133,135,23.92,190,12.9,15.2,192.78,104,255,24)
>
> # Make a countour - only as exerciseW <- convexhull.xy(xp,yp)
> #Create a ppp objectp_xy<-cbind(xp,yp)syn.ppp<-ppp(x=coordinates(p_xy)[,1],y=coordinates(p_xy)[,2],window=W, marks=area)syn.ppp <- as.ppp(syn.ppp)plot(syn.ppp, main=" ")
> First, I've like to study CSR of market point process (my hypothesis is that different size have a same spatial distribution) when area >= 0, area < 25 and area >=25, area < 55, for this I make:
> # Area 0-25env.formi1<-envelope(syn.ppp,nsim=99,fun=Kest, area >= 0, area < 25)plot(env.formi1,lwd=list(3,1,1,1), main="")
> # Area 25-55env.formi2<-envelope(syn.ppp,nsim=99,fun=Kest, area >=25, area < 55)plot(env.formi2,lwd=list(3,1,1,1), main="")
> My first problem is that I have the same plot in both conditions and I don't know why.
> Second, if I try to estimate the market intensity observed pattern
> est.int <- density(syn.ppp)est_xy <-  rmpoispp(est.int)plot(est_xy, main=" ")
> My output is only my points position without marked area in my ppp object created.
> My question is what is the problem with this Ripley's reduced second moment function approach? There are any way for study my point process when the area is a covariate of my point process?
> Thanks in advanced
> Alexandre
>
>
> -- ======================================================================Alexandre dos SantosProteção Florestal IFMT - Instituto Federal de Educação, Ciência e Tecnologia de Mato GrossoCampus CáceresCaixa Postal 244Avenida dos Ramires, s/nBairro: Distrito Industrial Cáceres - MT                      CEP: 78.200-000Fone: (+55) 65 99686-6970 (VIVO) (+55) 65 3221-2674 (FIXO)e-mails:[hidden email]         [hidden email] Lattes: http://lattes.cnpq.br/1360403201088680 OrcID: orcid.org/0000-0001-8232-6722   -   ResearcherID: A-5790-2016Researchgate: www.researchgate.net/profile/Alexandre_Santos10                       LinkedIn: br.linkedin.com/in/alexandre-dos-santos-87961635Mendeley:www.mendeley.com/profiles/alexandre-dos-santos6/
> ======================================================================
>
>
>
> [[alternative HTML version deleted]]
>
> _______________________________________________
> R-sig-Geo mailing list
> [hidden email]
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
> .
>

--
Marcelino de la Cruz Rot
Depto. de Biología y Geología
Física y Química Inorgánica
Universidad Rey Juan Carlos
Móstoles España

_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Reply | Threaded
Open this post in threaded view
|

DELDIR

PONS Frederic - CEREMA/DTerMed/DREC/SRILH
Dear all

Sometines we are using the deldir function for squeletisation.

Since the last review, deldir stops in some shapefile because polygons
are too narrow. In the version berfore, there was no problems.

We have readed the " Notes on error messages" and the problem of
anticlockwise order of triangle is listed.

In the trifnd R function , the code is
# Check that the vertices of the triangle listed in tau are
# in anticlockwise order.  (If they aren't then alles upgefucken
# ist; throw an error.)
call acchk(tau(1),tau(2),tau(3),anticl,x,y,ntot,eps)
if(!anticl) {
     call acchk(tau(3),tau(2),tau(1),anticl,x,y,ntot,eps)
     if(!anticl) {
         call fexit("Both vertex orderings are clockwise. See help for
deldir.")
     } else {
         ivtmp  = tau(3)
         tau(3) = tau(1)
         tau(1) = ivtmp
     }
}

We don't understand why do not order the bad triangles into the good
order. Perhaps, if this problem appears, the beginning of the deldir
function is not good.

If someone can explain us.

Best regards

*Frédéric Pons
*

        [[alternative HTML version deleted]]

_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Reply | Threaded
Open this post in threaded view
|

Re: DELDIR

Rolf Turner

This question should have been, in first instance, directed to me (the
maintainer of the deldir package), rather than to the r-sig-geo list.

On 24/07/19 1:47 AM, PONS Frederic - CEREMA/DTerMed/DREC/SRILH wrote:

> Dear all
>
> Sometines we are using the deldir function for squeletisation.

Had to Google this word. According to what I see it should be spelled
"squelettisation" (with a double "t").  The English word is
"skeletonisation", and this list is (for better or for worse) an English
language list.

> Since the last review, deldir stops in some shapefile because polygons
> are too narrow. In the version berfore, there was no problems.
>
> We have readed the " Notes on error messages" and the problem of
> anticlockwise order of triangle is listed.
>
> In the trifnd R function , the code is
> # Check that the vertices of the triangle listed in tau are
> # in anticlockwise order.  (If they aren't then alles upgefucken
> # ist; throw an error.)
> call acchk(tau(1),tau(2),tau(3),anticl,x,y,ntot,eps)
> if(!anticl) {
>       call acchk(tau(3),tau(2),tau(1),anticl,x,y,ntot,eps)
>       if(!anticl) {
>           call fexit("Both vertex orderings are clockwise. See help for
> deldir.")
>       } else {
>           ivtmp  = tau(3)
>           tau(3) = tau(1)
>           tau(1) = ivtmp
>       }
> }
>
> We don't understand why do not order the bad triangles into the good
> order. Perhaps, if this problem appears, the beginning of the deldir
> function is not good.
>
> If someone can explain us.

The error message you quote indicates that *both* orderings are bad so
something is toadally out of whack.  It would be rash (and arrogant) for
the software to fool around with the structure.  In such situations the
user should diagnose what is going wrong.

If you cannot diagnose the problem, please send me a reprex and I will
look into it.

cheers,

Rolf Turner

--
Honorary Research Fellow
Department of Statistics
University of Auckland
Phone: +64-9-373-7599 ext. 88276

_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Reply | Threaded
Open this post in threaded view
|

Re: [FORGED] Simple Ripley's CRS test for market point patters

Rolf Turner
In reply to this post by R-sig-geo mailing list

Thanks to Marcelino for chiming in on this.  I should have responded
earlier, but was "busy on a personal matter".

To add to what Marcelino said:

(1) Your post in HTML was horribly garbled and a struggle to read.
*PLEASE* set your email client so that it does *NOT* post in HTML when
posting to this list.

(2) A very minor point:  Your construction of "syn.ppp" was
unnecessarily complicated and convoluted.  You could simply do:

     syn.ppp  <- ppp(x=xp,y=yp,window=W,marks=area)

(3) You appear to be confused as to the distinction between "marks" and
"covariates".   These are superficially similar but conceptually
completely different.  What you are dealing with is *marks*.  There are
no covariates involved.  See [1], p. 147.

(4) Your calls to envelope() are mal-formed; the expression "area >= 0"
and "area < 25" will be taken as the values of "nrank" and ... who knows
what?  Did you not notice the warning messages you got about there being
something wrong with "nrank"?

You are being hopelessly naïve in expecting envelope() to interpret
"area >= 0" and "area < 25" in the way you want it to interpret them.
The spatstat package does not read minds.

Marcelino has shown you how to make a proper call.

(5) Since you are interested in categorising "area" into groups, rather
than being interested in the *numerical* value of area, you should do
the categorisation explicitly:

acat <- cut(area,breaks=c(-Inf,25,55,Inf),right=FALSE,
             labels=c("s","m","l")
# right=FALSE since you want "area < 25" rather than
# "area <= 25" --- although this makes no difference for the
# area values actually used.

syn.ppp <- ppp(x=xp,y=yp,marks=acat)

(6) It is not clear to me what your "research question" is.  Do you want
to test whether the intensities differ between the area categories?
Unless my thinking is hopelessly confused, this has nothing to do with
(or at least does not require the use of) the envelope() function:

f1 <- ppm(syn.ppp) # Same (constant) intensity for each area category.
f2 <- ppm(syn.ppp ~ marks) # Allows different (constant) intensity
                            # for each area category.
anova(f1,f2,test="Chi")

This test (not surprisingly) rejects the hypothesis that the intensities
are the same; p-value = 0.0020.

If this is not the hypothesis that you are interested in, please clarify
your thinking and your question, and get back to us.

cheers,

Rolf Turner

--
Honorary Research Fellow
Department of Statistics
University of Auckland
Phone: +64-9-373-7599 ext. 88276

[1] Spatial Point Patterns: Methodology and Applications with R
1st Edition, Adrian Baddeley, Ege Rubak, Rolf Turner
Chapman and Hall/CRC, 2015

P. S.  It is of course (!!!) highly recommended that you spend some time
reading the book cited above, if you are going to work in this area.

R. T.

On 23/07/19 9:36 AM, Alexandre Santos via R-sig-Geo wrote:

> Dear R-Sig-Geo Members,     I"ve like to find any simple way for apply CRS test for market point patters, for this I try to create a script below:
> #Packages require(spatstat)require(sp)
>
> # Create some points that represents ant nests in UTMxp<-c(371278.588,371250.722,371272.618,371328.421,371349.974,371311.95,371296.265,371406.46,371411.551,371329.041,371338.081,371334.182,371333.756,371299.818,371254.374,371193.673,371172.836,371173.803,371153.73,371165.051,371140.417,371168.279,371166.367,371180.575,371132.664,371129.791,371232.919,371208.502,371289.462,371207.595,371219.008,371139.921,371133.215,371061.467,371053.69,371099.897,371108.782,371112.52,371114.241,371176.236,371159.185,371159.291,371158.552,370978.252,371120.03,371116.993)
> yp<-c(8246507.94,8246493.176,8246465.974,8246464.413,8246403.465,8246386.098,8246432.144,8246394.827,8246366.201,8246337.626,8246311.125,8246300.039,8246299.594,8246298.072,8246379.351,8246431.998,8246423.913,8246423.476,8246431.658,8246418.226,8246400.161,8246396.891,8246394.225,8246400.391,8246370.244,8246367.019,8246311.075,8246255.174,8246255.085,8246226.514,8246215.847,8246337.316,8246330.197,8246311.197,8246304.183,8246239.282,8246239.887,8246241.678,8246240.361,8246167.364,8246171.581,8246171.803,8246169.807,8246293.57,8246183.194,8246189.926)
> # Then I create the size of each nest - my covariate used as marked processarea<-c(117,30,4,341,15,160,35,280,108,168,63,143,2,48,182,42,88,56,27,156,288,45,49,234,72,270,91,40,304,56,35,4,56.7,9,4.6,105,133,135,23.92,190,12.9,15.2,192.78,104,255,24)
>
> # Make a countour - only as exerciseW <- convexhull.xy(xp,yp)
> #Create a ppp objectp_xy<-cbind(xp,yp)syn.ppp<-ppp(x=coordinates(p_xy)[,1],y=coordinates(p_xy)[,2],window=W, marks=area)syn.ppp <- as.ppp(syn.ppp)plot(syn.ppp, main=" ")
> First, I've like to study CSR of market point process (my hypothesis is that different size have a same spatial distribution) when area >= 0, area < 25 and area >=25, area < 55, for this I make:
> # Area 0-25env.formi1<-envelope(syn.ppp,nsim=99,fun=Kest, area >= 0, area < 25)plot(env.formi1,lwd=list(3,1,1,1), main="")
> # Area 25-55env.formi2<-envelope(syn.ppp,nsim=99,fun=Kest, area >=25, area < 55)plot(env.formi2,lwd=list(3,1,1,1), main="")
> My first problem is that I have the same plot in both conditions and I don't know why.
> Second, if I try to estimate the market intensity observed pattern
> est.int <- density(syn.ppp)est_xy <-  rmpoispp(est.int)plot(est_xy, main=" ")
> My output is only my points position without marked area in my ppp object created.
> My question is what is the problem with this Ripley's reduced second moment function approach? There are any way for study my point process when the area is a covariate of my point process?
> Thanks in advanced
> Alexandre

_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Reply | Threaded
Open this post in threaded view
|

Re: Simple Ripley's CRS test for market point patters

R-sig-geo mailing list
Thanks for your help Marcelino e for the careful explanation Rolf Turner
and so sorry about my post script configuration, I expected that I
solved that in my new post.

First my variable area is a marked point (attribute or auxiliary
information about my point process - page 7 and not a spatial covariate,
effect in the outcome of my experimental area - page 50). Based in this
information, my hypophyses is that the *size of ant nests* a cause of
ecological intraspecific competition for resources (such as food and
territory) *have different patterns of spatial distribution*, for this:

#Packages
require(spatstat)
require(sp)

# Create some points that represents ant nests

xp<-c(371278.588,371250.722,371272.618,371328.421,371349.974,
371311.95,371296.265,371406.46,371411.551,371329.041,371338.081,
371334.182,371333.756,371299.818,371254.374,371193.673,371172.836,
371173.803,371153.73,371165.051,371140.417,371168.279,371166.367,
371180.575,371132.664,371129.791,371232.919,371208.502,371289.462,
371207.595,371219.008,371139.921,371133.215,371061.467,371053.69,
371099.897,371108.782,371112.52,371114.241,371176.236,371159.185,
371159.291,371158.552,370978.252,371120.03,371116.993)

yp<-c(8246507.94,8246493.176,8246465.974,8246464.413,8246403.465,
8246386.098,8246432.144,8246394.827,8246366.201,8246337.626,
8246311.125,8246300.039,8246299.594,8246298.072,8246379.351,
8246431.998,8246423.913,8246423.476,8246431.658,8246418.226,
8246400.161,8246396.891,8246394.225,8246400.391,8246370.244,
8246367.019,8246311.075,8246255.174,8246255.085,8246226.514,
8246215.847,8246337.316,8246330.197,8246311.197,8246304.183,
8246239.282,8246239.887,8246241.678,8246240.361,8246167.364,
8246171.581,8246171.803,8246169.807,8246293.57,8246183.194,8246189.926)

#Now I have the size of each nest (marked process)

area<-c(117,30,4,341,15,160,35,280,108,168,63,143,2,48,182,42,
88,56,27,156,288,45,49,234,72,270,91,40,304,56,35,4,56.7,9,4.6,
105,133,135,23.92,190,12.9,15.2,192.78,104,255,24)

# Make a countour for the window creation
W <- convexhull.xy(xp,yp)


# Class of nests size - large, medium and small
acat <- cut(area,breaks=c(-Inf,25,55,Inf),right=FALSE,
labels=c("s","m","l"))

#Create a ppp object

syn.ppp <- ppp(x=xp,y=yp,window=W, marks=acat)

# Test intensity hypothesis

f1 <- ppm(syn.ppp) # Same (constant) intensity for each area category.
f2 <- ppm(syn.ppp ~ marks) # Allows different (constant) intensity for
each area category.
anova(f1,f2,test="Chi")

#0.002015 ** OK, the hypothesis that the intensities are the same was
reject, but intensities is not the question.

Based in my initial hypothesis, I've like to show envelopes and observed
values of the use of some function for the three point patters (large,
medium and small ant nest size)under CSR. And for this I try to:

kS<-envelope(syn.ppp[syn.ppp$acat=="s"], nsim=99,fun=Kest)

plot(kS,lwd=list(3,1,1,1), main="")

kM<-envelope(syn.ppp[syn.ppp$acat=="m"], nsim=99,fun=Kest)

plot(kM,lwd=list(3,1,1,1), main="")

kL<-envelope(syn.ppp[syn.ppp$acat=="l"], nsim=99,fun=Kest)

plot(kL,lwd=list(3,1,1,1), main="")

But doesn't work yet. My approach now sounds correct?

Thanks in advanced,

Alexandre


--
======================================================================
Alexandre dos Santos
Prote????o Florestal
IFMT - Instituto Federal de Educa????o, Ci??ncia e Tecnologia de Mato Grosso
Campus C??ceres
Caixa Postal 244
Avenida dos Ramires, s/n
Bairro: Distrito Industrial
C??ceres - MT                      CEP: 78.200-000
Fone: (+55) 65 99686-6970 (VIVO) (+55) 65 3221-2674 (FIXO)

         [hidden email]
Lattes: http://lattes.cnpq.br/1360403201088680
OrcID: orcid.org/0000-0001-8232-6722
Researchgate: www.researchgate.net/profile/Alexandre_Santos10
LinkedIn: br.linkedin.com/in/alexandre-dos-santos-87961635
Mendeley:www.mendeley.com/profiles/alexandre-dos-santos6/
======================================================================

Em 23/07/2019 22:00, Rolf Turner escreveu:

>
> Thanks to Marcelino for chiming in on this.?? I should have responded
> earlier, but was "busy on a personal matter".
>
> To add to what Marcelino said:
>
> (1) Your post in HTML was horribly garbled and a struggle to read.
> *PLEASE* set your email client so that it does *NOT* post in HTML when
> posting to this list.
>
> (2) A very minor point:?? Your construction of "syn.ppp" was
> unnecessarily complicated and convoluted.?? You could simply do:
>
> ?????? syn.ppp?? <- ppp(x=xp,y=yp,window=W,marks=area)
>
> (3) You appear to be confused as to the distinction between "marks"
> and "covariates".???? These are superficially similar but conceptually
> completely different.?? What you are dealing with is *marks*.?? There are
> no covariates involved.?? See [1], p. 147.
>
> (4) Your calls to envelope() are mal-formed; the expression "area >=
> 0" and "area < 25" will be taken as the values of "nrank" and ... who
> knows what??? Did you not notice the warning messages you got about
> there being something wrong with "nrank"?
>
> You are being hopelessly na??ve in expecting envelope() to interpret
> "area >= 0" and "area < 25" in the way you want it to interpret them.
> The spatstat package does not read minds.
>
> Marcelino has shown you how to make a proper call.
>
> (5) Since you are interested in categorising "area" into groups,
> rather than being interested in the *numerical* value of area, you
> should do the categorisation explicitly:
>
> acat <- cut(area,breaks=c(-Inf,25,55,Inf),right=FALSE,
> ?????????????????????? labels=c("s","m","l")
> # right=FALSE since you want "area < 25" rather than
> # "area <= 25" --- although this makes no difference for the
> # area values actually used.
>
> syn.ppp <- ppp(x=xp,y=yp,marks=acat)
>
> (6) It is not clear to me what your "research question" is.?? Do you
> want to test whether the intensities differ between the area
> categories? Unless my thinking is hopelessly confused, this has
> nothing to do with (or at least does not require the use of) the
> envelope() function:
>
> f1 <- ppm(syn.ppp) # Same (constant) intensity for each area category.
> f2 <- ppm(syn.ppp ~ marks) # Allows different (constant) intensity
> ???????????????????????????????????????????????????? # for each area category.
> anova(f1,f2,test="Chi")
>
> This test (not surprisingly) rejects the hypothesis that the
> intensities are the same; p-value = 0.0020.
>
> If this is not the hypothesis that you are interested in, please
> clarify your thinking and your question, and get back to us.
>
> cheers,
>
> Rolf Turner
>

        [[alternative HTML version deleted]]

_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Reply | Threaded
Open this post in threaded view
|

Re: Simple Ripley's CRS test for market point patters

Rolf Turner

On 25/07/19 12:50 PM, ASANTOS wrote:

> Thanks for your help Marcelino e for the careful explanation Rolf Turner
> and so sorry about my post script configuration, I expected that I
> solved that in my new post.
>
> First my variable area is a marked point (attribute or auxiliary
> information about my point process - page 7 and not a spatial covariate,
> effect in the outcome of my experimental area - page 50). Based in this
> information, my hypophyses is that the *size of ant nests* a cause of
> ecological intraspecific competition for resources (such as food and
> territory) *have different patterns of spatial distribution*, for this:
>
> #Packages
> require(spatstat)
> require(sp)
>
> # Create some points that represents ant nests
>
> xp<-c(371278.588,371250.722,371272.618,371328.421,371349.974,
> 371311.95,371296.265,371406.46,371411.551,371329.041,371338.081,
> 371334.182,371333.756,371299.818,371254.374,371193.673,371172.836,
> 371173.803,371153.73,371165.051,371140.417,371168.279,371166.367,
> 371180.575,371132.664,371129.791,371232.919,371208.502,371289.462,
> 371207.595,371219.008,371139.921,371133.215,371061.467,371053.69,
> 371099.897,371108.782,371112.52,371114.241,371176.236,371159.185,
> 371159.291,371158.552,370978.252,371120.03,371116.993)
>
> yp<-c(8246507.94,8246493.176,8246465.974,8246464.413,8246403.465,
> 8246386.098,8246432.144,8246394.827,8246366.201,8246337.626,
> 8246311.125,8246300.039,8246299.594,8246298.072,8246379.351,
> 8246431.998,8246423.913,8246423.476,8246431.658,8246418.226,
> 8246400.161,8246396.891,8246394.225,8246400.391,8246370.244,
> 8246367.019,8246311.075,8246255.174,8246255.085,8246226.514,
> 8246215.847,8246337.316,8246330.197,8246311.197,8246304.183,
> 8246239.282,8246239.887,8246241.678,8246240.361,8246167.364,
> 8246171.581,8246171.803,8246169.807,8246293.57,8246183.194,8246189.926)
>
> #Now I have the size of each nest (marked process)
>
> area<-c(117,30,4,341,15,160,35,280,108,168,63,143,2,48,182,42,
> 88,56,27,156,288,45,49,234,72,270,91,40,304,56,35,4,56.7,9,4.6,
> 105,133,135,23.92,190,12.9,15.2,192.78,104,255,24)
>
> # Make a countour for the window creation
> W <- convexhull.xy(xp,yp)
>
>
> # Class of nests size - large, medium and small
> acat <- cut(area,breaks=c(-Inf,25,55,Inf),right=FALSE,
> labels=c("s","m","l"))
>
> #Create a ppp object
>
> syn.ppp <- ppp(x=xp,y=yp,window=W, marks=acat)
>
> # Test intensity hypothesis
>
> f1 <- ppm(syn.ppp) # Same (constant) intensity for each area category.
> f2 <- ppm(syn.ppp ~ marks) # Allows different (constant) intensity for
> each area category.
> anova(f1,f2,test="Chi")
>
> #0.002015 ** OK, the hypothesis that the intensities are the same was
> reject, but intensities is not the question.
>
> Based in my initial hypothesis, I've like to show envelopes and observed
> values of the use of some function for the three point patters (large,
> medium and small ant nest size)under CSR. And for this I try to:
>
> kS<-envelope(syn.ppp[syn.ppp$acat=="s"], nsim=99,fun=Kest)
>
> plot(kS,lwd=list(3,1,1,1), main="")
>
> kM<-envelope(syn.ppp[syn.ppp$acat=="m"], nsim=99,fun=Kest)
>
> plot(kM,lwd=list(3,1,1,1), main="")
>
> kL<-envelope(syn.ppp[syn.ppp$acat=="l"], nsim=99,fun=Kest)
>
> plot(kL,lwd=list(3,1,1,1), main="")
>
> But doesn't work yet. My approach now sounds correct?
>
> Thanks in advanced,

The formatting of your latest email was a great improvement and is now
admirably legible.

I am still not completely clear what your research question is, or what
hypotheses you wish to test.

Let me preface my remarks by saying that I have so far gone along with
your inclination to convert the numeric marks (area) of your pattern to
categorical (small, medium, large) marks.  I believe that this may be a
suboptimal approach.  (As someone said in an old post to R-help, on an
entirely  different topic, 'There is a reason that the speedometer in
your car doesn't just read "slow" and "fast".' :-) )

On the other hand there is not a lot on offer in the way of techniques
for handling numeric marks, so the "discretisation" approach may be the
best that one can do.

Now to proceed with trying to answer your question(s).  You say:

> my hypophyses is that the size of ant nests a cause of ecological
> intraspecific competition for resources (such as food and territory)
> have different patterns of spatial distribution

That's a (rather vague) "*alternative*" hypothesis.  The corresponding
null  hypothesis is that the three categories have the *same* spatial
distribution.   We reject that null hypothesis, since we reject the
hypothesis that they have the same intensities.

Now what?

You say:

> Based in my initial hypothesis, I've like to show envelopes and observed
> values of the use of some function for the three point patters (large,
> medium and small ant nest size) under CSR.

This is a bit incoherent and a demonstrates somewhat confused thinking.
  As the signature file of a frequent poster to R-help says "Tell me
what you want to do, not how you want to do it."

Perhaps the following will help to clarify.  For a categorical marked
point pattern (which is what we are currently dealing with) there are
infinitely many possibilities for the joint distribution of the patterns
corresponding to the marks.

The simplest of these is CSR (which means that the marks are essentially
irrelevant).  We have already rejected this hypothesis.

The next simplest is what is designated in [1; p. 584] as "CSRI"
(complete spatial randomness with independence).   This means that each
of the patterns is "individually CSR" (Poisson process, constant
intensity) and that the patterns are independent.

After that, the possibilities explode astronomically --- different
(non-Poisson) structures for each of the patterns, different dependence
structures, ....  The sky's the limit.

One "reasonably simple" possibility amongst this plethora of
possibilities is the multitype Strauss point process model, which could
be fitted (with some effort) using ppm() and profilepl().

However before we open that can of worms, we can test "CSRI":

fit <- ppm(syn.ppp ~ marks)
set.seed(42)
E <- envelope(fit,savefuns=TRUE)
dclf.test(E)
mad.test(E)

The dclf test gives a p-value of 0.14; the mad test gives a p-value of
0.18.  Thus there appears to be no evidence against the null hypothesis
of CSRI.

Therefore "CSRI" would appear to be an adequate/plausible model for
these data.

There may be more that needs to be done, but this is a start.

I would appreciate comments from Marcelino (and of course from Adrian
and Ege) especially if any of them notice something incorrect in my advice.

cheers,

Rolf

[1] Spatial Point Patterns: Methodology and Applications with R
1st Edition, Adrian Baddeley, Ege Rubak, Rolf Turner
Chapman and Hall/CRC, 2015

P. S.  You said of your efforts to continue your analysis: "But doesn't
work yet."  This is a very vague and unhelpful assertion!  One thing
that is wrong with your efforts is that in,  e.g.:

     kS<-envelope(syn.ppp[syn.ppp$acat=="s"], nsim=99,fun=Kest)

the object syn.ppp has no component named "acat"!!!  (You really need to
gain some understanding of how R works!)

You *could* say

    kS <- envelope(syn.ppp[syn.ppp$marks=="s"])

or (better)

     kS <- envelope(syn.ppp[marks(syn.ppp)=="s"])

to get the sort of result that you had hoped for.  Note that specifying
the arguments nsim=99 and fun=Kest, while harmless, is redundant (and
therefore tends to obfuscate).  These are the *default* values!!!  Read
the help, and learn to use R efficiently!

R.

--
Honorary Research Fellow
Department of Statistics
University of Auckland
Phone: +64-9-373-7599 ext. 88276

_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Reply | Threaded
Open this post in threaded view
|

Re: Simple Ripley's CRS test for market point patters --- addendum.

Rolf Turner

I have realised that there were a couple of oversights in my previous
posting on this issue.  One is a bit subtle; the other is a bit of a
blunder on my part.

First the "subtle" one. The test that I proposed for CSRI is a test done
using the estimated parameters of the proposed model to generate
realisations of data sets under the null hypothesis.  Such tests tend to
be conservative.  (See section 10.6.3, p. 388 ff., in [1].)

In the current instance (testing for CSRI) the conservatism can be
overcome by simulating data conditional on the numbers of points of each
type in the "real" data.  This can be done here via:

foo <- function(W){
s <- runifpoint(10,win=W)
m <- runifpoint(9,win=W)
l <- runifpoint(27,win=W)
superimpose(s=s,m=m,l=l)
}
simex <- expression(foo(W))

and then

set.seed(42)
E <- envelope(syn.ppp,simulate=simex,savefuns=TRUE)
dtst <- dclf.test(E)
mtst <- mad.test(E)

This gives p-values of 0.06 from the dclf test and 0.09 from the mad
test.  Thus there appears to be some slight evidence against the null
hypothesis.  (Evidence at the 0.10 significance level.)

That this should be so is *OBVIOUS* (!!!) if we turn to the unsubtle
point that I overlooked.  It is clear that the pattern of ants' nests
cannot be truly a realisation of a Poisson process since there must be a
bit of a "hard core" effect.  Two ants' nests cannot overlap.  Thus if
we approximate the shape of each nest by a disc, points i and j must be
a distance of at least r_i + r_j from each other, where r_i =
sqrt(area_i/pi), and similarly for r_j.

However I note that the data provided seem to violate this principle in
several instances.  E.g. points 41 and 42 are a distance of only 0.2460
metres apart but areas 41 and 42 are 12.9 and 15.2 square metres,
yielding putative radii of 3.5917 and 3.8987 metres, whence the closest
these points could possibly be (under the "disc-shaped assumption") is
7.4904 metres, far larger than 0.2460.   So something is a bit out of
whack here.  Perhaps these are made-up ("synthetic") data and the
process of making up the data did not take account of the minimum
distance constraint.

How to incorporate the "hard core" aspect of your (real?) data into the
modelling exercise, and what the impact of it is upon your research
question(s), is unclear to me and is likely to be complicated.

cheers,

Rolf

--
Honorary Research Fellow
Department of Statistics
University of Auckland
Phone: +64-9-373-7599 ext. 88276

[1] Spatial Point Patterns: Methodology and Applications with R
1st Edition, Adrian Baddeley, Ege Rubak, Rolf Turner
Chapman and Hall/CRC, 2015

_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Reply | Threaded
Open this post in threaded view
|

Re: Simple Ripley's CRS test for market point patters --- addendum.

R-sig-geo mailing list
Dr. Rolf Turner,

 ???? Thanks again and your explanations as a personal spatial course for
me :) You are very clear, but there are some questions yet. First, these
data set is not artificial, I have real coordinates and area of nests of
leaf-cutting ants in my example (in my complete data set I have more
than 3,000 coordinates and areas). The problem of overlap is because we
represent a ants nests area as a circle (for an easy comprehension and
modeling), but in the real world this areas was elliptical and some
times we have an eccentricity close to 1 (in this case I have two nests
close but not in overlap condition in the field). Based on our comments,
I used the CSRI test, but first I try to explore more the question about
intensities and size categories (s, m, and l), for this,could I test
intensities by pair-wise comparison? For the last, do you have any
solution for my overlap problem, if my objective is the study of my
problem as a realization of a Poisson process?

Please see my code below:

#Packages
require(spatstat)
require(sp)

# Create some points that represent ant nests

xp<-c(371278.588,371250.722,371272.618,371328.421,371349.974,
371311.95,371296.265,371406.46,371411.551,371329.041,371338.081,
371334.182,371333.756,371299.818,371254.374,371193.673,371172.836,
371173.803,371153.73,371165.051,371140.417,371168.279,371166.367,
371180.575,371132.664,371129.791,371232.919,371208.502,371289.462,
371207.595,371219.008,371139.921,371133.215,371061.467,371053.69,
371099.897,371108.782,371112.52,371114.241,371176.236,371159.185,
371159.291,371158.552,370978.252,371120.03,371116.993)

yp<-c(8246507.94,8246493.176,8246465.974,8246464.413,8246403.465,
8246386.098,8246432.144,8246394.827,8246366.201,8246337.626,
8246311.125,8246300.039,8246299.594,8246298.072,8246379.351,
8246431.998,8246423.913,8246423.476,8246431.658,8246418.226,
8246400.161,8246396.891,8246394.225,8246400.391,8246370.244,
8246367.019,8246311.075,8246255.174,8246255.085,8246226.514,
8246215.847,8246337.316,8246330.197,8246311.197,8246304.183,
8246239.282,8246239.887,8246241.678,8246240.361,8246167.364,
8246171.581,8246171.803,8246169.807,8246293.57,8246183.194,8246189.926)

#Now I have the size of each nest (marked process)

area<-c(117,30,4,341,15,160,35,280,108,168,63,143,2,48,182,42,
88,56,27,156,288,45,49,234,72,270,91,40,304,56,35,4,56.7,9,4.6,
105,133,135,23.92,190,12.9,15.2,192.78,104,255,24)

# Make a contour for the window creation
W <- convexhull.xy(xp,yp)


# Class of nests size - large, medium and small
acat <- cut(area,breaks=c(-Inf,25,55,Inf),right=FALSE,
labels=c("s","m","l"))

#Create a ppp object

syn.ppp <- ppp(x=xp,y=yp,window=W, marks=acat)

#Test initial hypothesis
f1 <- ppm(syn.ppp) # Same (constant) intensity for each area category.
f2 <- ppm(syn.ppp ~ marks) # Allows different (constant) intensity for
each area category.
anova(f1,f2,test="Chi") #The test?? rejects the hypothesis that the
intensities are the same - p =0.002015 **

*#First question, could I test intensities by pair-wise comparison?*

# Small vs medium sizes
acat2<-acat
levels(acat2)
levels(acat2)[1]<-"s_m"
levels(acat2)[2]<-"s_m"
levels(acat2)
syn.ppp2 <- ppp(x=xp,y=yp,window=W, marks=acat2)
f3 <- ppm(syn.ppp2 ~ marks)
f4 <- ppm(syn.ppp2)
anova(f4,f3,test="Chi") #Intensities of small and medium size are the
same - p=0.237

# Medium vs large sizes
acat3<-acat
levels(acat3)
levels(acat3)[2]<-"m_l"
levels(acat3)[3]<-"m_l"
levels(acat3)
syn.ppp3 <- ppp(x=xp,y=yp,window=W, marks=acat3)
f5 <- ppm(syn.ppp3 ~ marks)
f6 <- ppm(syn.ppp3)
anova(f5,f6,test="Chi") #Intensities of medium and large sizes are the
different - p=7.827e-05 ***

Finally small and medium sizes has the same intensity but is different
of large ants size!!! With this condition I will test the CSRI:

*#Now testing CSRI by simulation*

foo <- function(W){
s_m <- runifpoint(length(area<55),win=W)
l <- runifpoint(length(area>=55),win=W)
superimpose(s_m=s_m,l=l)
}
simex <- expression(foo(W))

#and then

set.seed(12)
E <- envelope(syn.ppp2,simulate=simex,nsim=999, savefuns=TRUE)
plot(E)
dtst <- dclf.test(E)
plot(dtst)
mtst <- mad.test(E)
plot(mtst)
# now gives p-value = 0.003 for dtst and p-value = 0.002 for mtst

Thanks in advanced,

Alexandre

--
======================================================================
Alexandre dos Santos
Prote????o Florestal
IFMT - Instituto Federal de Educa????o, Ci??ncia e Tecnologia de Mato Grosso
Campus C??ceres
Caixa Postal 244
Avenida dos Ramires, s/n
Bairro: Distrito Industrial
C??ceres - MT                      CEP: 78.200-000
Fone: (+55) 65 99686-6970 (VIVO) (+55) 65 3221-2674 (FIXO)

         [hidden email]
Lattes: http://lattes.cnpq.br/1360403201088680
OrcID: orcid.org/0000-0001-8232-6722
Researchgate: www.researchgate.net/profile/Alexandre_Santos10
LinkedIn: br.linkedin.com/in/alexandre-dos-santos-87961635
Mendeley:www.mendeley.com/profiles/alexandre-dos-santos6/
======================================================================

Em 26/07/2019 17:29, Rolf Turner escreveu:

>
> I have realised that there were a couple of oversights in my previous
> posting on this issue.?? One is a bit subtle; the other is a bit of a
> blunder on my part.
>
> First the "subtle" one. The test that I proposed for CSRI is a test
> done using the estimated parameters of the proposed model to generate
> realisations of data sets under the null hypothesis. Such tests tend
> to be conservative.?? (See section 10.6.3, p. 388 ff., in [1].)
>
> In the current instance (testing for CSRI) the conservatism can be
> overcome by simulating data conditional on the numbers of points of
> each type in the "real" data.?? This can be done here via:
>
> foo <- function(W){
> s <- runifpoint(10,win=W)
> m <- runifpoint(9,win=W)
> l <- runifpoint(27,win=W)
> superimpose(s=s,m=m,l=l)
> }
> simex <- expression(foo(W))
>
> and then
>
> set.seed(42)
> E <- envelope(syn.ppp,simulate=simex,savefuns=TRUE)
> dtst <- dclf.test(E)
> mtst <- mad.test(E)
>
> This gives p-values of 0.06 from the dclf test and 0.09 from the mad
> test.?? Thus there appears to be some slight evidence against the null
> hypothesis.?? (Evidence at the 0.10 significance level.)
>
> That this should be so is *OBVIOUS* (!!!) if we turn to the unsubtle
> point that I overlooked.?? It is clear that the pattern of ants' nests
> cannot be truly a realisation of a Poisson process since there must be
> a bit of a "hard core" effect.?? Two ants' nests cannot overlap.?? Thus
> if we approximate the shape of each nest by a disc, points i and j
> must be a distance of at least r_i + r_j from each other, where r_i =
> sqrt(area_i/pi), and similarly for r_j.
>
> However I note that the data provided seem to violate this principle
> in several instances.?? E.g. points 41 and 42 are a distance of only
> 0.2460 metres apart but areas 41 and 42 are 12.9 and 15.2 square
> metres, yielding putative radii of 3.5917 and 3.8987 metres, whence
> the closest
> these points could possibly be (under the "disc-shaped assumption") is
> 7.4904 metres, far larger than 0.2460.???? So something is a bit out of
> whack here.?? Perhaps these are made-up ("synthetic") data and the
> process of making up the data did not take account of the minimum
> distance constraint.
>
> How to incorporate the "hard core" aspect of your (real?) data into
> the modelling exercise, and what the impact of it is upon your
> research question(s), is unclear to me and is likely to be complicated.
>
> cheers,
>
> Rolf
>

        [[alternative HTML version deleted]]

_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo