Comment dessiner un graphique d'interaction avec des intervalles de confiance?

11

Mes tentatives:

  1. Je n'ai pas pu obtenir d'intervalles de confiance interaction.plot()

  2. et d'autre part plotmeans()du package 'gplot' n'afficherait pas deux graphiques. De plus, je ne pouvais pas imposer deux plotmeans()graphes l'un sur l'autre car par défaut les axes sont différents.

  3. J'ai eu un certain succès en utilisant le plotCI()package 'gplot' et en superposant deux graphiques mais la correspondance de l'axe n'était toujours pas parfaite.

Des conseils sur la façon de créer un graphique d'interaction avec des intervalles de confiance? Soit par une fonction, soit des conseils sur la superposition plotmeans()ou des plotCI()graphiques.

exemple de code

br=structure(list(tangle = c(140L, 50L, 40L, 140L, 90L, 70L, 110L, 
150L, 150L, 110L, 110L, 50L, 90L, 140L, 110L, 50L, 60L, 40L, 
40L, 130L, 120L, 140L, 70L, 50L, 140L, 120L, 130L, 50L, 40L, 
80L, 140L, 100L, 60L, 70L, 50L, 60L, 60L, 130L, 40L, 130L, 100L, 
70L, 110L, 80L, 120L, 110L, 40L, 100L, 40L, 60L, 120L, 120L, 
70L, 80L, 130L, 60L, 100L, 100L, 60L, 70L, 90L, 100L, 140L, 70L, 
100L, 90L, 130L, 70L, 130L, 40L, 80L, 130L, 150L, 110L, 120L, 
140L, 90L, 60L, 90L, 80L, 120L, 150L, 90L, 150L, 50L, 50L, 100L, 
150L, 80L, 90L, 110L, 150L, 150L, 120L, 80L, 80L), gtangles = c(141L, 
58L, 44L, 154L, 120L, 90L, 128L, 147L, 147L, 120L, 127L, 66L, 
118L, 141L, 111L, 59L, 72L, 45L, 52L, 144L, 139L, 143L, 73L,  
59L, 148L, 141L, 135L, 63L, 51L, 88L, 147L, 110L, 68L, 78L, 63L, 
64L, 70L, 133L, 49L, 129L, 100L, 78L, 128L, 91L, 121L, 109L, 
48L, 113L, 50L, 68L, 135L, 120L, 85L, 97L, 136L, 59L, 112L, 103L, 
62L, 87L, 92L, 116L, 141L, 70L, 121L, 92L, 137L, 85L, 117L, 51L, 
84L, 128L, 162L, 102L, 127L, 151L, 115L, 57L, 93L, 92L, 117L, 
140L, 95L, 159L, 57L, 65L, 130L, 152L, 90L, 117L, 116L, 147L, 
140L, 116L, 98L, 95L), up = c(-1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
-1L, -1L, 1L, 1L, 1L, 1L, -1L, -1L, -1L, -1L, 1L, 1L, -1L, -1L, 
1L, 1L, -1L, 1L, 1L, -1L, 1L, 1L, 1L, 1L, 1L, -1L, -1L, 1L, 1L, 
1L, 1L, -1L, -1L, 1L, 1L, -1L, -1L, -1L, -1L, -1L, -1L, -1L, 
1L, -1L, -1L, -1L, -1L, -1L, 1L, -1L, 1L, 1L, -1L, -1L, -1L, 
-1L, 1L, -1L, 1L, -1L, -1L, -1L, 1L, -1L, 1L, -1L, 1L, 1L, 1L, 
-1L, -1L, -1L, -1L, -1L, -1L, 1L, -1L, 1L, 1L, -1L, -1L, 1L, 
1L, 1L, -1L, 1L, 1L, 1L)), .Names = c("tangle", "gtangles", "up"
), class = "data.frame", row.names = c(NA, -96L))

plotmeans2 <- function(br, alph) {
dt=br;   tmp   <- split(br$gtangles, br$tangle);   
means <- sapply(tmp, mean);  stdev <- sqrt(sapply(tmp, var));  
n <- sapply(tmp,length);  
ciw   <- qt(alph, n) * stdev / sqrt(n)
plotCI(x=means, uiw=ciw, col="black", barcol="blue", lwd=1,ylim=c(40,150),  xlim=c(1,12)); 
par(new=TRUE) dt= subset(br,up==1);   
tmp   <- split(dt$gtangles, dt$tangle);  
means <- sapply(tmp, mean);  
stdev <- sqrt(sapply(tmp, var));  
n <- sapply(tmp,length); 
ciw  <- qt(0.95, n) * stdev / sqrt(n)
plotCI(x=means, uiw=ciw, type='l',col="black", barcol="red", lwd=1,ylim=c(40,150), xlim=c(1,12),pch='+');
abline(v=6);abline(h=90);abline(30,10); par(new=TRUE);
dt=subset(br,up==-1);   
tmp <- split(dt$gtangles, dt$tangle);  
means <- sapply(tmp, mean);  
stdev <- sqrt(sapply(tmp, var));  
n <- sapply(tmp,length); 
ciw <- qt(0.95, n) * stdev / sqrt(n)
plotCI(x=means, uiw=ciw, type='l', col="black", barcol="blue",   lwd=1,ylim=c(40,150), xlim=c(1,12),pch='-');abline(v=6);abline(h=90);
abline(30,10);
}

plotmeans2(br,.95)
Adam SA
la source

Réponses:

21

Si vous êtes prêt à utiliser ggplot , vous pouvez essayer le code suivant.

Avec un prédicteur continu

library(ggplot2)
gp <- ggplot(data=br, aes(x=tangle, y=gtangles)) 
gp + geom_point() + stat_smooth(method="lm", fullrange=T) + facet_grid(. ~ up)

pour un graphique d'interaction à facettes

entrez la description de l'image ici

Pour un tracé d'interaction standard (comme celui produit par interaction.plot()), il vous suffit de supprimer la facette.

gp <- ggplot(data=br, aes(x=tangle, y=gtangles, colour=factor(up))) 
gp + geom_point() + stat_smooth(method="lm")

entrez la description de l'image ici

Avec un prédicteur discret

En utilisant l' ToothGrowthensemble de données (voir help(ToothGrowth)),

ToothGrowth$dose.cat <- factor(ToothGrowth$dose, labels=paste("d", 1:3, sep=""))
df <- with(ToothGrowth , aggregate(len, list(supp=supp, dose=dose.cat), mean))
df$se <- with(ToothGrowth , aggregate(len, list(supp=supp, dose=dose.cat), 
              function(x) sd(x)/sqrt(10)))[,3]

opar <- theme_update(panel.grid.major = theme_blank(),
                     panel.grid.minor = theme_blank(),
                     panel.background = theme_rect(colour = "black"))
gp <- ggplot(df, aes(x=dose, y=x, colour=supp, group=supp))
gp + geom_line(aes(linetype=supp), size=.6) + 
     geom_point(aes(shape=supp), size=3) + 
     geom_errorbar(aes(ymax=x+se, ymin=x-se), width=.1)
theme_set(opar)

entrez la description de l'image ici

chl
la source
Merci beaucoup pour la réponse détaillée. Je voulais demander, existe-t-il un moyen de faire des intervalles de confiance verticaux à chaque niveau de la variable indépendante? Existe-t-il un moyen de supprimer l'arrière-plan et de revenir au graphique «à l'ancienne»?
Adam SA
1
@Adam J'ai mis à jour ma réponse avec le cas de 2 variables catégorielles + une variable de réponse continue - j'espère que c'est ce que vous vouliez dire. J'ai également ajouté du code pour montrer comment personnaliser le ggplotthème. En règle générale, vous pouvez gp + theme_bw()simplement supprimer le fond gris; ici, j'ai également supprimé la grille.
chl
12

Il y a aussi le package d' effets de Fox et Hong dans R. Voir le J. Stat. Doux. articles ici et ici pour des exemples avec des intervalles de confiance et la génération de code R.

Ce n'est pas aussi joli qu'une solution ggplot, mais un peu plus général et une bouée de sauvetage pour les GLM modérément complexes.

conjugateprior
la source
1
(+1) Je dois admettre que je préfère cette approche :-)
chl
@chl et / ou Conjugate, pouvez-vous nous expliquer pourquoi vous préférez cette approche? Cela aiderait des gens comme moi à décider dans quelle méthode investir du temps.
Michael Bishop
1
@MichaelBishop Essentiellement parce qu'il encapsule beaucoup de choses délicates (traçage sur lien vs échelle de réponse, affichage de l'IC à 95% pour GLMMM, marginalisation par rapport aux termes d'interaction, etc.) qui seraient difficiles à gérer dans quelques commandes R (et personnellement, J'aime beaucoup les latticegraphismes :)
chl