programmation / automatisation

[E.H. [compte supprimé]]

Bonjour à tous,

Serait il possible que vous m'expliquiez comment réaliser mes propres fonctions sous R.
A savoir, que je réalise un ensemble d'opérations identiques un grand nombre de fois, et ca j'ai pas envie de m'ennerver à chaque fois avec. Cette fonction devrait me permettre, en entrant :

LeNomDeMaFonction(Table_de_travail, Reponse_analysée, k -etant le nombre de variables-, n -etant le nombre d'individus de l'échantillon-)

pouvoir réaliser les actions suivantes :

1)je réalise des aov :
...j'ai un gd nb de modeles a tester avant...
aovfc62 <- aov(Reponse~ma+an+ty+co+ta+po, data = Table)
aovfc63 <- aov(Reponse~ma+an+ty+co+ta, data = Table)
aovfc64 <- aov(Reponse~an+ty+co+ta+po, data = Table)
aovfc65 <- aov(Reponse~ty+co+ta+po+ma, data = Table)
aovfc66 <- aov(Reponse~co+ta+po+ma+an, data = Table)
aovfc67 <- aov(Reponse~ta+po+ma+an+ty, data = Table)
aovfc68 <- aov(Reponse~po+ma+an+ty+co, data = Table)

2)je récupère l'AIC avec la fonction extractAIC qui me retourne les résultats df et AIC:
...
aic62 <- extractAIC(aovfc62)
aic63 <- extractAIC(aovfc63)
aic64 <- extractAIC(aovfc64)
aic65 <- extractAIC(aovfc65)
aic66 <- extractAIC(aovfc66)
aic67 <- extractAIC(aovfc67)
aic68 <- extractAIC(aovfc68)

---
ps : quelle est la différence entre la fonction AIC() et extractAIC(), j'ai pas bien saisi, puisque les 2 me retournent des résultats différent!

par exemple aic68 me retourne :

Code : Tout sélectionner

[1]   14.0000 -704.4646

et AIC(aovfc68) me retourne :

Code : Tout sélectionner

[1] 322.0091

Les deux sont bien sensées me donner l'AIC non ?
---


3)je calcule AICc

aicc1 <- aic1 + 2*k*(k+1)/(n-k-1)
aicc2 ... etc, etc...

4)j’exporte les résultats dans un tableau de la forme :

3 colonnes : N°_modele // df // AIC // AICc

5)exporter le tout pour pouvoir l'utiliser sous excel...


--
Si l'un d'entre vous avait la gentillesse de m'aider, ce serait super !
Merci merci merci :-)

[Renaud Lancelot]

Bonjour à tous,

Serait il possible que vous m'expliquiez comment réaliser mes propres fonctions sous R.

Lire la doc, suivre une formation... Désolé d'être si brutal mais il faut mettre les mains dans le cambouis. Un bon point de départ est la doc livrée avec R, en particulier "An introduction to R". Vous trouverez bcp d'autres docs y compris sur le site CRAN. La référence classique est:

Venables, W.N., Ripley, B.D., 2000. S programming. Springer-Verlag New York, Inc., New York, 264 p.

A savoir, que je réalise un ensemble d'opérations identiques un grand nombre de fois, et ca j'ai pas envie de m'ennerver à chaque fois avec. Cette fonction devrait me permettre, en entrant :

Code : Tout sélectionner

LeNomDeMaFonction([b]Table[/b]_de_travail, [b]Reponse[/b]_analysée, [b]k[/b] -etant le nombre de variables-, [b]n[/b] -etant le nombre d'individus de l'échantillon-)

pouvoir réaliser les actions suivantes :

1)je réalise des aov :
...j'ai un gd nb de modeles a tester avant...
aovfc62 <- aov(Reponse~ma+an+ty+co+ta+po, data = Table)
aovfc63 <- aov(Reponse~ma+an+ty+co+ta, data = Table)
aovfc64 <- aov(Reponse~an+ty+co+ta+po, data = Table)
aovfc65 <- aov(Reponse~ty+co+ta+po+ma, data = Table)
aovfc66 <- aov(Reponse~co+ta+po+ma+an, data = Table)
aovfc67 <- aov(Reponse~ta+po+ma+an+ty, data = Table)
aovfc68 <- aov(Reponse~po+ma+an+ty+co, data = Table)

2)je récupère l'AIC avec la fonction extractAIC qui me retourne les résultats df et AIC:
...
aic62 <- extractAIC(aovfc62)
aic63 <- extractAIC(aovfc63)
aic64 <- extractAIC(aovfc64)
aic65 <- extractAIC(aovfc65)
aic66 <- extractAIC(aovfc66)
aic67 <- extractAIC(aovfc67)
aic68 <- extractAIC(aovfc68)



---
ps : quelle est la différence entre la fonction AIC() et extractAIC(), j'ai pas bien saisi, puisque les 2 me retournent des résultats différent!

par exemple aic68 me retourne :

Code : Tout sélectionner

[1]   14.0000 -704.4646

et AIC(aovfc68) me retourne :

Code : Tout sélectionner

[1] 322.0091

Les deux sont bien sensées me donner l'AIC non ?

Il suffit de lire l'aide de extractAIC:

For linear models with unknown scale (i.e., for lm and aov), -2log L is computed from the deviance and uses a different additive constant to logLik and hence AIC. If RSS denotes the (weighted) residual sum of squares then extractAIC uses for - 2log L the formulae RSS/s - n (corresponding to Mallows' Cp) in the case of known scale s and n log (RSS/n) for unknown scale. AIC only handles unknown scale and uses the formula n log (RSS/n) - n + n log 2π - sum log w where w are the weights.

Ces différences sont sans importance dans la mesure où seules les différences entre critères comptent pour AIC and friends.

Exemple:

Code : Tout sélectionner

> data(npk, package = "MASS")
> m1 <- aov(yield ~  N * P * K, npk)
> m2 <- aov(yield ~  (N + P + K)^2, npk)
> m3 <- aov(yield ~  N + P + K, npk)
>
> (res1 <- sapply(c("m1", "m2", "m3"), function(x) extractAIC(get(x))))
          m1       m2       m3
[1,]  8.0000  7.00000  4.00000
[2,] 88.4697 88.21144 84.58295
> (res2 <- AIC(m1, m2, m3))
   df      AIC
m1  9 158.5788
m2  8 158.3205
m3  5 154.6920
> abs(diff(res1[2, ]))
      m2       m3
0.258257 3.628491
> abs(diff(res2$AIC))
[1] 0.258257 3.628491


3)je calcule AICc

aicc1 <- aic1 + 2*k*(k+1)/(n-k-1)
aicc2 ... etc, etc...

4)j’exporte les résultats dans un tableau de la forme :

3 colonnes : N°_modele // df // AIC // AICc

5)exporter le tout pour pouvoir l'utiliser sous excel...

Voici un exemple pour vous mettre le pied à l'étrier. Tout d'abord une fonction pour calculer AIC et AICc à partir de modèles aov:

Code : Tout sélectionner

AIC.aov <- function (object, ...) {
  Call <- match.call()
  nobs <- object$df.residual + object$rank
  object <- list(object, ...)
  val <- lapply(object, function(x){
    if("aov" != class(x)[1])
      stop("Fitted objects are not all of class aov.")
    LL <- logLik(x)
    ll <- c(LL)
    n <- x$df.residual + x$rank
    if (n != nobs)
      stop("The number of observations was not the same for all the fitted objects.")
    k <- attr(LL, "df")
    aic <- -2 * ll + 2 * k
    data.frame(k = k,
               AIC = aic,
               AICc = aic + 2 * k * (k + 1)/(n - k - 1))
    })
  val <- do.call("rbind", val)
  row.names(val) <- as.character(Call[-1])
  val
  }


Utilisation:

Code : Tout sélectionner

> data(npk, package = "MASS")
> m1 <- aov(yield ~  N * P * K, npk)
> m2 <- aov(yield ~  (N + P + K)^2, npk)
> m3 <- aov(yield ~  N + P + K, npk)
>
> AIC.aov(m1, m2, m3)
   k      AIC     AICc
m1 9 158.5788 171.4359
m2 8 158.3205 167.9205
m3 5 154.6920 158.0253


Le résultat étant un data.frame, il suffit d'utiliser write.table pour l'écrire dans un fichier texte... mais quel dommage d'utiliser Excel à ce stade !

Pour répondre à vos besoins, je viens de mettre à jour le package metomet disponible sur ce site GuR, en ajoutant une méthode sic (some information criteria) pour les objets de classe aov:

Code : Tout sélectionner

> m4 <- lm(yield ~  N * P * K, npk)
> m5 <- lm(yield ~  (N + P + K)^2, npk)
> m6 <- lm(yield ~  N + P + K, npk)
>
> m7 <- glm(yield ~  N * P * K, data = npk)
> m8 <- glm(yield ~  (N + P + K)^2, data = npk)
> m9 <- glm(yield ~  N + P + K, data = npk)
>
> sic(m1, m2, m3)
    n k        LL      AIC     AICc      BIC
m1 24 9 -70.28938 158.5788 171.4359 169.1812
m2 24 8 -71.16025 158.3205 167.9205 167.7449
m3 24 5 -72.34600 154.6920 158.0253 160.5823
> sic(m4, m5, m6)
    n k        LL      AIC     AICc      BIC
m4 24 9 -70.28938 158.5788 171.4359 169.1812
m5 24 8 -71.16025 158.3205 167.9205 167.7449
m6 24 5 -72.34600 154.6920 158.0253 160.5823
> sic(m7, m8, m9)
    n k        LL      AIC     AICc      BIC
m7 24 9 -70.28938 158.5788 171.4359 169.1812
m8 24 8 -71.16025 158.3205 167.9205 167.7449
m9 24 5 -72.34600 154.6920 158.0253 160.5823


Renaud

[E.H. [compte supprimé]]

Bonjour, et merci de votre réponse.

Lire la doc, suivre une formation... Désolé d'être si brutal mais il faut mettre les mains dans le cambouis. Un bon point de départ est la doc livrée avec R, en particulier "An introduction to R". Vous trouverez bcp d'autres docs y compris sur le site CRAN. La référence classique est:

Venables, W.N., Ripley, B.D., 2000. S programming. Springer-Verlag New York, Inc., New York, 264 p.

Pas de mal pour la brutalité, je suis prêt à mettre les mains dans le cambouis, mais je ne suis pas programmeur dans l'ame et pas bcp plus statisticien. J'essaie simplement (avec grande difficulté) de répondre à la demande d'une entreprise avec ma formation d'ingénieur agri/agro et dans les quelques semaines qui me sont laissées, donc pas de possibilité de formation (en temps et moyens)... Déchiffrer des cas généraux basiques pour en faire une application spécifique n'est pas donné au premier venu (moi en l'occurence).

Le résultat étant un data.frame, il suffit d'utiliser write.table pour l'écrire dans un fichier texte... mais quel dommage d'utiliser Excel à ce stade !

Désolé, mais comme je l'ai dit, je ne suis pas programmeur, et j'ai l'avantage de bien maitriser excel... je fais donc avec les moyens du bord, même si c'ets évident qu'ils ne sont pas optimaux! mais pour ce que je cherche à en faire par la suite, je n'ai pas vraiment besoin de plus.

Merci encore à vous

[E.H. [compte supprimé]]

Bonjour,

Je ne trouve pas la réponse... comment connaitre le contenu d'une fonction ?

Par exemple, pour la fonction AIC.aov, avoir le résultat :

Code : Tout sélectionner

function (object, ...) {
  Call <- match.call()
  nobs <- object$df.residual + object$rank
  object <- list(object, ...)
  val <- lapply(object, function(x){
    if("aov" != class(x)[1])
      stop("Fitted objects are not all of class aov.")
    LL <- logLik(x)
    ll <- c(LL)
    n <- x$df.residual + x$rank
    if (n != nobs)
      stop("The number of observations was not the same for all the fitted objects.")
    k <- attr(LL, "df")
    aic <- -2 * ll + 2 * k
    data.frame(k = k,
               AIC = aic,
               AICc = aic + 2 * k * (k + 1)/(n - k - 1))
    })
  val <- do.call("rbind", val)
  row.names(val) <- as.character(Call[-1])
  val
  }

J'ai vu qu'on pouvais le faire par :
>AIC.aov
ou
>body(AIC.aov)

Mais par exemple, je cherche a connaitre le contenu de extractAIC

et ca ne me retourne que :

Code : Tout sélectionner

function (fit, scale, k = 2, ...)
UseMethod("extractAIC")
<environment: namespace:stats>

ca m'aide pas beaucoup...

[Logez Maxime]

Bonjour,

Si tu regardes extractAIC tu obtiens le résultat suivant :

Code : Tout sélectionner

function (fit, scale, k = 2, ...)
UseMethod("extractAIC")
<environment: namespace:stats>


Ce résultat te dit que selon la classe de l'objet que tu rentres en paramètre de cette fonction, la fonction qui sera utilisé sera différente. De plus si tu tapes ?extractAIC tu vois dans la section Details que c'est une fonction générique pour des objets de classes : "aov", "coxph", "glm", "lm", "negbin" et "survreg". Dans ce genre de cas il te suffit de taper "nom du package dans lequel se trouve ka fonction":::"nom de la fonction"."nom de la classe de l'objet" :

Code : Tout sélectionner

stats:::extractAIC.aov
function (fit, scale = 0, k = 2, ...)
{
    n <- length(fit$residuals)
    edf <- n - fit$df.residual
    RSS <- deviance.lm(fit)
    dev <- if (scale > 0)
        RSS/scale - n
    else n * log(RSS/n)
    c(edf, dev + k * edf)
}

# un autre exemple :
stats:::extractAIC.glm
function (fit, scale = 0, k = 2, ...)
{
    n <- length(fit$residuals)
    edf <- n - fit$df.residual
    aic <- fit$aic
    c(edf, aic + (k - 2) * edf)
}

Maxime

[Eric Pagot]

pourquoi ne pas utiliser une fenêtre de script pour écrire toutes ces instructions ? En faisant des copier-coller, il suffit de changer les n°. Cela ne prend pas pas beaucoup de temps une fois que la fonction définie fonctionne.

[E.H. [compte supprimé]]

bug ...

[E.H. [compte supprimé]]

Bonjour,

Merci pour ces réponse.
J'ai essayé d'éditer la fonction AIC, le fonctionnement est encore différent apparement, je n'ai pas réussi...

pourquoi ne pas utiliser une fenêtre de script pour écrire toutes ces instructions ? En faisant des copier-coller, il suffit de changer les n°. Cela ne prend pas pas beaucoup de temps une fois que la fonction définie fonctionne.

Parce que je fais des analyses en série sur de multiples réponses et de multiples bases de données... 8 bases, 11 réponses... le faire a chaque fois, c'est long... et j'imaginais que coder serait long, mais moins que de modifier sur de nombreuses lignes de code un chiffre. Toute facon, si j'y arrive pas à coder, ben j'aurias plus le choix


Du coup j'ai essayé d'avancer un peu, mais ca marche pas :'(

VOilà ce que j'essaye de faire :

Code : Tout sélectionner

Rep <-  // il faut assigner une réponse étudiée
Table <-  // il faut assigner une table

var1 <-  // il faut assigner une variable étudiée
var2 <-  // il faut assigner une variable étudiée
var3 <-  // il faut assigner une variable étudiée
var4 <-  // il faut assigner une variable étudiée

i<- 0

/// la syntaxe m(i) ne fonctionne pas puisqu'il comprend une fonction, m$i non plus... j(ai donc fais les tests avec les valeurs i = 0, i+1 = 1

m(i) <- var1*var2*var3*var4
M(i) <- Rep ~ m(i)
aov(i) <- aov(M(i), data=Table)
AIC(oav(i))

m(i+1) <- m(i) - var1:var2:var3:var4
M(i+1) <- Rep ~ m(i+1)
aov(i+1) <- aov(M(i+1), data=Table)
AIC(oav(i+1))

si AIC(aov(i+1)) < AIC(aov(i))
alors
m(i+2) <- m(i+1) - var1:var2:var3
sinon
m(i+2) <- m(i) - var1:var2:var3
M(i+2) <- Rep ~ m(i+2)
aov(i+2) <- aov(M(i+2), data=Table)
AIC(oav(i+2))

etc, etc...

Mais ca ne marche pas, car si j'essaie d'attribuer :
var1 <- an
var2 <- ty
etc...

il me dit que an n'est pas un objet trouvé (normal c'est du texte qui correspond au titre d'un colonne)

j'ai donc essayé
var1 <- 'an'
var2 <- 'ty'
etc...

mais la, du coup, c'est vraiment du texte, et je ne peux plus rien en faire :

m0 <- var1*var2*var3*var4
impossible car var1, 2, etc... ne sont pas des nombres, mais je veux pas que ce soit des nombres, ce sont des composantes d'une formule...

du coup je suis bloqué.

Si j'essaye :

Code : Tout sélectionner

> m0 <- 'an*ty*co*ta'
> m0
[1] "an*ty*co*ta"

> M0 <- 'ASRF3~m0'
> M0
[1] "ASRF3~m0"
> aov0<-aov(M0, data=fc)
Erreur dans terms.default(formula, "Error", data = data) :
        pas de composante terms

Je suppos equ'il le prend comme tu texte et pas comme des composantes de formule...

Alors que :

Code : Tout sélectionner

> M0<-ASRF3~an*ty*co*ta
> M0
ASRF3 ~ an * ty * co * ta

et là je peux faire des régressions dessus...


Cela veut il dire que faut que j'abandonne la programmation et que je me colle tout à la main ? Dites moi franchement, qu'au moins j'arrete de me prendre la tete à me faciliter le travail si ca sert a rien ou si je suis vraiment trop loin du but...

[Logez Maxime]

Bonjour,

Pour éditer la fonction AIC :

Code : Tout sélectionner

stats:::AIC.default
function (object, ..., k = 2)
{
    if (length(list(...))) {
        object <- list(object, ...)
        val <- lapply(object, logLik)
        val <- as.data.frame(t(sapply(val, function(el) c(attr(el,
            "df"), AIC(el, k = k)))))
        names(val) <- c("df", "AIC")
        Call <- match.call()
        Call$k <- NULL
        row.names(val) <- as.character(Call[-1])
        val
    }
    else AIC(logLik(object), k = k)
}


Pour ce qui est de ta question je pense que tu devrais mettre tes formules de modèle directement comme paramètre de la fonction, ainsi que le nom du tableau qui contient tes variables de façon à ce que tu utilises directement tes modèles :

Code : Tout sélectionner

BOB <- y~x
table <- data.frame(x=rnorm(100),y=rnorm(100))
test <- function(formula,tableau){
  AIC(aov(formula,data=tableau))
}
test(BOB,table)
[1] 286.8519

C'est peut être une piste à étudier.

Maxime

[E.H. [compte supprimé]]

Merci, effectivement ca peut etre une piste.
Simplement, mon modele 0 (m0) permet de construire les modeles successifs selon le résultats des AIC des modeles suivants... donc ca ne me permet pas de retrancher les interactions au meilleur modele et je dois donc qd meme rentrer mes modeles manuellement selon le résultat

je crois que je vais abandonner cette voie.

Question subsidiaire :

Code : Tout sélectionner

> m0<-ASRF3~an*ty*co*ta
> aov0<-aov(m0, data= fcb)
> summary(aov0)
             Df Sum Sq Mean Sq F value    Pr(>F)   
an            1  4.027   4.027 19.4000 1.786e-05 ***
ty            3  2.685   0.895  4.3124 0.0057435 **
co            5  1.433   0.287  1.3813 0.2331596   
ta            1  0.054   0.054  0.2612 0.6098867   
an:ty         3  4.533   1.511  7.2796 0.0001213 ***
an:co         5  2.507   0.501  2.4159 0.0376743 * 
ty:co         6  0.605   0.101  0.4862 0.8181343   
an:ta         1  0.092   0.092  0.4444 0.5058280   
ty:ta         3  0.223   0.074  0.3588 0.7827927   
co:ta         5  1.160   0.232  1.1181 0.3522322   
an:ty:co      1  0.362   0.362  1.7430 0.1883836   
an:ty:ta      2  0.469   0.235  1.1307 0.3250035   
an:co:ta      1  0.093   0.093  0.4457 0.5052098   
ty:co:ta      2  0.005   0.003  0.0131 0.9870310   
Residuals   186 38.605   0.208

> AIC(aov0)
[1] 323.982                     

Ici, nombre de paramètre (somme des df) = 39

Par application de la fonction de Renaud AIC.aov() :

Code : Tout sélectionner

> AIC.aov(aov0)
     k     AIC     AICc
aov0 41 323.982 342.6994

ce qui me donne nombre de param = 41

De meme, la fonction AIC:

Code : Tout sélectionner

>AIC(aov0,aov1)
     df      AIC
aov0 41 323.9820
aov1  3 318.8704

renvoie aussi 41 paramètres

Bon ok, jsuis pas doué, mais je comprend pas 41<>39 ! Pourquoi en les comptant on trouve 39, et les scripts 41 ?


---
Derniere question :

Comment, après réalisation d'une anova, afficher les coefficients estimés.
Par exemple sous Splus, je peux sortir poutr le modele ASRF3 ~ an + co + an:co (meilleur modele dans ce cas) :

Estimated K Coefficients for K-level Factor:
$"(Intercept)":
(Intercept)
0.6258806

$an:
A2005 B2006
0.1247261 -0.1247261

$co:
CoE CoF CoFJ CoM CoMF CoMFJ
-0.173574 -0.074912 -0.05842164 0.08122266 0.08631494 0.13937

$"an:co":
A2005CoE B2006CoE A2005CoF B2006CoF A2005CoFJ B2006CoFJ
0.1178722 -0.1178722 -0.271041 0.271041 0.3581173 -0.3581173

A2005CoM B2006CoM A2005CoMF B2006CoMF A2005CoMFJ B2006CoMFJ
0.02606008 -0.02606008 -0.3314506 0.3314506 0.100442 -0.100442

qui me permet de connaitre le sens des variations... indispensable ! mais j'aimerais éviter de passer sous Splus rien que pour ca.


Merci pour vos réponses, toujours bienvenues!

[Logez Maxime]

Bonjour,

je ne vois rien de bien simple a ton problème il te faudrait faire des boucles récursives pour exploiter un modèle dans son ensemble mais je ne sais pas trop comment faire. Tu peux toujours regarder aussi du côté de la fonction drop1 qui peut te servir en ayant au préalable fait des modèles de type glm plutôt que aov, les modèles linéaires étant un cas particulier de glm. De plus en passant par des glm tu n'aurais plus de problème concernant l'estimation du critère d'AIC.

Code : Tout sélectionner

data(swiss)
glm1 <- glm(Fertility~.^5,data=swiss)
drop1(glm1)
Single term deletions

Model:
Fertility ~ (Agriculture + Examination + Education + Catholic +
    Infant.Mortality)^5
                                                            Df Deviance    AIC
<none>                                                           770.02 330.80
Agriculture:Examination:Education:Catholic:Infant.Mortality  1   770.84 328.86

# tu as ainsi le AIC pour un modèle avec interaction pour toutes les variables et le AIC pour un modèle avec les interactions sauf celles de l'interaction entre toutes les variables.

glm2 <- glm(Fertility~.^4,data=swiss)
drop1(glm2)
Single term deletions

Model:
Fertility ~ (Agriculture + Examination + Education + Catholic +
    Infant.Mortality)^4
                                                   Df Deviance    AIC
<none>                                                  770.84 328.86
Agriculture:Examination:Education:Catholic          1   868.99 332.49
Agriculture:Examination:Education:Infant.Mortality  1   823.97 329.99
Agriculture:Examination:Catholic:Infant.Mortality   1   791.37 328.09
Agriculture:Education:Catholic:Infant.Mortality     1   792.41 328.15
Examination:Education:Catholic:Infant.Mortality     1   787.26 327.85

# la tu as les AIC pour les différents modèles moins une des interactions.


Tu peux aussi jouer sur le paramètre scope pour savoir quel terme tu enlèves. Avec une boucle récusrvie tu devrais pouvoir enlever au fur et a mesure des termes suivant si le critère de AIC est inférieure ou non.

En espérant t'aider un peu.

Maxime

[E.H. [compte supprimé]]

Je viens de tomber sur ce topic par hasard en faisant une recherche de cours sur l'aic... et d'exemples (que je n'ai pas trouvé, ou du moins qui ne me conviennent pas, si 'ets trop théorique, je suis perdu)

viewtopic.php?p=503

ils parlent ici de la fonction step()

quel bonheur! tout mes problemes seraient réglés avec ca...
c'est ettonant qu'en vous parlant de sélection de modele vous n'ayez pas pensé à m'orienter la dessus...

La fonction step n'est pas la panacée. Ne permet pas d'utiliser l'AICc comme critère, par exemple. De plus, en cas de modèle complexe, avec des interactions par exemple, rien n'assure que le modèle avec l'AIC le plus faible a bien été sélectionné: la fct peut s'arrêter à des minimum "locaux". Il faut souvent vérifier à la main.

par contre, le top du top aurait été qu'une fonction similaire existe non pas en employant la méthode extractAIC(), mais AIC(), voire même la fonction proposée par Renaud AIC.aov(), qui me sortirait au passage l'AICc, n, k, n/k.

existe t elle cette fonction ? sinon y' apas une bonne volonté pour m'aider à la modifier dans le sens souhaité ?

Regarder et modifier le code dans le package metomet selon vos besoins (voir sic.R). Le fichier source est disponible à ftp://ftp.cirad.fr/pub/group-r/groupe-r/Packages/Sources/metomet_0.6-5.tar.gz

et concrèetement, je n'ai pas saisi la différence entre la méthode de calcul de extractAIC et AIC. certes c'est la différence entre deux AIC qui est importante, mais je vois pas pourquoi dans un cas j'ai AIC négatif et l'autre positif... l'extrait souligné par renaud dans son premier post ne m'interpele pas du tout!

Une histoire de constante utilisée dans le calcul qui n'est pas la même dans un cas ou dans l'autre. Tout ça est expliqué dans les références citées dans l'aide des fct, ainsi que dans la fiche sur AIC disponible sur le forum.

Renaud

[Renaud Lancelot]

Désolé, j'ai édité le mail en croyant y répondre. Voir donc le mail précédent pour mes réponses.

Renaud

[E.H. [compte supprimé]]

Bonjour,
pas grave pour le message.

1)
par contre, je n'arrive pas a installer le package, c'ets la premiere fois que j'ai a faire cela. je me suis basé sur cet article : http://pbil.univ-lyon1.fr/R/fichestd/tdr11.pdf

je l'ai désarchivé et mis en .zip pour l'importer depuis le menu packages de R.

Voici le résultat :

Code : Tout sélectionner

> utils:::menuInstallLocal()
updating HTML package descriptions
> library()
//ouverture de la fenetre des librairies :
metomet                 ** No title available (pre-2.0.0 install?) **
> library(metomet)
Erreur dans library(metomet) : 'metomet' n'est pas un package valide -- a-t-il été installé < 2.0.0 ?

merci de bien vouloir m'aider (désolé d'etre si chiant). je suppose qu'il y a une autre manipulation a faire

La fonction step n'est pas la panacée. Ne permet pas d'utiliser l'AICc comme critère, par exemple.

ok pour l'AICc, j'avais remarqué... d'autant que c'ets cet indicateur qui m'interresse.

2)

De plus, en cas de modèle complexe, avec des interactions par exemple, rien n'assure que le modèle avec l'AIC le plus faible a bien été sélectionné: la fct peut s'arrêter à des minimum "locaux".

Je ne comprend pas bien...
J'ai testé la fonction, et elle prend en compte les interactions pourtant ?!
cf resultat :

Code : Tout sélectionner

> fcbm0
ASRF3 ~ an + ta + co + ty + an:ta + an:co + an:ty + ta:co + ta:ty +
    co:ty + an:ta:co + an:ta:ty + an:co:ty + ta:co:ty + an:ta:co:ty

> fcba0
Call:
   aov(formula = fcbm0, data = fcb)

Terms:
                      an       ta       co       ty    an:ta    an:co    an:ty
Sum of Squares   4.02651  0.27382  2.07052  1.82845  0.04614  5.81686  1.22930
Deg. of Freedom        1        1        5        3        1        5        3
                   ta:co    ta:ty    co:ty an:ta:co an:ta:ty ta:co:ty Residuals
Sum of Squares   0.82643  0.01560  1.18703  0.58536  0.33828  0.00542  38.60473
Deg. of Freedom        5        3        6        3        1        2       186

Residual standard error: 0.4555791
56 out of 96 effects not estimable
Estimated effects may be unbalanced

> step(fcba0)
Start:  AIC= -319.38
 ASRF3 ~ an + ta + co + ty + an:ta + an:co + an:ty + ta:co + ta:ty + 
    co:ty + an:ta:co + an:ta:ty + an:co:ty + ta:co:ty + an:ta:co:ty


Step:  AIC= -319.38
 ASRF3 ~ an + ta + co + ty + an:ta + an:co + an:ty + ta:co + ta:ty + 
    co:ty + an:ta:co + an:ta:ty + an:co:ty + ta:co:ty


Step:  AIC= -319.38
 ASRF3 ~ an + ta + co + ty + an:ta + an:co + an:ty + ta:co + ta:ty + 
    co:ty + an:ta:co + an:ta:ty + ta:co:ty

           Df Sum of Sq     RSS     AIC
- ta:co:ty  2      0.01   38.61 -323.35
- an:ta:co  2      0.21   38.82 -322.13
- an:ta:ty  1      0.26   38.86 -319.87
<none>                    38.60 -319.38

Step:  AIC= -323.35
 ASRF3 ~ an + ta + co + ty + an:ta + an:co + an:ty + ta:co + ta:ty + 
    co:ty + an:ta:co + an:ta:ty

           Df Sum of Sq     RSS     AIC
- an:ta:co  2      0.25   38.86 -325.86
- co:ty     5      1.48   40.09 -324.85
- an:ta:ty  1      0.34   38.95 -323.38
<none>                    38.61 -323.35

Step:  AIC= -325.86
 ASRF3 ~ an + ta + co + ty + an:ta + an:co + an:ty + ta:co + ta:ty + 
    co:ty + an:ta:ty

           Df Sum of Sq     RSS     AIC
- ta:co     4      1.04   39.90 -327.92
- co:ty     5      1.47   40.33 -327.48
- an:ta:ty  2      0.67   39.53 -326.00
<none>                    38.86 -325.86
- an:co     4      2.11   40.98 -321.90

Step:  AIC= -327.92
 ASRF3 ~ an + ta + co + ty + an:ta + an:co + an:ty + ta:ty + co:ty + 
    an:ta:ty

           Df Sum of Sq     RSS     AIC
- co:ty     6      1.08   40.98 -333.90
- an:ta:ty  3      0.79   40.69 -329.47
<none>                    39.90 -327.92
- an:co     4      2.34   42.24 -323.02

Step:  AIC= -333.9
 ASRF3 ~ an + ta + co + ty + an:ta + an:co + an:ty + ta:ty + an:ta:ty

           Df Sum of Sq     RSS     AIC
- an:ta:ty  3      0.54   41.51 -336.97
<none>                    40.98 -333.90
- an:co     5      2.23   43.21 -331.92

Step:  AIC= -336.97
 ASRF3 ~ an + ta + co + ty + an:ta + an:co + an:ty + ta:ty

        Df Sum of Sq     RSS     AIC
- ta:ty  3      0.05   41.56 -342.69
- an:ty  3      0.73   42.24 -339.02
- an:ta  1      0.06   41.57 -338.64
<none>                 41.51 -336.97
- an:co  5      2.43   43.94 -334.12

Step:  AIC= -342.69
 ASRF3 ~ an + ta + co + ty + an:ta + an:co + an:ty

        Df Sum of Sq     RSS     AIC
- an:ta  1      0.05   41.62 -344.41
<none>                 41.56 -342.69
- an:ty  3      1.23   42.79 -342.10
- an:co  5      2.53   44.09 -339.36

Step:  AIC= -344.41
 ASRF3 ~ an + ta + co + ty + an:co + an:ty

        Df Sum of Sq     RSS     AIC
- ta     1      0.22   41.83 -345.24
<none>                 41.62 -344.41
- an:ty  3      1.21   42.83 -343.93
- an:co  5      2.51   44.12 -341.19

Step:  AIC= -345.24
 ASRF3 ~ an + co + ty + an:co + an:ty

        Df Sum of Sq     RSS     AIC
<none>                 41.83 -345.24
- an:ty  3      1.16   42.99 -345.08
- an:co  5      2.42   44.25 -342.53

Call:
   aov(formula = ASRF3 ~ an + co + ty + an:co + an:ty, data = fcb)

Terms:
                      an       co       ty    an:co    an:ty Residuals
Sum of Squares   4.02651  2.23061  1.88796  5.72219  1.15603  41.83116
Deg. of Freedom        1        5        3        5        3       208

Residual standard error: 0.4484544
Estimated effects may be unbalanced


Je constate que c'est bien le meilleur modèle qui est sélectionné a chaque fois, et qu'il teste bien l'effet des interactions. je n'ai donc pas compris ta remarque...
encore une fois, désolé d'etre si chiant...

je commence a désespérer, mais ca me confirme que je me lancerais pas dans une thèse a la fin de l'année :-)

Merci d'avance

[Renaud Lancelot]

Bonjour,
pas grave pour le message.

1)
par contre, je n'arrive pas a installer le package, c'ets la premiere fois que j'ai a faire cela. je me suis basé sur cet article : http://pbil.univ-lyon1.fr/R/fichestd/tdr11.pdf

je l'ai désarchivé et mis en .zip pour l'importer depuis le menu packages de R.

Voici le résultat :

Code : Tout sélectionner

> utils:::menuInstallLocal()
updating HTML package descriptions
> library()
//ouverture de la fenetre des librairies :
metomet                 ** No title available (pre-2.0.0 install?) **
> library(metomet)
Erreur dans library(metomet) : 'metomet' n'est pas un package valide -- a-t-il été installé < 2.0.0 ?

merci de bien vouloir m'aider (désolé d'etre si chiant). je suppose qu'il y a une autre manipulation a faire

Je vous ai indiqué les sources du package afin que vous puissiez regarder et éditer plus facilement le code des fonctions qui s'y trouvent et voir comment ces fonctions sont organisées. Si vous voulez installer le package, il faut passer par la commande

R CMD INSTALL blabla

à exécuter depuis une fenêtre DOS (si vous êtes sous Windows). Voir le manuel "R installation and administration". Cela nécessite d'avoir installé Perl. Si vous ne voyez pas du tout de quoi il s'agit, ne vous lancez pas là-dedans tout de suite.

Si vous souhaitez simplement installer le package à partir de sa version binaire, télécharger la version .zip (dispo sur GuR) et utiliser l'option "Installer le(s) package(s) depuis des fichiers zip" du menu "Packages".

La fonction step n'est pas la panacée. Ne permet pas d'utiliser l'AICc comme critère, par exemple.

ok pour l'AICc, j'avais remarqué... d'autant que c'ets cet indicateur qui m'interresse.

2)

De plus, en cas de modèle complexe, avec des interactions par exemple, rien n'assure que le modèle avec l'AIC le plus faible a bien été sélectionné: la fct peut s'arrêter à des minimum "locaux".

Je ne comprend pas bien...

La fonction step (et stepAIC dans le package MASS) essaient de retirer, ou d'ajouter, les termes du modèle un par un. Elle s'arrête quand elle arrive à un minimum. Il s'avère parfois, mais pas toujours, qu'on parvient à diminuer l'AIC en enlevant plus d'un terme à la fois, par exemple une interaction et un des effets principaux impliqués. Ce ne serait pas surprenant avec le modèle que vous utilisez.

J'ai testé la fonction, et elle prend en compte les interactions pourtant ?!
cf resultat :

Code : Tout sélectionner

> fcbm0
ASRF3 ~ an + ta + co + ty + an:ta + an:co + an:ty + ta:co + ta:ty +
    co:ty + an:ta:co + an:ta:ty + an:co:ty + ta:co:ty + an:ta:co:ty

> fcba0
Call:
   aov(formula = fcbm0, data = fcb)

Terms:
                      an       ta       co       ty    an:ta    an:co    an:ty
Sum of Squares   4.02651  0.27382  2.07052  1.82845  0.04614  5.81686  1.22930
Deg. of Freedom        1        1        5        3        1        5        3
                   ta:co    ta:ty    co:ty an:ta:co an:ta:ty ta:co:ty Residuals
Sum of Squares   0.82643  0.01560  1.18703  0.58536  0.33828  0.00542  38.60473
Deg. of Freedom        5        3        6        3        1        2       186

Residual standard error: 0.4555791
56 out of 96 effects not estimable
Estimated effects may be unbalanced

> step(fcba0)
Start:  AIC= -319.38
 ASRF3 ~ an + ta + co + ty + an:ta + an:co + an:ty + ta:co + ta:ty + 
    co:ty + an:ta:co + an:ta:ty + an:co:ty + ta:co:ty + an:ta:co:ty


Step:  AIC= -319.38
 ASRF3 ~ an + ta + co + ty + an:ta + an:co + an:ty + ta:co + ta:ty + 
    co:ty + an:ta:co + an:ta:ty + an:co:ty + ta:co:ty


Step:  AIC= -319.38
 ASRF3 ~ an + ta + co + ty + an:ta + an:co + an:ty + ta:co + ta:ty + 
    co:ty + an:ta:co + an:ta:ty + ta:co:ty

           Df Sum of Sq     RSS     AIC
- ta:co:ty  2      0.01   38.61 -323.35
- an:ta:co  2      0.21   38.82 -322.13
- an:ta:ty  1      0.26   38.86 -319.87
<none>                    38.60 -319.38

Step:  AIC= -323.35
 ASRF3 ~ an + ta + co + ty + an:ta + an:co + an:ty + ta:co + ta:ty + 
    co:ty + an:ta:co + an:ta:ty

           Df Sum of Sq     RSS     AIC
- an:ta:co  2      0.25   38.86 -325.86
- co:ty     5      1.48   40.09 -324.85
- an:ta:ty  1      0.34   38.95 -323.38
<none>                    38.61 -323.35

Step:  AIC= -325.86
 ASRF3 ~ an + ta + co + ty + an:ta + an:co + an:ty + ta:co + ta:ty + 
    co:ty + an:ta:ty

           Df Sum of Sq     RSS     AIC
- ta:co     4      1.04   39.90 -327.92
- co:ty     5      1.47   40.33 -327.48
- an:ta:ty  2      0.67   39.53 -326.00
<none>                    38.86 -325.86
- an:co     4      2.11   40.98 -321.90

Step:  AIC= -327.92
 ASRF3 ~ an + ta + co + ty + an:ta + an:co + an:ty + ta:ty + co:ty + 
    an:ta:ty

           Df Sum of Sq     RSS     AIC
- co:ty     6      1.08   40.98 -333.90
- an:ta:ty  3      0.79   40.69 -329.47
<none>                    39.90 -327.92
- an:co     4      2.34   42.24 -323.02

Step:  AIC= -333.9
 ASRF3 ~ an + ta + co + ty + an:ta + an:co + an:ty + ta:ty + an:ta:ty

           Df Sum of Sq     RSS     AIC
- an:ta:ty  3      0.54   41.51 -336.97
<none>                    40.98 -333.90
- an:co     5      2.23   43.21 -331.92

Step:  AIC= -336.97
 ASRF3 ~ an + ta + co + ty + an:ta + an:co + an:ty + ta:ty

        Df Sum of Sq     RSS     AIC
- ta:ty  3      0.05   41.56 -342.69
- an:ty  3      0.73   42.24 -339.02
- an:ta  1      0.06   41.57 -338.64
<none>                 41.51 -336.97
- an:co  5      2.43   43.94 -334.12

Step:  AIC= -342.69
 ASRF3 ~ an + ta + co + ty + an:ta + an:co + an:ty

        Df Sum of Sq     RSS     AIC
- an:ta  1      0.05   41.62 -344.41
<none>                 41.56 -342.69
- an:ty  3      1.23   42.79 -342.10
- an:co  5      2.53   44.09 -339.36

Step:  AIC= -344.41
 ASRF3 ~ an + ta + co + ty + an:co + an:ty

        Df Sum of Sq     RSS     AIC
- ta     1      0.22   41.83 -345.24
<none>                 41.62 -344.41
- an:ty  3      1.21   42.83 -343.93
- an:co  5      2.51   44.12 -341.19

Step:  AIC= -345.24
 ASRF3 ~ an + co + ty + an:co + an:ty

        Df Sum of Sq     RSS     AIC
<none>                 41.83 -345.24
- an:ty  3      1.16   42.99 -345.08
- an:co  5      2.42   44.25 -342.53

Call:
   aov(formula = ASRF3 ~ an + co + ty + an:co + an:ty, data = fcb)

Terms:
                      an       co       ty    an:co    an:ty Residuals
Sum of Squares   4.02651  2.23061  1.88796  5.72219  1.15603  41.83116
Deg. of Freedom        1        5        3        5        3       208

Residual standard error: 0.4484544
Estimated effects may be unbalanced


Je constate que c'est bien le meilleur modèle qui est sélectionné a chaque fois, et qu'il teste bien l'effet des interactions. je n'ai donc pas compris ta remarque...
encore une fois, désolé d'etre si chiant...

je commence a désespérer, mais ca me confirme que je me lancerais pas dans une thèse a la fin de l'année :-)

Merci d'avance

Faut pas désespérer... Ces questions ne sont pas simples. Vous auriez intérêt à lire un bon bouquin sur la sélection des modèles, ce qui est une question très importante - et difficile - en stats. Voir par exemple:

Burnham, K.P., Anderson, D.R., 2002. Model selection and multimodel inference: a practical information-theoretic approach., 2nd Edition. Springer-Verlag, New-York, 496 p.

Vous y verrez que l'application aveugle de méthodes pas-à-pas n'est pas recommandée, car on peut facilement tomber sur un artefact. Il est préférable d'établir une liste pas trop grande de modèles plausibles et de les comparer tous à l'aide d'un critère d'information. Même les experts ne sont pas d'accord entre eux sur les modèles pouvant être comparés. Pour avoir le point de vue d'un puriste (B Ripley, auteur des fct step et stepAIC), voir http://web.maths.unsw.edu.au/~inge/statlearn/ripley1.pdf. Le package BMA (Bayesian model averaging) met d'ailleurs en application pas mal d'idées de ce doc avec une approche plus "Ripleyenne" que celle expliquée dans B&A. Le package BMA a été présenté dans un numéro assez récent de R-News: voir http://cran.r-project.org/doc/Rnews/Rnews_2005-2.pdf