Comment effectuer une transformation isométrique de log-ratio

J'ai des données sur les comportements de mouvement (temps passé à dormir, sédentaire et faire de l'activité physique) qui totalisent environ 24 (comme en heures par jour). Je veux créer une variable qui capture le temps relatif passé dans chacun de ces comportements - on m'a dit qu'une transformation isométrique de log-ratio accomplirait cela.

Il semble que je devrais utiliser la fonction ilr dans R, mais je ne trouve aucun exemple réel avec du code. Où est-ce que je commence?

Les variables dont je dispose sont le temps passé à dormir, le temps sédentaire moyen, l'activité physique moyenne moyenne légère, l'activité physique moyenne modérée et l'activité physique vigoureuse moyenne. Le sommeil a été auto-déclaré, tandis que les autres sont des moyennes de jours valides de données d'accéléromètre. Donc, pour ces variables, les cas ne totalisent pas exactement 24.

Ma conjecture: je travaille en SAS, mais il semble que R sera beaucoup plus facile à utiliser pour cette partie. Importez donc d'abord les données avec uniquement les variables d'intérêt. Utilisez ensuite la fonction acomp (). Ensuite, je ne peux pas comprendre la syntaxe de la fonction ilr (). Toute aide serait très appréciée.

La transformation ILR (Isometric Log-Ratio) est utilisée dans l'analyse des données de composition. Toute observation donnée est un ensemble de valeurs positives résumant à l'unité, telles que les proportions de produits chimiques dans un mélange ou les proportions du temps total passé dans diverses activités. L'invariant somme-à-unité implique que bien qu'il puisse y avoir composantes à chaque observation, il n'y a que fonctionnellement indépendantes. (Géométriquement, les observations se trouvent sur un simplexe dimensionnel dans l' espace euclidien dimensionnel \ mathbb {R} ^ k . Cette nature simplicial se manifeste dans les formes triangulaires des diagrammes de dispersion des données simulées illustrées ci-dessous.)$k\ge 2$$k-1$$k-1$$k$$\mathbb{R}^k$

En règle générale, les distributions des composants deviennent «plus agréables» lorsque le journal est transformé. Cette transformation peut être mise à l'échelle en divisant toutes les valeurs d'une observation par leur moyenne géométrique avant de prendre les journaux. (De manière équivalente, les journaux des données dans toute observation sont centrés en soustrayant leur moyenne.) Ceci est connu sous le nom de transformation "Centré Log-Ratio", ou CLR. Les valeurs résultantes se trouvent toujours dans un hyperplan dans , car la mise à l'échelle fait que la somme des journaux est nulle. L'ILR consiste à choisir n'importe quelle base orthonormée pour cet hyperplan: les coordonnées de chaque observation transformée deviennent ses nouvelles données. De manière équivalente, l'hyperplan est tourné (ou réfléchi) pour coïncider avec le plan avec la disparition de$\mathbb{R}^k$$k-1$$k^\text{th}$coordonner et on utilise les premières coordonnées . (Parce que les rotations et les réflexions préservent la distance, ce sont des isométries , d'où le nom de cette procédure.)$k-1$

Tsagris, Preston et Wood affirment qu '"un choix standard de [la matrice de rotation] est la sous-matrice Helmert obtenue en supprimant la première ligne de la matrice Helmert."$H$

La matrice Helmert d'ordre est construite de manière simple (voir Harville p. 86 par exemple). Sa première ligne est de s. La ligne suivante est l'une des plus simples pouvant être orthogonales à la première ligne, à savoir . La ligne est parmi les plus simples orthogonales à toutes les lignes précédentes: ses premières entrées sont s, ce qui garantit qu'elle est orthogonale aux lignes et son entrée est définie sur pour la rendre orthogonale à la première ligne (c'est-à-dire que ses entrées doivent être égales à zéro). Toutes les lignes sont ensuite redimensionnées en unités de longueur.$k$$1$$(1, -1, 0, \ldots, 0)$$j$$j-1$$1$$2, 3, \ldots, j-1$$j^\text{th}$$1-j$

Voici, pour illustrer le modèle, la matrice Helmert avant que ses lignes aient été redimensionnées:$4\times 4$

(Edit ajouté en août 2017) Un aspect particulièrement agréable de ces «contrastes» (qui sont lus ligne par ligne) est leur interprétabilité. La première ligne est supprimée, laissant lignes restantes pour représenter les données. La deuxième ligne est proportionnelle à la différence entre la deuxième variable et la première. La troisième ligne est proportionnelle à la différence entre la troisième variable et les deux premières. Généralement, la ligne ( ) reflète la différence entre la variable et toutes celles qui la précèdent, les variables . Cela laisse la première variable$k-1$$j$$2\le j \le k$$j$$1, 2, \ldots, j-1$$j=1$comme "base" pour tous les contrastes. J'ai trouvé ces interprétations utiles en suivant l'ILR par analyse en composantes principales (PCA): il permet d'interpréter les chargements, au moins grossièrement, en termes de comparaisons entre les variables d'origine. J'ai inséré une ligne dans l' Rimplémentation de ilrci-dessous qui donne aux variables de sortie des noms appropriés pour aider à cette interprétation. (Fin de l'édition.)

Depuis Rfournit une fonction contr.helmertpour créer de telles matrices (bien que sans mise à l'échelle, et avec des lignes et des colonnes négatives et transposées), vous n'avez même pas besoin d'écrire le code (simple) pour le faire. En utilisant cela, j'ai implémenté l'ILR (voir ci-dessous). Pour l'exercer et le tester, j'ai généré tirages indépendants à partir d'une distribution de Dirichlet (avec les paramètres 1, 2, , ) et tracé leur matrice de nuage de points. Ici, .$1000$$1,2,3,4$$k=4$

Les points s'agglutinent tous près des coins inférieurs gauches et remplissent des zones triangulaires de leurs zones de traçage, comme cela est caractéristique des données de composition.

Leur ILR n'a que trois variables, à nouveau tracées sous forme de matrice de nuage de points:

Cela semble en effet plus agréable: les diagrammes de dispersion ont acquis des formes de «nuage elliptique» plus caractéristiques, mieux adaptées à des analyses de second ordre telles que la régression linéaire et l'ACP.

Tsagris et al. généraliser le CLR en utilisant une transformation de Box-Cox, qui généralise le logarithme. (Le journal est une transformation de Box-Cox avec le paramètre ) Il est utile parce que, comme le soutiennent les auteurs (correctement à mon humble avis), dans de nombreuses applications, les données doivent déterminer leur transformation. Pour ces données de Dirichlet, un paramètre de (qui est à mi-chemin entre aucune transformation et une transformation logarithmique) fonctionne à merveille:$0$$1/2$

"Beau" se réfère à la description simple que cette image permet: au lieu d'avoir à spécifier l'emplacement, la forme, la taille et l'orientation de chaque nuage de points, nous devons seulement observer que (pour une excellente approximation) tous les nuages ​​sont circulaires avec des rayons similaires . En effet, le CLR a simplifié une description initiale nécessitant au moins 16 chiffres en un qui ne nécessite que 12 chiffres et l'ILR l'a réduit à seulement quatre chiffres (trois emplacements univariés et un rayon), au prix de la spécification du paramètre ILR de - un cinquième numéro. Lorsque de telles simplifications spectaculaires se produisent avec des données réelles, nous pensons généralement que nous sommes sur quelque chose: nous avons fait une découverte ou obtenu un aperçu.$1/2$

Cette généralisation est implémentée dans la ilrfonction ci-dessous. La commande pour produire ces variables "Z" était simplement

z <- ilr(x, 1/2)

Un avantage de la transformation de Box-Cox est son applicabilité aux observations qui incluent des vrais zéros: elle est toujours définie à condition que le paramètre soit positif.

Les références

Michail T. Tsagris, Simon Preston et Andrew TA Wood, Une transformation de puissance basée sur les données pour les données de composition . arXiv: 1106.1451v2 [stat.ME] 16 juin 2011.

David A. Harville, Algèbre matricielle du point de vue d'un statisticien . Springer Science & Business Media, 27 juin 2008.

Voici le Rcode.

#
# ILR (Isometric log-ratio) transformation.
# `x` is an `n` by `k` matrix of positive observations with k >= 2.
#
ilr <- function(x, p=0) {
  y <- log(x)
  if (p != 0) y <- (exp(p * y) - 1) / p       # Box-Cox transformation
  y <- y - rowMeans(y, na.rm=TRUE)            # Recentered values
  k <- dim(y)[2]
  H <- contr.helmert(k)                       # Dimensions k by k-1
  H <- t(H) / sqrt((2:k)*(2:k-1))             # Dimensions k-1 by k
  if(!is.null(colnames(x)))                   # (Helps with interpreting output)
    colnames(z) <- paste0(colnames(x)[-1], ".ILR")
  return(y %*% t(H))                          # Rotated/reflected values
}
#
# Specify a Dirichlet(alpha) distribution for testing.
#
alpha <- c(1,2,3,4)
#
# Simulate and plot compositional data.
#
n <- 1000
k <- length(alpha)
x <- matrix(rgamma(n*k, alpha), nrow=n, byrow=TRUE)
x <- x / rowSums(x)
colnames(x) <- paste0("X.", 1:k)
pairs(x, pch=19, col="#00000040", cex=0.6)
#
# Obtain the ILR.
#
y <- ilr(x)
colnames(y) <- paste0("Y.", 1:(k-1))
#
# Plot the ILR.
#
pairs(y, pch=19, col="#00000040", cex=0.6)

Pour votre cas d'utilisation, il est probablement correct de tout réduire à un. Le fait que les chiffres ne correspondent pas exactement à 24 ajoutera un peu de bruit supplémentaire aux données, mais cela ne devrait pas gâcher les choses.

Comme @whuber l'a correctement dit, étant donné qu'il s'agit de proportions, nous devons tenir compte des dépendances entre les variables (puisqu'elles s'additionnent à une). La transformation ilr traite cela de manière appropriée, car elle transforme les variables en pour les proportions$\mathbb{R}^{D-1}$$D$

Tous les détails techniques mis à part, il est important de savoir comment interpréter correctement les données transformées. En fin de compte, la transformation ilr fait simplement référence aux rapports logiques des groupes. Mais il le définit par rapport à une hiérarchie prédéfinie. Si vous définissez une hiérarchie comme suit

chaque variable transformée peut être calculée comme

$b_i = \sqrt{\frac{rs}{r + s}}\ln \frac{g(R_i)}{g(S_i)}$

où représente un nœud interne dans la hiérarchie, définit une partition de variables correspondant à , définit l'autre partition de variables correspondant à et fait référence à la moyenne géométrique. Ces variables transformées sont également appelées soldes.$i$$R_i$$i$$S_i$$i$$g(...)$

La question suivante est donc la suivante: comment définissez-vous votre hiérarchie de variables? Cela dépend vraiment de vous, mais si vous avez trois variables, il n'y a pas trop de combinaisons avec lesquelles jouer. Par exemple, vous pouvez simplement définir la hiérarchie à

                        /-A
            /(A|B)-----|
-(AB|C)----|            \-B
           |
            \-C

où Areprésente le temps passé à dormir, Breprésente le temps passé avec des sédentaires, Creprésente le temps consacré à l'activité physique (A|B)représente le rapport de log normalisé entre (c'est-à-dire ), et fait référence au rapport de log normalisé entre , et (c'est-à-dire ) . S'il y a beaucoup de variables, je vérifie une partie du travail effectué avec les soldes principaux$A$ $B$$\frac{1}{\sqrt{2}}\ln \frac{A}{B}$$(AB|C)$$A$$B$$C$$\frac{\sqrt{2}}{\sqrt{3}} \ln \frac{AB}{C}$

Mais pour en revenir à votre question initiale, comment pouvez-vous utiliser ces informations pour effectuer la transformation ilr?

Si vous utilisez R, je vérifie le package de compositions

Pour utiliser ce package, vous devrez comprendre comment créer une partition binaire séquentielle (SBP), c'est-à-dire comment définir la hiérarchie. Pour la hiérarchie définie ci-dessus, vous pouvez représenter le SBP avec la matrice suivante.

        A  B  C
(A|B)   1 -1  0
(AB|C)  1  1 -1

où les valeurs positives représentent les variables du numérateur, les valeurs négatives représentent les variables du dénominateur et les zéros représentent l'absence de cette variable dans le solde. Vous pouvez créer la base orthonormée à l'aide balanceBasedu SBP que vous avez défini.
Une fois que vous avez cela, vous devriez pouvoir passer dans votre tableau des proportions avec la base que vous avez calculée ci-dessus.

Je vérifierais cette référence pour la définition originale des soldes

Les messages ci-dessus répondent à la question de savoir comment construire une base ILR et obtenir vos soldes ILR. Pour ajouter à cela, le choix de quelle base peut faciliter l'interprétation de vos résultats.

Vous pouvez être intéressé par une partition la partition suivante:

(1) (sommeil, sédentaire | activité physique) (2) (sommeil | sédentaire).

Puisque vous avez trois parties dans votre composition, vous obtiendrez deux soldes ILR à analyser. En configurant la partition comme ci-dessus, vous pouvez obtenir des soldes correspondant à "actif ou non" (1) et "quelle forme d'inactivité" (2).

Si vous analysez chaque solde ILR séparément, par exemple en effectuant une régression par rapport à l'heure du jour ou à la période de l'année pour voir s'il y a des changements, vous pouvez interpréter les résultats en termes de changements dans "actif ou non" et de changements dans "quelle forme d'inactivité".

Si, en revanche, vous exécuterez des techniques comme PCA qui obtiennent une nouvelle base dans l'espace ILR, vos résultats ne dépendront pas de votre choix de partition. En effet, vos données existent dans l'espace CLR, le plan D-1 orthogonal au vecteur unique et les équilibres ILR sont différents choix d'axes normalisés pour décrire la position des données sur le plan CLR.