2016-10-23 4 views
1

J'essaye de construire un simple graphique XY avec la production de lait (appelée FCM) de deux groupes différents de vaches (de la sortie que j'ai obtenue du modèle mixte, en utilisant les lsmeans et SE) . j'ai pu construire l'intrigue afficher les lsmeans en utilisant la fonction xyplot en treillis:Dessin SE dans xyplot avec les barres d'erreur

library(lattice)  
xyplot(lsmean~Time, type="b", group=Group, data=lsmeans2[order(lsmeans2$Time),], 
     pch=16, ylim=c(10,35), col=c("darkorange","darkgreen"), 
     ylab="FCM (kg/day)", xlab="Week", lwd=2, 
     key=list(space="top", 
       lines=list(col=c("darkorange","darkgreen"),lty=c(1,1),lwd=2), 
       text=list(c("Confinement Group","Pasture Group"), cex=0.8))) 

Je veux maintenant ajouter les barres d'erreur. J'ai essayé certaines choses avec la fonction panel.arrow, en copiant et en collant d'autres exemples, mais je n'ai rien trouvé de plus.

J'apprécierais vraiment de l'aide!

Mon jeu de données lsmeans2:

Group Time lsmean SE df lower.CL upper.CL 
Stall wk1 26.23299 0.6460481 59 24.19243 28.27356 
Weide wk1 25.12652 0.6701080 58 23.00834 27.24471 
Stall wk10 21.89950 0.6460589 59 19.85890 23.94010 
Weide wk10 18.45845 0.6679617 58 16.34705 20.56986 
Stall wk2 25.38004 0.6460168 59 23.33957 27.42050 
Weide wk2 22.90409 0.6679617 58 20.79269 25.01549 
Stall wk3 25.02474 0.6459262 59 22.98455 27.06492 
Weide wk3 24.05886 0.6679436 58 21.94751 26.17020 
Stall wk4 23.91630 0.6456643 59 21.87694 25.95565 
Weide wk4 22.23608 0.6678912 58 20.12490 24.34726 
Stall wk5 23.97382 0.6493483 59 21.92283 26.02481 
Weide wk5 18.14550 0.6677398 58 16.03480 20.25620 
Stall wk6 24.48899 0.6456643 59 22.44963 26.52834 
Weide wk6 19.40022 0.6697394 58 17.28319 21.51724 
Stall wk7 24.98107 0.6459262 59 22.94089 27.02126 
Weide wk7 19.71200 0.6677398 58 17.60129 21.82270 
Stall wk8 22.65167 0.6460168 59 20.61120 24.69214 
Weide wk8 19.35759 0.6678912 58 17.24641 21.46877 
Stall wk9 22.64381 0.6460481 59 20.60324 24.68438 
Weide wk9 19.26869 0.6679436 58 17.15735 21.38004 
+1

Pourriez-vous s'il vous plaît montrer 'dput (lsmeans2)'? – Christoph

Répondre

1

Y at-il une raison particulière que vous voulez utiliser xplot? ggplot2 est beaucoup plus facile à travailler et plus joli. Voici un exemple de ce que je pense que vous voulez.

#load ggplot2 
library(ggplot2) 

#load data 
d = structure(list(Group = structure(c(1L, 2L, 1L, 2L, 1L, 2L, 1L, 
2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L), .Label = c("Stall", 
"Weide"), class = "factor"), Time = structure(c(1L, 1L, 2L, 2L, 
3L, 3L, 4L, 4L, 5L, 5L, 6L, 6L, 7L, 7L, 8L, 8L, 9L, 9L, 10L, 
10L), .Label = c("wk1", "wk10", "wk2", "wk3", "wk4", "wk5", "wk6", 
"wk7", "wk8", "wk9"), class = "factor"), lsmean = c(26.23299, 
25.12652, 21.8995, 18.45845, 25.38004, 22.90409, 25.02474, 24.05886, 
23.9163, 22.23608, 23.97382, 18.1455, 24.48899, 19.40022, 24.98107, 
19.712, 22.65167, 19.35759, 22.64381, 19.26869), SE = c(0.6460481, 
0.670108, 0.6460589, 0.6679617, 0.6460168, 0.6679617, 0.6459262, 
0.6679436, 0.6456643, 0.6678912, 0.6493483, 0.6677398, 0.6456643, 
0.6697394, 0.6459262, 0.6677398, 0.6460168, 0.6678912, 0.6460481, 
0.6679436), df = c(59L, 58L, 59L, 58L, 59L, 58L, 59L, 58L, 59L, 
58L, 59L, 58L, 59L, 58L, 59L, 58L, 59L, 58L, 59L, 58L), lower.CL = c(24.19243, 
23.00834, 19.8589, 16.34705, 23.33957, 20.79269, 22.98455, 21.94751, 
21.87694, 20.1249, 21.92283, 16.0348, 22.44963, 17.28319, 22.94089, 
17.60129, 20.6112, 17.24641, 20.60324, 17.15735), upper.CL = c(28.27356, 
27.24471, 23.9401, 20.56986, 27.4205, 25.01549, 27.06492, 26.1702, 
25.95565, 24.34726, 26.02481, 20.2562, 26.52834, 21.51724, 27.02126, 
21.8227, 24.69214, 21.46877, 24.68438, 21.38004)), .Names = c("Group", 
"Time", "lsmean", "SE", "df", "lower.CL", "upper.CL"), class = "data.frame", row.names = c(NA, 
-20L)) 

#fix week 
library(stringr) 
library(magrittr) 
d$Time %<>% as.character() %>% str_replace(pattern = "wk", replacement = "") %>% as.numeric() 

#plot 
ggplot(d, aes(Time, lsmean, color = Group, group = Group)) + 
    geom_point() + 
    geom_errorbar(aes(ymin = lower.CL, ymax = upper.CL), width = .2) + 
    geom_line() + 
    ylim(10, 35) + 
    scale_x_continuous(name = "Week", breaks = 1:10) + 
    ylab("FCM (kg/day)") + 
    scale_color_discrete(label = c("Confinement Group","Pasture Group")) 

enter image description here

+0

Merci beaucoup! L'heure n'est pas dans le bon ordre, mais "useR" a résolu cela dans son exemple. Je suis très reconnaissant pour la réponse rapide! –

+0

J'ai mis à jour le code pour corriger la commande de la semaine. Notez que je viens d'utiliser les données que vous m'avez données et que vos données ont été hors service. J'ai enlevé "wk" des chaînes parce que c'était redondant; l'unité de temps est déjà indiquée par l'étiquette X. – Deleet

+0

Vous avez raison. Je vous remercie! –

0

Pour être complet, voici une solution en utilisant xyplot:

# Reproducible data 
lsmeans2 = structure(list(Group = structure(c(1L, 2L, 1L, 2L, 1L, 2L, 1L, 
2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L), .Label = c("Stall", 
"Weide"), class = "factor"), Time = structure(c(1L, 1L, 2L, 2L, 
3L, 3L, 4L, 4L, 5L, 5L, 6L, 6L, 7L, 7L, 8L, 8L, 9L, 9L, 10L, 
10L), .Label = c("wk1", "wk10", "wk2", "wk3", "wk4", "wk5", "wk6", 
"wk7", "wk8", "wk9"), class = "factor"), lsmean = c(26.23299, 
25.12652, 21.8995, 18.45845, 25.38004, 22.90409, 25.02474, 24.05886, 
23.9163, 22.23608, 23.97382, 18.1455, 24.48899, 19.40022, 24.98107, 
19.712, 22.65167, 19.35759, 22.64381, 19.26869), SE = c(0.6460481, 
0.670108, 0.6460589, 0.6679617, 0.6460168, 0.6679617, 0.6459262, 
0.6679436, 0.6456643, 0.6678912, 0.6493483, 0.6677398, 0.6456643, 
0.6697394, 0.6459262, 0.6677398, 0.6460168, 0.6678912, 0.6460481, 
0.6679436), df = c(59L, 58L, 59L, 58L, 59L, 58L, 59L, 58L, 59L, 
58L, 59L, 58L, 59L, 58L, 59L, 58L, 59L, 58L, 59L, 58L), lower.CL = c(24.19243, 
23.00834, 19.8589, 16.34705, 23.33957, 20.79269, 22.98455, 21.94751, 
21.87694, 20.1249, 21.92283, 16.0348, 22.44963, 17.28319, 22.94089, 
17.60129, 20.6112, 17.24641, 20.60324, 17.15735), upper.CL = c(28.27356, 
27.24471, 23.9401, 20.56986, 27.4205, 25.01549, 27.06492, 26.1702, 
25.95565, 24.34726, 26.02481, 20.2562, 26.52834, 21.51724, 27.02126, 
21.8227, 24.69214, 21.46877, 24.68438, 21.38004)), .Names = c("Group", 
"Time", "lsmean", "SE", "df", "lower.CL", "upper.CL"), class = "data.frame", row.names = c(NA, 
-20L)) 

xyplot(lsmean~Time, type="b", group=Group, data=lsmeans2[order(lsmeans2$Time),], 
     panel = function(x, y, ...){ 
     panel.arrows(x, y, x, lsmeans2$upper.CL, length = 0.15, 
         angle = 90, col=c("darkorange","darkgreen")) 
     panel.arrows(x, y, x, lsmeans2$lower.CL, length = 0.15, 
         angle = 90, col=c("darkorange","darkgreen")) 
     panel.xyplot(x,y, ...) 
     }, 
     pch=16, ylim=c(10,35), col=c("darkorange","darkgreen"), 
     ylab="FCM (kg/day)", xlab="Week", lwd=2, 
     key=list(space="top", 
       lines=list(col=c("darkorange","darkgreen"),lty=c(1,1),lwd=2), 
       text=list(c("Confinement Group","Pasture Group"), cex=0.8))) 

L'argument longueur panel.arrows change la largeur des têtes d'erreur. Vous pouvez jouer avec ce paramètre pour obtenir une largeur que vous aimez.

xyplot with wrong ordering

Notez que même si vous aviez lsmeans2[order(lsmeans2$Time),] lorsque vous spécifiez le data =, l'ordre du temps est toujours mauvais. C'est parce que le temps est un facteur, et que R ne sait pas que vous voulez le commander par le suffixe numérique de wk. Cela signifie, qu'il triera wk10 avant WK2, parce que 1 est inférieur à 2. Vous pouvez utiliser cette petite astuce ci-dessous pour commander correctement:

# Order first by the character lenght, then by Time 
Timelevels = levels(lsmeans2$Time) 
Timelevels = Timelevels[order(nchar(Timelevels), Timelevels)] 

# Reorder the levels 
lsmeans2$Time = factor(lsmeans2$Time, levels = Timelevels) 

# Create Subset 
lsmeansSub = lsmeans2[order(lsmeans2$Time),] 

xyplot(lsmean~Time, type="b", group=Group, data=lsmeansSub, 
     panel = function(x, y, yu, yl, ...){ 
     panel.arrows(x, y, x, lsmeansSub$upper.CL, length = 0.15, 
         angle = 90, col=c("darkorange","darkgreen")) 
     panel.arrows(x, y, x, lsmeansSub$lower.CL, length = 0.15, 
         angle = 90, col=c("darkorange","darkgreen")) 
     panel.xyplot(x, y, ...) 
     }, 
     pch=16, ylim=c(10,35), col=c("darkorange","darkgreen"), 
     ylab="FCM (kg/day)", xlab="Week", lwd=2, 
     key=list(space="top", 
       lines=list(col=c("darkorange","darkgreen"),lty=c(1,1),lwd=2), 
       text=list(c("Confinement Group","Pasture Group"), cex=0.8))) 

Notez que même après réordonner les niveaux de « temps », J'ai encore besoin d'utiliser les données triées pour l'argument data =. En effet, xyplot trace les points dans l'ordre qui apparaît dans l'ensemble de données, et non dans l'ordre des niveaux de facteur.

xyplot with correct ordering

+0

Merci beaucoup! J'apprends encore et cela m'a beaucoup aidé! –

+0

@MelanieS Si vous pensez qu'une réponse est utile, veuillez l'accepter afin que les autres puissent la voir. – useR