Bifactor Model in 3D

I was playing around with a Bifactor Model and found no elegant way to do it in 2D. So here is my attempt to do it in 3D:

My code in R:

library(rgl)
library(heplots)
vNorm<-function(x){sqrt(t(x)%*%x)}
vUnit<-function(a,b){(b-a)/vNorm(b-a)}

r3dDefaults$windowRect <- c(10, 40, 700, 700)

open3d()
nBarbs<-20
s1<-c(3,0,6)
s2<-c(12,0,6)
s3<-c(21,0,6)
g<-c(12,0,-6)
o<-c(12,-12, 6)
iDist<- c(0,0,-6)
iSpace<- c(3,0,0)
i1<-cbind(s1-iSpace+iDist,s1+iDist,s1+iSpace+iDist)
for (i in 1:3){shade3d( translate3d( cube3d(col="gray80"), i1[1,i],i1[2,i],i1[3,i]))}

i2<-cbind(s2-iSpace+iDist,s2+iDist,s2+iSpace+iDist)
for (i in 1:3){shade3d( translate3d( cube3d(col="gray60"), i2[1,i],i2[2,i],i2[3,i]))}

i3<-cbind(s3-iSpace+iDist,s3+iDist,s3+iSpace+iDist)
for (i in 1:3){shade3d( translate3d( cube3d(col="gray40"), i3[1,i],i3[2,i],i3[3,i]))}

spheres3d(s1,col="gray80",point_antialias=TRUE,smooth=TRUE)
spheres3d(s2,col="gray60",point_antialias=TRUE,smooth=TRUE)
spheres3d(s3,col="gray40",point_antialias=TRUE,smooth=TRUE)
for (i in 1:3){
arrow3d(s1,i1[,i]+c(0,0,1),color='gray80',n=nBarbs,barblen=0.2,lwd=2)
arrow3d(s2,i2[,i]+c(0,0,1),color='gray60',n=nBarbs,barblen=0.2,lwd=2)
arrow3d(s3,i3[,i]+c(0,0,1),color='gray40',n=nBarbs,barblen=0.2,lwd=2)
arrow3d(o,io[,i]+c(0,0,1),color='gray90',n=nBarbs,barblen=0.2,lwd=2)
}

spheres3d(g,col="black",point_antialias=T,smooth=T)
for (i in 0:8){
arrow3d(g,c(i*3,0,-1),color="black",barlen=0.05,n=nBarbs,barblen=0.15,lwd=2)
#   text3d(x=i*3,y=-1.3,0,paste0("T",i))
}

spheres3d(o,col="gray90",point_antialias=TRUE,smooth=TRUE)
io<-cbind(o-iSpace+iDist,o+iDist,o+iSpace+iDist)
for (i in 1:3){shade3d( translate3d( cube3d(col="gray80"), io[1,i],io[2,i],io[3,i]))}

arrow3d(s1,o-vUnit(s1,o),color='gray80',n=nBarbs,barblen=0.1,lwd=2)
arrow3d(s2,o-vUnit(s2,o),color='gray60',n=nBarbs,barblen=0.1,lwd=2)
arrow3d(s3,o-vUnit(s3,o),color='gray40',n=nBarbs,barblen=0.1,lwd=2)
arrow3d(g,o-vUnit(g,o),color='black',n=nBarbs,barblen=0.1,lwd=2)

text3d(s1+c(0,0,2),texts="S1",font=1,family="serif")
text3d(s2+c(0,0,2),texts="S2",font=1,family="serif")
text3d(s3+c(0,0,2),texts="S3",font=1,family="serif")
text3d(o+c(0,0,2),texts="Outcome",font=1,family="serif")
text3d(g+c(0,0,-2),texts="g",font=3,family="serif")
if (!rgl.useNULL())
play3d(spin3d(axis=c(0,0,1), rpm=10), duration=6)