Visualizing data spread along both y- and x-axis in 2D scatterplot I'm currently in the process of doing a literature review of a brain atrophy measurement, called BPF, commonly used in research regarding neurodegenerative disorders.
I have searched for publications having investigated BPF and for each study population I have come across I have noted population size, mean BPF, mean age and the standard deviations of these two means. I do not have access to source data.
I would like to visualize all these populations in a 2D bubble-plot (scatterplot with variable dot-sizes), with BPF on one axis and age on the other, each data point representing one population. I would like the size of the points (bubbles) to be representing population size. I would like to include error bars for each point representing +-SD of age (bar running along the age axis) and +-SD of BPF (bar running along the BPF-axis).  
I am currently using SPSS21 and have myself come so far as to be able to create the bubble-plot (using size(PopulationSize) for the dots) with error bar for one of the variables (using SPSS high-low-close graph), but I can not seem to get a second error bar along the other axis in there.
Any suggestions on how to accomplish making this graph, using SPSS or other software, would be greatly appreciated.
Thank you for your help!
 A: Here is a solution using R:

R Code:
#Make up data
age<-runif(10, 30,70)
agesd<-runif(10, 0.1,5)
bpf<-runif(10, 0,1)
bpfsd<-runif(10, 0.01,.2)    
pop.size<-runif(10,5,50)

#The plot
plot(age,bpf, pch=16, cex=log(pop.size), col=rainbow(length(pop.size)), 
ylim=c(0,1),xlim=c(20,90))
segments(age+agesd,bpf,age-agesd,bpf, lwd=2)
segments(age,bpf+bpfsd,age,bpf-bpfsd, lwd=2)    
legend("topright", legend=paste("Study",1:10), 
col=rainbow(length(pop.size)), pt.cex=1.5, pch=16)

A: If your errors are really orthogonal in both directions I would probably suggest using boxes instead of error bars (if they aren't ellipses (Friendly et al. 2013) are probably a better choice). Below the image is SPSS code illustrating how to plot boxes in GPL.

data list free / x y pop.
begin data
1 4 5
2 3 10
3 2 15
4 1 20
end data.
compute xlow = x - (x*.4)/sqrt(pop).
compute xhigh = x + (x*.4)/sqrt(pop).
compute ylow = y - (y*.4)/sqrt(pop).
compute yhigh = y + (y*.4)/sqrt(pop).
formats ALL (F1.0).

*Boxes.
GGRAPH
  /GRAPHDATASET NAME="graphdataset" VARIABLES=x y xlow xhigh ylow yhigh pop
  /GRAPHSPEC SOURCE=INLINE.
BEGIN GPL
 SOURCE: s=userSource(id("graphdataset"))
 DATA: x=col(source(s), name("x"))
 DATA: y=col(source(s), name("y"))
 DATA: pop=col(source(s), name("pop"))
 DATA: xlow=col(source(s), name("xlow"))
 DATA: xhigh=col(source(s), name("xhigh"))
 DATA: ylow=col(source(s), name("ylow"))
 DATA: yhigh=col(source(s), name("yhigh"))
 TRANS: casenum = index() 
 GUIDE: axis(dim(1), label("x"))
 GUIDE: axis(dim(2), label("y"))
 SCALE: pow(aesthetic(aesthetic.size)) 
 ELEMENT: polygon(position(link.hull((xlow + xhigh)*(ylow + yhigh))), color.interior(color.grey), split(casenum), transparency.interior(transparency."0.5"))
 ELEMENT: point(position(x*y), color.interior(color.black), size(pop))
END GPL.

I can show you how you would plot error bars in SPSS, but I think this is much clearer a picture. For scatterplots of more and overlapping cases things to keep in mind are to make the box elements semi-transparent so overlap is more clear. Also in the first graph plotting the point at the center of the rectangle may seem superflous, but this becomes more necessary as the area of the box shrinks, places that have higher certainty in the estimates shrink to nothing in the plot otherwise.
Below is a similar example extended to 246 observations.

set seed 10.
input program.
loop #i = 5 to 250.
  compute x = RV.UNIFORM(5,50).
  compute y = -.5*x + RV.NORMAL(0,10).
  compute pop = #i.
  end case.
end loop.
end file.
end input program.
compute xlow = x - (x*.4)/sqrt(pop).
compute xhigh = x + (x*.4)/sqrt(pop).
compute ylow = y - (y*.4)/sqrt(pop).
compute yhigh = y + (y*.4)/sqrt(pop).
formats ALL (F1.0).

*Boxes.
GGRAPH
  /GRAPHDATASET NAME="graphdataset" VARIABLES=x y xlow xhigh ylow yhigh pop
  /GRAPHSPEC SOURCE=INLINE.
BEGIN GPL
 SOURCE: s=userSource(id("graphdataset"))
 DATA: x=col(source(s), name("x"))
 DATA: y=col(source(s), name("y"))
 DATA: pop=col(source(s), name("pop"))
 DATA: xlow=col(source(s), name("xlow"))
 DATA: xhigh=col(source(s), name("xhigh"))
 DATA: ylow=col(source(s), name("ylow"))
 DATA: yhigh=col(source(s), name("yhigh"))
 TRANS: casenum = index() 
 GUIDE: axis(dim(1), label("x"))
 GUIDE: axis(dim(2), label("y"))
 SCALE: pow(aesthetic(aesthetic.size)) 
 ELEMENT: polygon(position(link.hull((xlow + xhigh)*(ylow + yhigh))), color.interior(color.grey), split(casenum), transparency.interior(transparency."0.5"))
 ELEMENT: point(position(x*y), color.interior(color.black), size(pop), transparency.interior(transparency."0.1"), transparency.exterior(transparency."1"))
END GPL.


Below is some code showing how to plot error bars instead of boxes in SPSS. Ellipses is more difficult, as you need some estimate of the correlation between the errors. Here I've posted some code that provides a function to draw the points for an ellipse in SPSS given the centroid, major and minor axis, and the degree of the slope. See this post on how to calculate an ellipse given a covariance matrix.
*Error Bars.
GGRAPH
  /GRAPHDATASET NAME="graphdataset" VARIABLES=x y xlow xhigh ylow yhigh
  /GRAPHSPEC SOURCE=INLINE.
BEGIN GPL
 SOURCE: s=userSource(id("graphdataset"))
 DATA: x=col(source(s), name("x"))
 DATA: y=col(source(s), name("y"))
 DATA: xlow=col(source(s), name("xlow"))
 DATA: xhigh=col(source(s), name("xhigh"))
 DATA: ylow=col(source(s), name("ylow"))
 DATA: yhigh=col(source(s), name("yhigh"))
 TRANS: casenum = index() 
 GUIDE: axis(dim(1), label("x"))
 GUIDE: axis(dim(2), label("y"))
 ELEMENT: edge(position((xlow + xhigh)*y)))
 ELEMENT: edge(position(x*(ylow + yhigh)))
 ELEMENT: point(position(x*y), color.interior(color.black))
END GPL.

