Help! utilisation boxcox, probleme de script

Questions sur les fonctions statistiques de R

Modérateur : Groupe des modérateurs

Orianne Liet
Messages : 9
Enregistré le : 31 Aoû 2006, 03:14

Help! utilisation boxcox, probleme de script

Messagepar Orianne Liet » 08 Sep 2006, 04:22

bonjour,
je realise une analyse necessitant un changement de variable ( en effet, les donnees telles quelles ne me permettent pas de valider le test de levesne). je cherche donc a utiliser la fonction boxplot. helas, sans resiultat
le message d erreur est : "response variable must be positive"

merci davance pour votre aide

Renaud Lancelot
Messages : 2484
Enregistré le : 16 Déc 2004, 08:01
Contact :

Messagepar Renaud Lancelot » 08 Sep 2006, 12:36

Je pense que vous voulez dire boxcox et non boxplot ?
La transformation de Box et Cox consiste à rechercher dans un intervalle de valeurs de lambda, la fonction des données observées de la forme:
f(y) = (y^lambda - 1) / lambda si lambda != 0 et log(y) si lambda == 0,
maximisant la vraisemblance de cette fonction.
Il y a donc un pb si lambda prend la valeur 0 alors que y prend des valeurs négatives ou nulles.
Or, dans la fonction boxcox disponible dans le package MASS, l'argument lambda prend par défaut les valeurs: lambda = seq(-2, 2, 1/10) (voir l'aide de ?boxcox).

Code : Tout sélectionner

> seq(-2, 2, 1/10)
 [1] -2.0 -1.9 -1.8 -1.7 -1.6 -1.5 -1.4 -1.3 -1.2
[10] -1.1 -1.0 -0.9 -0.8 -0.7 -0.6 -0.5 -0.4 -0.3
[19] -0.2 -0.1  0.0  0.1  0.2  0.3  0.4  0.5  0.6
[28]  0.7  0.8  0.9  1.0  1.1  1.2  1.3  1.4  1.5
[37]  1.6  1.7  1.8  1.9  2.0


Donc 0 fait partie des valeurs balayées, et vous rencontrez alors des pbs si y <= 0.

Je suggère que vous modifiez les valeurs prises par lambda de manière à éviter le 0, par exemple:

Code : Tout sélectionner

> seq(-2, 2, length = 50)
 [1] -2.00000000 -1.91836735 -1.83673469
 [4] -1.75510204 -1.67346939 -1.59183673
 [7] -1.51020408 -1.42857143 -1.34693878
[10] -1.26530612 -1.18367347 -1.10204082
[13] -1.02040816 -0.93877551 -0.85714286
[16] -0.77551020 -0.69387755 -0.61224490
[19] -0.53061224 -0.44897959 -0.36734694
[22] -0.28571429 -0.20408163 -0.12244898
[25] -0.04081633  0.04081633  0.12244898
[28]  0.20408163  0.28571429  0.36734694
[31]  0.44897959  0.53061224  0.61224490
[34]  0.69387755  0.77551020  0.85714286
[37]  0.93877551  1.02040816  1.10204082
[40]  1.18367347  1.26530612  1.34693878
[43]  1.42857143  1.51020408  1.59183673
[46]  1.67346939  1.75510204  1.83673469
[49]  1.91836735  2.00000000


Renaud

Orianne Liet
Messages : 9
Enregistré le : 31 Aoû 2006, 03:14

Messagepar Orianne Liet » 11 Sep 2006, 07:28

Merci pour votre reponse ( oui cest bien boxocx et non boxplot, excusez moi!!)
effectivement ma variable prend a plusieurs reprises la valeur 0, ce qui exclut la possibilite d avoir lambda=0.
j'ai donc utilise la meme plage te le meme pas que vous m aviez conseille pour lamda:
script = boxcox(B~P*T,lambda=seq(-2,2,length=50),data=baby).

(B~P*T etant le modele utilise pour executer mon ANOVA. B=nombre de trous de crabes, P= distance a la mer (en pourcentage), T=position sur la cote; baby etant le tableau importe sous R).

j ai aussi essaye d avoir lamda>0, avec:
script = boxcox(B~P*T,lambda=seq(1,3,length=50),data=baby).
sans aucun resultat.


j obtiens toujours le meme message d erreur:
"response variable must be positive"

ce qui, apres votre explication aurait du etre corrige.

auriez vous une idee sur l origine de ce probleme??

Renaud Lancelot
Messages : 2484
Enregistré le : 16 Déc 2004, 08:01
Contact :

Messagepar Renaud Lancelot » 11 Sep 2006, 11:15

Difficile à dire comme ça. Il me faudrait le fichier ou un extrait permettant de reproduire le pb que vous avez.

Renaud

Tillard
Messages : 87
Enregistré le : 17 Déc 2004, 10:32

boxcox et réponse négative

Messagepar Tillard » 12 Sep 2006, 07:23

Bonjour
la reponse est dans la fonction

Code : Tout sélectionner

> library(MASS)
> methods(boxcox)
[1] boxcox.default* boxcox.formula* boxcox.lm*     

   Non-visible functions are asterisked

> MASS:::boxcox.default

function (object, lambda = seq(-2, 2, 1/10), plotit = TRUE, interp = (plotit &&
    (m < 100)), eps = 1/50, xlab = expression(lambda), ylab = "log-Likelihood",
    ...)
{
    if (is.null(object$y) || is.null(object$qr))
        stop(paste(deparse(substitute(object)), "does not have both 'qr' and 'y' components"))
    y <- object$y
    n <- length(y)
    if (any(y <= 0))
        stop("response variable must be positive")
....



des qu'une valeur de y est <=0, la fonction boxcox affiche le message.
Vous ne pouvez donc pas utiliser boxcox sur des données négatives ou nulles.

Cordialement
Emmanuel Tillard
UMR ERRC (Elevage des Ruminants en Regions Chaudes)
CIRAD - St PIERRE (La Réunion)
tel: 02 62 49 92 54

Orianne Liet
Messages : 9
Enregistré le : 31 Aoû 2006, 03:14

Messagepar Orianne Liet » 12 Sep 2006, 07:33

premierement, voici le jeu de donnees
S T P B
A 0000m 00% 0
B 0000m 00% 0
C 0000m 00% 0
A 0500m 00% 0
B 0500m 00% 3
C 0500m 00% 0
A 1000m 00% 0
B 1000m 00% 0
C 1000m 00% 0
A 1500m 00% 0
B 1500m 00% 1
C 1500m 00% 0
A 2000m 00% 3
B 2000m 00% 2
C 2000m 00% 3
A 2500m 00% 0
B 2500m 00% 0
C 2500m 00% 0
A 3000m 00% 0
B 3000m 00% 0
C 3000m 00% 0
A 3500m 00% 0
B 3500m 00% 4
C 3500m 00% 6
A 4000m 00% 0
B 4000m 00% 0
C 4000m 00% 0
A 4500m 00% 0
B 4500m 00% 0
C 4500m 00% 0
A 0000m 025% 4
B 0000m 025% 2
C 0000m 025% 5
A 0500m 025% 1
B 0500m 025% 1
C 0500m 025% 3
A 1000m 025% 1
B 1000m 025% 1
C 1000m 025% 0
A 1500m 025% 4
B 1500m 025% 0
C 1500m 025% 1
A 2000m 025% 1
B 2000m 025% 1
C 2000m 025% 2
A 2500m 025% 0
B 2500m 025% 0
C 2500m 025% 0
A 3000m 025% 0
B 3000m 025% 0
C 3000m 025% 0
A 3500m 025% 0
B 3500m 025% 0
C 3500m 025% 0
A 4000m 025% 0
B 4000m 025% 1
C 4000m 025% 1
A 4500m 025% 0
B 4500m 025% 0
C 4500m 025% 0
A 0000m 050% 4
B 0000m 050% 2
C 0000m 050% 0
A 0500m 050% 2
B 0500m 050% 0
C 0500m 050% 1
A 1000m 050% 0
B 1000m 050% 0
C 1000m 050% 0
A 1500m 050% 9
B 1500m 050% 0
C 1500m 050% 0
A 2000m 050% 0
B 2000m 050% 0
C 2000m 050% 1
A 2500m 050% 0
B 2500m 050% 0
C 2500m 050% 0
A 3000m 050% 3
B 3000m 050% 4
C 3000m 050% 2
A 3500m 050% 2
B 3500m 050% 0
C 3500m 050% 2
A 4000m 050% 1
B 4000m 050% 1
C 4000m 050% 2
A 4500m 050% 0
B 4500m 050% 0
C 4500m 050% 0
A 0000m 075% 4
B 0000m 075% 2
C 0000m 075% 3
A 0500m 075% 1
B 0500m 075% 2
C 0500m 075% 3
A 1000m 075% 0
B 1000m 075% 0
C 1000m 075% 0
A 1500m 075% 1
B 1500m 075% 0
C 1500m 075% 0
A 2000m 075% 1
B 2000m 075% 0
C 2000m 075% 0
A 2500m 075% 2
B 2500m 075% 0
C 2500m 075% 1
A 3000m 075% 3
B 3000m 075% 1
C 3000m 075% 2
A 3500m 075% 0
B 3500m 075% 1
C 3500m 075% 1
A 4000m 075% 0
B 4000m 075% 0
C 4000m 075% 1
A 4500m 075% 0
B 4500m 075% 0
C 4500m 075% 0
A 0000m 100% 0
B 0000m 100% 0
C 0000m 100% 0
A 0500m 100% 0
B 0500m 100% 3
C 0500m 100% 0
A 1000m 100% 0
B 1000m 100% 0
C 1000m 100% 1
A 1500m 100% 2
B 1500m 100% 0
C 1500m 100% 0
A 2000m 100% 0
B 2000m 100% 0
C 2000m 100% 0
A 2500m 100% 0
B 2500m 100% 0
C 2500m 100% 0
A 3000m 100% 0
B 3000m 100% 3
C 3000m 100% 1
A 3500m 100% 0
B 3500m 100% 2
C 3500m 100% 2
A 4000m 100% 1
B 4000m 100% 0
C 4000m 100% 0
A 4500m 100% 0
B 4500m 100% 0
C 4500m 100% 0

Orianne Liet
Messages : 9
Enregistré le : 31 Aoû 2006, 03:14

Messagepar Orianne Liet » 12 Sep 2006, 07:39

merci pour vos reponses
cela me permet d avancer e tsur mon projet, et sur l utilisation de R!

je viens jsute de poster le jeu de donnees.
sinon, selon ce que m a dit renaud, je ne comprends pourquoi boxcox ne peut etre utilise avec y<=0, car la fonction utilisee passe pour tout y.
mais au vu du script donnee par "tillard", je comprends le message d erreur.

alors:
1- est il possible d'utiliser une sorte d "offset" pour mes valeurs ?
2- si il m est impossible d utiliser boxcox, quelle autre methode de transformation est possible? faut il que je cherche une fonction par moi mem, au vu du graphe des residus?


merci de votre aide, en esperant ne pas abuser de votre temps et de votre patience

Renaud Lancelot
Messages : 2484
Enregistré le : 16 Déc 2004, 08:01
Contact :

Messagepar Renaud Lancelot » 12 Sep 2006, 08:21

Hey, nice shot ! En vérifiant dans un bouquin sur la régression:

Draper, N.R. and H. Smith, Applied Regression Analysis. Wiley Series in Probability and Statistics, ed. V. Barnett, et al. 1998, New York (USA): John Wiley & Sons, Inc. 706 p.

les auteurs expliquent que le calcul fait appel à la moyenne géométrique ce qui explique que les valeurs de y doivent être > 0. Le code de boxcox est d'ailleurs explicite:

Code : Tout sélectionner

> MASS:::boxcox.default
function (object, lambda = seq(-2, 2, 1/10), plotit = TRUE, interp = (plotit &&
    (m < 100)), eps = 1/50, xlab = expression(lambda), ylab = "log-Likelihood",
    ...)
{
    if (is.null(object$y) || is.null(object$qr))
        stop(paste(deparse(substitute(object)), "does not have both 'qr' and 'y' components"))
    y <- object$y
    n <- length(y)
    if (any(y <= 0))
        stop("response variable must be positive")
    xqr <- object$qr
    logy <- log(y)
    ydot <- exp(mean(logy))
[...]


Donc, il faut trouver autre chose. Difficile d'aller plus loin sans connaître la nature des données. Une approche possible est d'utiliser un modèle de régression des moindres carrés pondérés, ou des moindre carrés généralisés. Voir un bouquin sur la régression linéaire. Dans le premier cas (réponses indépendantes), lm est utilisable en spécifiant un vecteur de poids (inverse du carré des résidus). Dans le second cas (réponses répétées par sujet), voir la fct gls dans le package nlme.

Renaud

Renaud Lancelot
Messages : 2484
Enregistré le : 16 Déc 2004, 08:01
Contact :

Messagepar Renaud Lancelot » 12 Sep 2006, 08:25

Orianne Liet a écrit :premierement, voici le jeu de donnees
S T P B
A 0000m 00% 0
B 0000m 00% 0
C 0000m 00% 0
A 0500m 00% 0
B 0500m 00% 3
[...]


Il faudrait une description des données: réponse, variables explicatives,...

Renaud

Orianne Liet
Messages : 9
Enregistré le : 31 Aoû 2006, 03:14

Messagepar Orianne Liet » 12 Sep 2006, 09:07

Bonjour,

en fait, j avais psote la description plus haut:
B=nombre de trous de crabes
P= distance a la mer (en pourcentage)
T=position sur la cote
S=repetition

la reponse est donc B, variables explicatives = P et T

modele utilise en anova: B~P*T


j espere que mes explications suffiront, je ne suis pas sur d avoir tout decrit

Renaud Lancelot
Messages : 2484
Enregistré le : 16 Déc 2004, 08:01
Contact :

Messagepar Renaud Lancelot » 12 Sep 2006, 09:17

OK. Donc clairement, c'est pas de la régression linéaire mais plutôt de la régression de Poisson (voire Poisson surdispersé si variance >> moyenne). A voir avec glm dans un premier temps. J'essaierai de jeter un coup d'oeil mais pas le temps tt de suite.

Renaud

Renaud Lancelot
Messages : 2484
Enregistré le : 16 Déc 2004, 08:01
Contact :

Messagepar Renaud Lancelot » 12 Sep 2006, 21:49

Aïe. Il ne semble pas y avoir bcp de structure dans les données, ou il manque des variables explicatives importantes:

Code : Tout sélectionner

## récup des données: sans objet pour vous
Crabs <- read.table("d:/analyses/travail/data/crabs.txt", header = TRUE)

## transformation et changement de nom des variables
Crabs$distance <- as.numeric(substring(Crabs$T, 1, 4))
Crabs$position <- Crabs$P
Crabs$count <- Crabs$B

# exploration graphique
library(lattice)

# données par répétition
xyplot(count ~ distance | position * S, data = Crabs,
  panel = function(x, y){
    panel.xyplot(x, y, cex = 1.1, type = "o", pch = 19)
    })

# données groupées
xyplot(count ~ distance | position, data = Crabs,
  layout = c(1, 5), aspect = "xy",
  panel = function(x, y){
    panel.xyplot(x, jitter(y))
    Y <- tapply(y, x, mean)
    X <- unique(x)
    llines(X, Y)
    })

# conclusion: pas de structure visible, ou cohérente entre
# les différentes combinasions de variables explicatives.
# Je ne sais pas ce qu'on peut faire d'utile

# exemple de modèle de Poisson
fm1 <- glm(count ~ distance * position, data = Crabs, family = poisson)
summary(fm1)
 
# package metomet sur le site GuR
library(metomet)
# calcul du Chi2 de Pearson
gof(fm)
# ==> pb de surdispersion... Mais pas la peine
# d'aller plus loin sans covariable plus intéressante.


Renaud

Orianne Liet
Messages : 9
Enregistré le : 31 Aoû 2006, 03:14

Messagepar Orianne Liet » 13 Sep 2006, 02:06

qunad vous dites "il manque des variables explicatives importantes"
cela signifie donc que l on ne peut rien tirer des resultats, et que les explications sont ailleurs?
je n ai donc rien a nalyser en qq sorte??

sinon, j ai rentre mes donnes sous la forme:
T= 0000m, 0500m, 1000m,...,4500m
P= 0%,..., 100%

mais comment etre sur de la notation?? car en rentrant
T= 0, 500, 1000, 1500....
p= 0, 0.25, 0.5

l anova ne renvoie pas la meme chose
d ailleurs, vous avez utilise cette derniere notation...?? qu en pensez?

Orianne Liet
Messages : 9
Enregistré le : 31 Aoû 2006, 03:14

Messagepar Orianne Liet » 13 Sep 2006, 02:07

oups j ai oublie qqchoses...

impossible de downloader le package metomet qqsoit le lien, le telechargement se lance, mais apres rien ne bouge...
et aussi
j ai tape votre script, et pour la loi de poisson, certains trucs m etonnent.
on a seulement une ligne pour la distance T (que vous avez passe en numerique)
et pour les positons p, on n obtient aucun resutat pour P=0%

Renaud Lancelot
Messages : 2484
Enregistré le : 16 Déc 2004, 08:01
Contact :

Messagepar Renaud Lancelot » 14 Sep 2006, 10:19

Orianne Liet a écrit :oups j ai oublie qqchoses...

impossible de downloader le package metomet qqsoit le lien, le telechargement se lance, mais apres rien ne bouge...
et aussi
j ai tape votre script, et pour la loi de poisson, certains trucs m etonnent.
on a seulement une ligne pour la distance T (que vous avez passe en numerique)
et pour les positons p, on n obtient aucun resutat pour P=0%


1) Je viens de tester de mon côté: cela fonctionne normalement. Merci de réessayer et de me dire si le problème persiste: je verrai avec l'administrateur du site ftp.

2) Concernant le script: ce qui est important, c'est surtout l'analyse exploratoire qui montre qu'il n'y a pas de relation simple et homogène, pour les différents panels, entre les comptages et la distance, conditionnellement à la position. C'est pour faciliter cette analyse que j'ai transformé la distance en variable numérique. Si, par chance, on avait pu mettre en évidence une relation simple, on aurait pu modéliser l'influence de la distance avec une relation linéaire ou quadratique (ou autre...), nécessitant le calcul de peu de paramètres: recherche d'un modèle parsimonieux. Dans l'exemple de modèle que je vous ai donné, la distance est codée sous forme numérique, et la ligne que vous voyez dans le résultat pour cette variable vous donne le coefficient de la pente de l'équation log(count) = b0 + b1 * distance + (b2,...) * position + termes d'interaction. C'était juste un exemple...

Renaud

Concernant l'absence de résultats pour P = 0%, c'est parce cette variable est un facteur, et que la modalité de référence (0%) est intégrée dans la constante du modèle. C'est normal. Si vous voulez la voir apparaître, vous pouvez spécifier un modèle sans intercept, par exemple:

Code : Tout sélectionner

glm(count ~ -1 + distance * position, data = Crabs, family = poisson)


Renaud


Retourner vers « Archives : Fonctions statistiques »

Qui est en ligne

Utilisateurs parcourant ce forum : Aucun utilisateur enregistré et 1 invité