fonction predict avec option newdata sur betabin

[Vladimir Grosbois]

Bonjour
Je rencontre un problème pour appliquer la fonction predict avec option newdata dans la cas d'un modèle construit avec la fonction betabin du package AOD

J'ai opéré exactement comme dans les exemples donnés dans la documentation de la fonction betabin (il y a un exemple dans lequel la fonction predict est appliquée à un modèle béta-binomial et une table qui contient les combinaisons des modalités des variables explicatives pour lesquelles on veut des prédictions)

Et j'obtiens le message d'erreur suivant
Error in X %*% b : non-conformable arguments


A noter que la fonction predict sans l'option newdata fonctionne bien, mais ce n'est pas satisfaisant étant donné que cela fournit les prédiction uniquement pour les combinaisons de modalités des variables explicatives présentes dans le jeu de donné sur lequel le modèle est ajusté. Or, j'ai besoin de prédictions pour d'autres combinaisons...

A noter aussi que parmi les variables explicatives du modèle, certaines sont catégorielles et d'autres quantitatives.

Quelqu'un a-t'il déjà rencontré ce genre de problème et pourrait il m'orienter vers une solution

Merci d'avance

[Logez Maxime]

Bonjour,

as-tu fait attention que le data.frame qui contient les données que tu veux prédire a bien les mêmes noms de variable ou de modalités que celui utilisé pour la prédiction ?

Maxime

[Vladimir Grosbois]

Bonjour,

as-tu fait attention que le data.frame qui contient les données que tu veux prédire a bien les mêmes noms de variable ou de modalités que celui utilisé pour la prédiction ?

Maxime

Bonsoir
Oui j'ai bien vérifié et je viens de vérifier à nouveau, pas de souci à ce niveau là.
Mais merci pour la suggestion, ça aurait pu être ça....
Bonne soirée

[Matthieu Lesnoff]

Valdimir, difficile d'identifier le pb sans connaitre ton data.frame, le modèle que tu as ajusté et ton newdata

Ton plan est il complet ou incomplet ? Dans le dernier cas, certaines combinaisons de modalités peuvent ne pas être estimables. Mais je pense que le pb est ailleurs.

Avec un newdata, le coeur de predict pour betabin fait simplement :

Code : Tout sélectionner

f <- object@formula[-2]
X <- model.matrix(object = f, data = newdata)
eta <- as.vector(X %*% b)


(avec "object" étant le modèle ajusté, X la matrice plan, b le coefficient des variables fixes, et eta=g^(-1)(p)).

En remarque, tu peux toi même facilement construire n'importe quelle prédiction (estimable) en définissant correctement une matrice X par rapport à b ou une partie de b.

Pour ton pb, vérifie que ton X est correctement construite par rapport à b (c'est le point central). Il faut que chaque ligne de X soit multipliable par b, et que chaque colonne corresponde au bon coefficient (attention à l'ordre des varibales). Tu peux généer X à partir de model.matrix().

J'ai testé un exemple ci-dessous en incluant une variable continue :

Code : Tout sélectionner

> tmp <- orob2
> tmp$x <- rnorm(n = nrow(tmp), mean = 3, sd = 6)
> tmp
   seed     root  n  y           x
1   O75     BEAN 39 10  3.50723490
2   O75     BEAN 62 23 -7.19829544
3   O75     BEAN 81 23  0.04066238
4   O75     BEAN 51 26  4.16604617
5   O75     BEAN 39 17  8.14147741
6   O75 CUCUMBER  6  5 10.44778199
7   O75 CUCUMBER 74 53  0.55930554
8   O75 CUCUMBER 72 55  4.92978772
9   O75 CUCUMBER 51 32 -0.51976112
10  O75 CUCUMBER 79 46 14.82109790
11  O75 CUCUMBER 13 10 -0.72310558
12  O73     BEAN 16  8  6.07792514
13  O73     BEAN 30 10  1.24799338
14  O73     BEAN 28  8  5.73092162
15  O73     BEAN 45 23  1.81535844
16  O73     BEAN  4  0  2.13382292
17  O73 CUCUMBER 12  3 12.72234807
18  O73 CUCUMBER 41 22  0.94993457
19  O73 CUCUMBER 30 15  4.22959492
20  O73 CUCUMBER 51 32  9.35200937
21  O73 CUCUMBER  7  3  1.30378450

> fm <- betabin(formula = cbind(y, n - y) ~ seed * root + x, random = ~ 1, data = tmp)

> fm
Beta-binomial model
-------------------
betabin(formula = cbind(y, n - y) ~ seed * root + x, random = ~1,
    data = tmp)

Convergence was obtained after 571 iterations.

Fixed-effect coefficients:
                       Estimate Std. Error    z value Pr(> |z|)
(Intercept)          -4.166e-01  2.255e-01 -1.848e+00 6.463e-02
seedO75              -1.171e-01  2.754e-01 -4.253e-01 6.706e-01
rootCUCUMBER          5.426e-01  2.986e-01  1.817e+00 6.918e-02
x                    -8.254e-03  1.875e-02 -4.403e-01 6.597e-01
seedO75:rootCUCUMBER  8.097e-01  3.764e-01  2.151e+00 3.147e-02

Overdispersion coefficients:
                 Estimate Std. Error   z value   Pr(> z)
phi.(Intercept) 1.183e-02  1.125e-02 1.052e+00 1.465e-01

Log-likelihood statistics
   Log-lik      nbpar    df res.   Deviance        AIC       AICc
-5.367e+01          6         15  3.075e+01  1.193e+02  1.253e+02
 

> New <- expand.grid(seed = levels(orob2$seed), root = levels(orob2$root))
> New$x <-  rnorm(n = nrow(New), mean = 3, sd = 6)
> New
  seed     root         x
1  O73     BEAN -1.505671
2  O75     BEAN  2.449098
3  O73 CUCUMBER  6.790455
4  O75 CUCUMBER  7.462190

> model.matrix(formula(cbind(y, n - y) ~ seed * root + x)[-2], New)
  (Intercept) seedO75 rootCUCUMBER         x seedO75:rootCUCUMBER
1           1       0            0 -1.505671                    0
2           1       1            0  2.449098                    0
3           1       0            1  6.790455                    0
4           1       1            1  7.462190                    1

> coef(fm)
         (Intercept)              seedO75         rootCUCUMBER                    x seedO75:rootCUCUMBER
        -0.416630013         -0.117144209          0.542591485         -0.008254265          0.809666958

> predict(fm, New)
[1] 0.4003032 0.3649393 0.5174707 0.6806780


M.

PS: tu peux aussi comparer les résultats avec le modèle beta-binomial ajusté avec gamlss.

[Vladimir Grosbois]

Bonjour tout le monde,

Merci beaucoup Matthieu pour cet éclaircissement qui m'a permis d'identifier mon problème. C'était bien un problème de correspondance entre le vecteur contenant les estimations des coefficients du modèle et la matrice de design contenant les combinaisons des facteurs pour lesquelles je veux des estimations.

Ce que je retiens c'est que le nombre de colonnes de la matrice de design doit être égal au nombre de paramètres estimé par le modèle.

En l'occurrence, mon problème résultait du fait que j'utilisais la formule suivante pour construire la table contenant les combinaisons des variables explicatives pour lesquelles j'avais besoin d'une valeur prédite

New <- expand.grid(
year = levels(duck$year),
sample = levels(duck$sample),
genus = levels(duck$genus),
d_tot = seq(-1,1,0.2),
m_ann_tp= 0,
ts_beg_pws=0,
ts_1sep=0,
ts_1sep2=0)

duck étant la table de données sur laquelle le modèle a été ajusté.

Or cette table duck ne correspond pas aux données brutes, on a éliminé certains enregistrements, notamment tous les enregistrements de trois des 7 années. De ce fait, year étant un facteur, duck$year contient 7 niveau donc trois n'ont aucun enregistrement et donc aucun coefficient associé dans le modèle.

La modif suivante a résolu le problème

New <- expand.grid(
year = c("d","e","f","g"),
sample = levels(duck$sample),
genus = levels(duck$genus),
d_tot = seq(-1,1,0.2),
m_ann_tp= 0,
ts_beg_pws=0,
ts_1sep=0,
ts_1sep2=0)

[Vladimir Grosbois]

Au fait, et pour en finir avec les problèmes de prédiction
Concernant la construction "à la main" des prédictions avec intervalles de confiance, on doit utiliser, pour obtenir des erreurs standard des prédictions, la matrice de variance-covariance des estimations des coefs du modèle.
Quelqu'un saurait il comment on obtient cette matrice de variance-covariance ?
Merci

[Matthieu Lesnoff]

Quelqu'un saurait il comment on obtient cette matrice de variance-covariance ?

C'est le slot @varparam

Code : Tout sélectionner

> fm <- betabin(formula = cbind(y, n - y) ~ seed * root, random = ~ 1, data = orob2)
> round(fm@varparam, digits = 4)
                     (Intercept) seedO75 rootCUCUMBER seedO75:rootCUCUMBER phi.(Intercept)
(Intercept)               0.0477 -0.0479      -0.0469               0.0469          -2e-04
seedO75                  -0.0479  0.0749       0.0467              -0.0734           4e-04
rootCUCUMBER             -0.0469  0.0467       0.0881              -0.0881          -2e-04
seedO75:rootCUCUMBER      0.0469 -0.0734      -0.0881               0.1428           2e-04
phi.(Intercept)          -0.0002  0.0004      -0.0002               0.0002           1e-04

Tu peux voir tous les slots disponibles en faisant, par exemple, attributes(fm).

Au fait, et pour en finir avec les problèmes de prédiction
Concernant la construction "à la main" des prédictions avec intervalles de confiance, on doit utiliser, pour obtenir des erreurs standard des prédictions.

Tu peux les avoir directement avec predict(..., se = TRUE).

Les IC calculés avec des "se" supposent la normalité asympotique des estimateurs ML. Si tu veux ce type d'IC sur p, il vaut mieux d'abort calculer des IC sur l'echelle du lien (c'est adire sur X*b) à partir des "se" données par predict(, se = TRUE, type = "link"), puis calculer les IC sur p en inversant le lien. Les "se" calculées directement sur p sont approximatives (delta method, comme dans la fonction glm), et l'hypothèse de normalité est souvent moins fiable sur p que sur b.

M.

[Renaud Lancelot]

Il y a une méthode vcov pour betabin. Il est en général préférable d'utiliser ces utilitaires plutôt que d'accéder directement aux slots.

Code : Tout sélectionner

> library(aod)
Package aod, version 1.2
> fm <- betabin(formula = cbind(y, n - y) ~ seed * root, random = ~ 1, data = orob2)
> vcov(fm)
                     (Intercept)     seedO75 rootCUCUMBER seedO75:rootCUCUMBER
(Intercept)           0.04765561 -0.04787638  -0.04685347           0.04685561
seedO75              -0.04787638  0.07491924   0.04670550          -0.07342921
rootCUCUMBER         -0.04685347  0.04670550   0.08810137          -0.08809994
seedO75:rootCUCUMBER  0.04685561 -0.07342921  -0.08809994           0.14283430

[Vladimir Grosbois]

Merci Matthieu, merci Renaud
Tout est clair
A bientôt