# How to replicate a plot of means for a 2 by 3 by 4 design in R?

I have a graph that I did in SPSS

I would like to replicate this in R.

The data file is here: http://dl.dropbox.com/u/22681355/sendergraphR.csv

UPDATE:

I've figured out how to to it for the other graph:

graph2<-read.csv(file="SendergraphR.csv")
ddr <- recast(graph2,no.GREEN+no.RED+Senderidentity~variable,
fun.aggregate=mean,id.var=c("no.GREEN","no.RED","Senderidentity"))

qplot(x=no.GREEN,y=Mean_Message,data=ddr,colour=Senderidentity,
group=Senderidentity,geom="line")+facet_wrap(~no.RED,ncol=1)


Now what I would like to do is to separate this second graph into two columns one looking at the cases where the variable urn1 is blue and the other where its red: here's a graph of how it should look:

I was thinking of doing it the following way:

ddr <- recast(graph2,no.GREEN+no.RED+Senderidentity+urn1~variable,
fun.aggregate=mean,id.var=c("no.GREEN","no.RED","Senderidentity", "urn1"))
qplot(x=no.GREEN,y=Mean_Message,data=ddr,colour=Senderidentity, group=Senderidentity,geom="line")+facet_wrap(~no.RED,ncol=1) +
facet_wrap(~urn1, ncol=2)


But this doesn't seem to work. What am I doing wrong?

• There are two problems here. One to calculate the mean probabilities and means, another to plot the results. First problem requires some model, which you do not specify. So it is hard to help you, since we need to guess the model. The second should be easy to do, using ggplot2. Please supply the response variables of your models, then it would be possible to recreate these graphs. – mpiktas Aug 22 '11 at 10:41
• @mpiktas on the graph Mean probability of guessing RED is just a bad way of looking at the variable DecisionasReceiver. Basically we're interested in how many chose 1 under different conditions. Similarly on the second graph Mean message just means how many people chose 1 in the DecisionasSender variable. Both datasets contain all the information. Here's how I did it SPSS: dl.dropbox.com/u/22681355/graphforreceiver.tiff dl.dropbox.com/u/22681355/sender-graph.tiff These should give you an idea of how its set up. The dataset is the same. THanks! – upabove Aug 22 '11 at 10:47
• at this point this seems not to be a statistics question--- migrate to SO R tag? – John Aug 22 '11 at 11:52
• I would suggest to keep this question here as it already got valiant follow-up and an accepted answer (not that I don't want to close or migrate the 4th Rish post of the day). – chl Aug 22 '11 at 15:26
• @Daniel, please revert this to the previous question, and ask the new question separately, linking to this one. Also please see the FAQ and properly format the code, now it has an error in it. The answer to your second question is very simple, but you must figure it out by yourself. This site is for learning, not for getting your job done by somebody else. If somebody does not show inclination to learn (as you do), there is no incentive to help. – mpiktas Aug 23 '11 at 7:19

This is more or less standard ggplot2 problem. So I will only give an idea how to reproduce this graph. First get the data

dd<-structure(list(no.GREEN=c(1L,1L,2L,2L,3L,3L,4L,4L,1L,1L,2L,2L,3L,3L,4L,4L,1L,1L,2L,2L,3L,3L,4L,4L,1L,1L,2L,2L,3L,3L,4L,4L,1L,1L,2L,2L,3L,3L,4L,4L,1L,1L,2L,2L,3L,3L,4L,4L,1L,1L,2L,2L,3L,3L,4L,4L,1L,1L,2L,2L,3L,3L,4L,4L,1L,1L,2L,2L,3L,3L,4L,4L,1L,1L,2L,2L,3L,3L,4L,4L,1L,1L,2L,2L,3L,3L,4L,4L,1L,1L,2L,2L,3L,3L,4L,4L,1L,1L,2L,2L,3L,3L,4L,4L,1L,1L,2L,2L,3L,3L,4L,4L,1L,1L,2L,2L,3L,3L,4L,4L,1L,1L,2L,2L,3L,3L,4L,4L,1L,1L,2L,2L,3L,3L,4L,4L,1L,1L,2L,2L,3L,3L,4L,4L,1L,1L,2L,2L,3L,3L,4L,4L,1L,1L,2L,2L,3L,3L,4L,4L,1L,1L,2L,2L,3L,3L,4L,4L,1L,1L,2L,2L,3L,3L,4L,4L,1L,1L,2L,2L,3L,3L,4L,4L,1L,1L,2L,2L,3L,3L,4L,4L,1L,1L,2L,2L,3L,3L,4L,4L,1L,1L,2L,2L,3L,3L,4L,4L,1L,1L,2L,2L,3L,3L,4L,4L,1L,1L,2L,2L,3L,3L,4L,4L,1L,1L,2L,2L,3L,3L,4L,4L,1L,1L,2L,2L,3L,3L,4L,4L,1L,1L,2L,2L,3L,3L,4L,4L,1L,1L,2L,2L,3L,3L,4L,4L,1L,1L,2L,2L,3L,3L,4L,4L,1L,1L,2L,2L,3L,3L,4L,4L,1L,1L,2L,2L,3L,3L,4L,4L,1L,1L,2L,2L,3L,3L,4L,4L,1L,1L,2L,2L,3L,3L,4L,4L,1L,1L,2L,2L,3L,3L,4L,4L,1L,1L,2L,2L,3L,3L,4L,4L,1L,1L,2L,2L,3L,3L,4L,4L,1L,1L,2L,2L,3L,3L,4L,4L),no.RED=c(5L,5L,5L,5L,5L,5L,5L,5L,5L,5L,5L,5L,5L,5L,5L,5L,5L,5L,5L,5L,5L,5L,5L,5L,5L,5L,5L,5L,5L,5L,5L,5L,5L,5L,5L,5L,5L,5L,5L,5L,5L,5L,5L,5L,5L,5L,5L,5L,5L,5L,5L,5L,5L,5L,5L,5L,5L,5L,5L,5L,5L,5L,5L,5L,5L,5L,5L,5L,5L,5L,5L,5L,5L,5L,5L,5L,5L,5L,5L,5L,5L,5L,5L,5L,5L,5L,5L,5L,5L,5L,5L,5L,5L,5L,5L,5L,5L,5L,5L,5L,5L,5L,5L,5L,5L,5L,5L,5L,5L,5L,5L,5L,5L,5L,5L,5L,5L,5L,5L,5L,5L,5L,5L,5L,5L,5L,5L,5L,5L,5L,5L,5L,5L,5L,5L,5L,5L,5L,5L,5L,5L,5L,5L,5L,5L,5L,5L,5L,5L,5L,5L,5L,7L,7L,7L,7L,7L,7L,7L,7L,7L,7L,7L,7L,7L,7L,7L,7L,7L,7L,7L,7L,7L,7L,7L,7L,7L,7L,7L,7L,7L,7L,7L,7L,7L,7L,7L,7L,7L,7L,7L,7L,7L,7L,7L,7L,7L,7L,7L,7L,7L,7L,7L,7L,7L,7L,7L,7L,7L,7L,7L,7L,7L,7L,7L,7L,7L,7L,7L,7L,7L,7L,7L,7L,7L,7L,7L,7L,7L,7L,7L,7L,9L,9L,9L,9L,9L,9L,9L,9L,9L,9L,9L,9L,9L,9L,9L,9L,9L,9L,9L,9L,9L,9L,9L,9L,9L,9L,9L,9L,9L,9L,9L,9L,9L,9L,9L,9L,9L,9L,9L,9L,9L,9L,9L,9L,9L,9L,9L,9L,9L,9L,9L,9L,9L,9L,9L,9L,9L,9L,9L,9L,9L,9L,9L,9L,9L,9L,9L,9L,9L,9L,9L,9L,9L,9L,9L,9L,9L,9L,9L,9L,9L,9L,9L,9L,9L,9L,9L,9L,9L,9L,9L,9L,9L,9L,9L,9L),Messagereceived=structure(c(2L,1L,2L,1L,2L,1L,2L,1L,2L,1L,2L,1L,2L,1L,2L,1L,2L,1L,2L,1L,2L,1L,2L,1L,2L,1L,2L,1L,2L,1L,2L,1L,2L,1L,2L,1L,2L,1L,2L,1L,2L,1L,2L,1L,2L,1L,2L,1L,2L,1L,2L,1L,2L,1L,2L,1L,2L,1L,2L,1L,2L,1L,2L,1L,2L,1L,2L,1L,2L,1L,2L,1L,2L,1L,2L,1L,2L,1L,2L,1L,2L,1L,2L,1L,2L,1L,2L,1L,2L,1L,2L,1L,2L,1L,2L,1L,2L,1L,2L,1L,2L,1L,2L,1L,2L,1L,2L,1L,2L,1L,2L,1L,2L,1L,2L,1L,2L,1L,2L,1L,2L,1L,2L,1L,2L,1L,2L,1L,2L,1L,2L,1L,2L,1L,2L,1L,2L,1L,2L,1L,2L,1L,2L,1L,2L,1L,2L,1L,2L,1L,2L,1L,2L,1L,2L,1L,2L,1L,2L,1L,2L,1L,2L,1L,2L,1L,2L,1L,2L,1L,2L,1L,2L,1L,2L,1L,2L,1L,2L,1L,2L,1L,2L,1L,2L,1L,2L,1L,2L,1L,2L,1L,2L,1L,2L,1L,2L,1L,2L,1L,2L,1L,2L,1L,2L,1L,2L,1L,2L,1L,2L,1L,2L,1L,2L,1L,2L,1L,2L,1L,2L,1L,2L,1L,2L,1L,2L,1L,2L,1L,2L,1L,2L,1L,2L,1L,2L,1L,2L,1L,2L,1L,2L,1L,2L,1L,2L,1L,2L,1L,2L,1L,2L,1L,2L,1L,2L,1L,2L,1L,2L,1L,2L,1L,2L,1L,2L,1L,2L,1L,2L,1L,2L,1L,2L,1L,2L,1L,2L,1L,2L,1L,2L,1L,2L,1L,2L,1L,2L,1L,2L,1L,2L,1L,2L,1L,2L,1L,2L,1L,2L,1L,2L,1L,2L,1L,2L,1L,2L,1L,2L,1L,2L,1L,2L,1L,2L,1L,2L,1L,2L,1L,2L,1L,2L,1L,2L,1L),.Label=c("blue","red"),class="factor"),Decisionasreceivercode=c(0L,1L,0L,1L,1L,0L,1L,0L,1L,1L,1L,1L,1L,0L,1L,0L,0L,1L,0L,1L,1L,0L,1L,0L,0L,0L,0L,0L,0L,0L,0L,0L,0L,0L,0L,0L,0L,0L,0L,0L,0L,1L,0L,1L,1L,0L,1L,0L,0L,1L,0L,1L,1L,0L,1L,0L,0L,1L,0L,1L,1L,0L,1L,0L,0L,0L,0L,0L,1L,0L,1L,0L,0L,1L,0L,0L,1L,0L,1L,0L,0L,1L,0L,0L,1L,0L,1L,0L,0L,1L,1L,0L,1L,0L,1L,0L,0L,1L,1L,1L,1L,0L,1L,0L,0L,1L,0L,1L,1L,0L,1L,0L,0L,1L,0L,1L,1L,0L,1L,0L,1L,0L,1L,0L,1L,0L,1L,0L,1L,1L,0L,0L,1L,1L,0L,0L,0L,1L,0L,0L,0L,1L,1L,1L,0L,1L,1L,0L,1L,0L,1L,0L,0L,1L,0L,1L,1L,0L,1L,0L,0L,1L,0L,1L,1L,0L,1L,0L,0L,1L,0L,1L,1L,0L,1L,0L,0L,1L,0L,0L,1L,1L,0L,1L,1L,0L,1L,0L,1L,0L,1L,0L,1L,1L,0L,0L,1L,0L,1L,0L,1L,1L,0L,0L,1L,1L,1L,0L,1L,0L,0L,1L,1L,0L,1L,0L,0L,1L,0L,1L,0L,0L,1L,1L,1L,1L,1L,0L,1L,0L,1L,0L,1L,1L,1L,1L,1L,1L,1L,0L,1L,1L,1L,1L,1L,1L,1L,1L,1L,1L,1L,1L,1L,1L,1L,1L,0L,1L,0L,1L,1L,0L,1L,0L,1L,1L,1L,1L,1L,0L,1L,0L,0L,0L,1L,0L,1L,0L,1L,0L,1L,1L,1L,0L,1L,0L,1L,0L,1L,1L,1L,1L,1L,1L,0L,1L,1L,0L,0L,1L,1L,0L,0L,1L,0L,1L,0L,0L,1L,0L,0L,0L,1L,0L,0L,0L,1L,0L,0L,1L,0L,1L,1L,1L,1L,0L,1L,0L),OptimalResponse=c(1L,0L,1L,0L,1L,0L,1L,0L,1L,0L,1L,0L,1L,0L,1L,0L,1L,0L,1L,0L,1L,0L,1L,0L,1L,0L,1L,0L,1L,0L,1L,0L,1L,0L,1L,0L,1L,0L,1L,0L,1L,0L,1L,0L,1L,0L,1L,0L,1L,0L,1L,0L,1L,0L,1L,0L,1L,0L,1L,0L,1L,0L,1L,0L,1L,0L,1L,0L,1L,0L,1L,0L,1L,0L,1L,0L,1L,0L,1L,0L,1L,0L,1L,0L,1L,0L,1L,0L,1L,0L,1L,0L,1L,0L,1L,0L,1L,0L,1L,0L,1L,0L,1L,0L,1L,0L,1L,0L,1L,0L,1L,0L,1L,0L,1L,0L,1L,0L,1L,0L,1L,0L,1L,0L,1L,0L,1L,0L,1L,0L,1L,0L,1L,0L,1L,0L,1L,0L,1L,0L,1L,0L,1L,0L,1L,0L,1L,0L,1L,0L,1L,0L,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,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("no.GREEN","no.RED","Messagereceived","Decisionasreceivercode","OptimalResponse"),class="data.frame",row.names=c(NA,-328L))


Then transform it to get the desired statistics:

ddr<-recast(dd,no.GREEN+no.RED+Messagereceived~variable,fun.aggregate=mean,id.var=c("no.GREEN","no.RED","Messagereceived"))


And finaly plot it:

qplot(x=no.GREEN,y=Decisionasreceivercode,data=ddr,colour=no.RED,group=no.RED,geom="line")+facet_wrap(~Messagereceived,ncol=1)


The result will be the following:

I will leave the second graph and the cleaning up of the first one as an exercise :)

Update: There is alternative cleaner way to get the aggregation:

ddr<-aggregate(dd\$Decision,by=as.list(dd[,1:3]),mean)


We lose the name of the variable this way, but save some time with not typing three names two times.

• What? You're not doing all his work for him? +1 for that (well, that and the excellent answer)! – Nick Sabbe Aug 22 '11 at 11:51
• @mpiktas thank you very much! there definitely remains quite a few things to clean up, but this is an excellent lead. My only remaining question is how I can add an additional line with the Optimal variable plotted. thanks again – upabove Aug 22 '11 at 12:12
• @mpiktas I'm trying to reproduce the graph to have only 1,2,3,4 in the bottom. I've added xlim(1,4) but it won't change it, what am I missing? – upabove Aug 22 '11 at 13:05
• @Daniel, I suggest to read ggplot2 book. xlim only sets the limits, they are fine here. What you want is to change the type of variable. qplot assumes that no.GREEN is numeric, hence it produces numeric axis. Try changing it to factor and then inspect the results. – mpiktas Aug 22 '11 at 13:19
• @mpiktas I do want to read the ggplot2 book in the future, but currently I do not have it. I could only download the chapter relating to qplot but it didn't really detail what I should do. I've tried changing it to numeric by adding .class="numeric"but an error always comes up like Error: unexpected ',' in "L," – upabove Aug 22 '11 at 15:12

You could also use the lattice library:

A simple version (not shown):

library(lattice)
group=no.RED, type="a", data=dd)


Or with better labeling:

names(dd)[3] <- 'Message Received'
xyplot(Decisionasreceivercode ~ factor(no.GREEN) | Message Received,
group=no.RED, type="a", xlab="No. Green",

• See import.csv to import your file into R. The L means only that it's an integer; import.csv will automatically decide how best to store it (the defaults can be changed if you like). mpitkas used dput to output the data in a way that could easily be used by others, but that's usually something you would do only to help others recreate your code. – Aaron - Reinstate Monica Aug 22 '11 at 14:59
• @Aaron, you probably meant read.csv? Quick search on rseek.org did not return any hits for import.csv. – mpiktas Aug 23 '11 at 7:09
• Yes, read.csv. Thanks for the fix -- perhaps I was thinking about PROC IMPORT and got mixed up. – Aaron - Reinstate Monica Aug 23 '11 at 15:24