Overlap entre deux ellipses représentant la niche dans ade4

[Nicolas Romillac]

Bonjour à tous,
j'utilise la fonction 'niche' de ade4 pour dessiner les niches écologiques de plusieurs espèces (calculant l'Outlying Mean Index, OMI). J'obtiens un graphique avec des ellipses représentant la position des niches des espèces dans l'espace de mes variables environnementales.
Est-ce que quelqu'un saurait comment faire pour calculer l'overlap entre deux ellipse? ou a minima, obtenir les coordonnées des ellipses?

Ci-dessous un exemple simplifié:

Code : Tout sélectionner

#dataset environnement
a<-c(0.77,0.18,0.25,1.26)
b<-c(0.63,2.91,0.29,1.04)
env<-cbind(a,b)
dudi1 <- dudi.pca(env, scale = TRUE, scan = FALSE, nf = 3)

#dataset abondances d'espèces
sp1<-c(0,33,30,22)
sp2<-c(74,40,27,94)
sp<-cbind.data.frame(sp1,sp2)

#calcul de la niche
library(ade4)
nic1 <- niche(dudi1, sp, scann = FALSE)
plot(nic1)

#alternativement, pour représenter juste les ellipses
s.distri(dfxy = nic1$ls, dfdistri = sp,cstar=0,label=names(sp),clabel=1,cpoint=0,ylim=c(-3,3))


(je ne sais pas comment faire pour insérer l'image produite)

Bonne journée à tous

[Logez Maxime]

Bonjour,

Il te faut décortiquer le code des fonctions de ade4 :

Code : Tout sélectionner

util.ellipse <- function (mx, my, vx, cxy, vy, coeff)
{
    lig <- 1000
    epsi <- 1e-10
    x <- 0
    y <- 0
    if (vx < 0)
        vx <- 0
    if (vy < 0)
        vy <- 0
    if (vx == 0 && vy == 0)
        return(NULL)
    delta <- (vx - vy) * (vx - vy) + 4 * cxy * cxy
    delta <- sqrt(delta)
    l1 <- (vx + vy + delta)/2
    l2 <- vx + vy - l1
    if (l1 < 0)
        l1 <- 0
    if (l2 < 0)
        l2 <- 0
    l1 <- sqrt(l1)
    l2 <- sqrt(l2)
    test <- 0
    if (vx == 0) {
        a0 <- 0
        b0 <- 1
        test <- 1
    }
    if ((vy == 0) && (test == 0)) {
        a0 <- 1
        b0 <- 0
        test <- 1
    }
    if (((abs(cxy)) < epsi) && (test == 0)) {
        if (vx > vy) {
            a0 <- 1
            b0 <- 0
        }
        else {
            a0 <- 0
            b0 <- 1
        }
        test <- 1
    }
    if (test == 0) {
        a0 <- 1
        b0 <- (l1 * l1 - vx)/cxy
        norm <- sqrt(a0 * a0 + b0 * b0)
        a0 <- a0/norm
        b0 <- b0/norm
    }
    a1 <- 2 * pi/lig
    c11 <- coeff * a0 * l1
    c12 <- (-coeff) * b0 * l2
    c21 <- coeff * b0 * l1
    c22 <- coeff * a0 * l2
    angle <- 0
    for (i in 1:lig) {
        cosinus <- cos(angle)
        sinus <- sin(angle)
        x[i] <- mx + c11 * cosinus + c12 * sinus
        y[i] <- my + c21 * cosinus + c22 * sinus
        angle <- angle + a1
    }
    return(list(x = x, y = y, seg1 = c(mx + c11, my + c21, mx -
        c11, my - c21), seg2 = c(mx + c12, my + c22, mx - c12,
        my - c22)))
}

w1 <- t(t(sp)/colSums(sp))
mns <- t(nic1$ls) %*% w1
v1 <- (t(nic1$ls)-mns[,1])^2 %*% w1[,1,drop=F]
cxy <- sum(apply((t(nic1$ls)-mns[,1]),2,prod)*w1[,1])

res <- util.ellipse(mns[1,1], mns[2,1], v1[1,1], cxy, v1[2,1], 1.5)


Pour le reste, je sais qu'il existe dans les packages "spatiaux" des fonctions intersecter deux polygones et calculer l'aire de cette surface, mais je ne me souviens pas comment faire. On peut aussi s'en sortir avec le package geometry :

Code : Tout sélectionner

library(geometry)
v2 <- (t(nic1$ls)-mns[,2])^2 %*% w1[,2,drop=F]
cxy2 <- sum(apply((t(nic1$ls)-mns[,2]),2,prod)*w1[,2])
res2 <- util.ellipse(mns[1,2], mns[2,2], v2[1,1], cxy2, v2[2,1], 1.5)

intersectn(do.call(cbind, res[1:2]), do.call(cbind, res2[1:2]))$ch$area


Cordialement,
Maxime

[Nicolas Romillac]

Merci, je vais tenter de me plonger dans le code!

[Nicolas Romillac]

Pour info, le code bricolé pour obtenir les coordonnées de l'ellipse de l'espèce 1

Code : Tout sélectionner

library(ade4)
#dataset environnement
a<-c(0.77,0.18,0.25,1.26)
b<-c(0.63,2.91,0.29,1.04)
env<-cbind(a,b)
dudi1 <- dudi.pca(env, scale = TRUE, scan = FALSE, nf = 3)

#dataset abondances d'espèces
sp1<-c(0,33,30,22)
sp2<-c(74,40,27,94)
sp<-cbind.data.frame(sp1,sp2)

#calcul de la niche
nic1 <- niche(dudi1, sp, scann = FALSE)

#coefficient de l'ellipse
cellipse = 1.5
csp<-1#numéro de colonne correspondant à l'espèce dont on veut les coordonnées

dfxy = nic1$ls
dfdistri = sp
label=names(sp)



w1 <- unlist(lapply(dfdistri, sum))
dfdistri <- t(t(dfdistri)/w1)

coo <- scatterutil.base(dfxy = dfxy, xax = 1, yax = 2,
        xlim = NULL, ylim = c(-3,3), grid = TRUE, addaxes = TRUE,
        cgrid = 1, include.origin = TRUE, origin = c(0, 0),
        sub = "", csub = 1, possub = "bottomleft", pixmap = NULL,
        contour = NULL, area = NULL, add.plot = FALSE)
x<-coo$x
y<-coo$y
z<-dfdistri[, csp]#le numéro de la colonne indique l'epsèce considérée
coul = rep(1, length(x))

####fonction util.ellipse#####
    util.ellipse <- function(mx, my, vx, cxy, vy, coeff) {
        lig <- 100
        epsi <- 1e-10
        x <- 0
        y <- 0
        if (vx < 0)
            vx <- 0
        if (vy < 0)
            vy <- 0
        if (vx == 0 && vy == 0)
            return(NULL)
        delta <- (vx - vy) * (vx - vy) + 4 * cxy * cxy
        delta <- sqrt(delta)
        l1 <- (vx + vy + delta)/2
        l2 <- vx + vy - l1
        if (l1 < 0)
            l1 <- 0
        if (l2 < 0)
            l2 <- 0
        l1 <- sqrt(l1)
        l2 <- sqrt(l2)
        test <- 0
        if (vx == 0) {
            a0 <- 0
            b0 <- 1
            test <- 1
        }
        if ((vy == 0) && (test == 0)) {
            a0 <- 1
            b0 <- 0
            test <- 1
        }
        if (((abs(cxy)) < epsi) && (test == 0)) {
            if (vx > vy) {
                a0 <- 1
                b0 <- 0
            }
            else {
                a0 <- 0
                b0 <- 1
            }
            test <- 1
        }
        if (test == 0) {
            a0 <- 1
            b0 <- (l1 * l1 - vx)/cxy
            norm <- sqrt(a0 * a0 + b0 * b0)
            a0 <- a0/norm
            b0 <- b0/norm
        }
        a1 <- 2 * pi/lig
        c11 <- coeff * a0 * l1
        c12 <- (-coeff) * b0 * l2
        c21 <- coeff * b0 * l1
        c22 <- coeff * a0 * l2
        angle <- 0
        for (i in 1:lig) {
            cosinus <- cos(angle)
            sinus <- sin(angle)
            x[i] <- mx + c11 * cosinus + c12 * sinus
            y[i] <- my + c21 * cosinus + c22 * sinus
            angle <- angle + a1
        }
        return(list(x = x, y = y, seg1 = c(mx + c11, my + c21,
            mx - c11, my - c21), seg2 = c(mx + c12, my + c22,
            mx - c12, my - c22)))
    }
####fin de la fonction#####

z <- z/sum(z)
m1 <- sum(x * z)
m2 <- sum(y * z)
v1 <- sum((x - m1) * (x - m1) * z)
v2 <- sum((y - m2) * (y - m2) * z)
cxy <- sum((x - m1) * (y - m2) * z)
ell <- util.ellipse(m1, m2, v1, cxy, v2, cellipse)
ell#coordonnées de l'ellipse

#plot de l'ellipse
plot(ell$y~ell$x)
polygon(ell$x, ell$y, border = coul)
segments(ell$seg1[1], ell$seg1[2], ell$seg1[3], ell$seg1[4],lty = 2, col = coul)
segments(ell$seg2[1], ell$seg2[2], ell$seg2[3], ell$seg2[4],lty = 2, col = coul)




[Logez Maxime]

Bonjour,

Si je peux me permettre tu t'es un peu embêté pour rien, je t'avais fait le travail :

Code : Tout sélectionner

cellipse = 1.5
csp<-1#numéro de colonne correspondant à l'espèce dont on veut les coordonnées

dfxy = nic1$ls
dfdistri = sp
label=names(sp)

w1 <- unlist(lapply(dfdistri, sum))
dfdistri <- t(t(dfdistri)/w1)

# ne sert à rien sauf si tu veux représenter graphiquement les valeurs
coo <- scatterutil.base(dfxy = dfxy, xax = 1, yax = 2,
        xlim = NULL, ylim = c(-3,3), grid = TRUE, addaxes = TRUE,
        cgrid = 1, include.origin = TRUE, origin = c(0, 0),
        sub = "", csub = 1, possub = "bottomleft", pixmap = NULL,
        contour = NULL, area = NULL, add.plot = FALSE)
x<-coo$x # idem
y<-coo$y # idem

z<-dfdistri[, csp]#le numéro de la colonne indique l'epsèce considérée
coul = rep(1, length(x))

z <- z/sum(z)
m1 <- sum(x * z)
m2 <- sum(y * z)
v1 <- sum((x - m1) * (x - m1) * z)
v2 <- sum((y - m2) * (y - m2) * z)
cxy <- sum((x - m1) * (y - m2) * z)
ell <- util.ellipse(m1, m2, v1, cxy, v2, cellipse)
ell


Tout la fin de ce code je te l'avais donné :

Code : Tout sélectionner


w1 <- t(t(sp)/colSums(sp)) # équivalent de z mais pour toutes les espèces
mns <- t(nic1$ls) %*% w1 # équivalent de m1 et m2 mais pour toutes les espèces
v1 <- (t(nic1$ls)-mns[,1])^2 %*% w1[,1,drop=F] # équivalent de v1 et v2
cxy <- sum(apply((t(nic1$ls)-mns[,1]),2,prod)*w1[,1]) # équivalent de cxy

# equivalent de ell
res <- util.ellipse(mns[1,1], mns[2,1], v1[1,1], cxy, v1[2,1], 1.5)


D'ailleurs il faut être vigilant sur ce que représente ces ellipses : https://pbil.univ-lyon1.fr/R/pdf/qr3.pdf
Cordialement,
Maxime

[Nicolas Romillac]

Merci!
je n'avais pas vu la fin de ton code derrière la fonction