@mkitti - great suggestion - here is my R code and results now:
# Make an array
n = 20 ; nx = 400 ; v = 40 ; d = 3
array4d <- array( runif(n*nx*v*d ,-1,0),
dim = c(n, nx, v, d) )
#Assign to Julia
julia_assign("A", array4d)
microbenchmark:: microbenchmark(
res1 <- which_max_4D(array4d),
res2 <- julia_eval("whichmax4D1(A)", "R"),
res3 <- julia_eval("whichmax4D2(A)", "R"),
res4 <- julia_eval("whichmax4Dmap(A)", "R"),
res5 <- julia_eval("whichmax4Dreduce(A)", "R"),
res6 <- julia_eval("whichmax4Dtullio(A)", "R"),
res7 <- julia_eval("whichmax4D5_leandro(A)", "R"),
res8 <- julia_eval("whichmax4D5_elrod(A)", "R"),
res9 <- julia_eval("whichmax4D6(A)", "R"),
res10 <- whichmaxC(array4d), times=500)
Wow!! 
I also ran a benchmark on the julia_assign command:
Looks like you are right @mkitti - the slow part is moving the data from R to Julia
Edit: I ran it on an even bigger array also:
