You could do a monte-carlo test of I>0:
First lets create a very correlated raster:
> r = raster(matrix(1:(50*50),50,50))
> Moran(r)
[1] 0.9694908
And now do 99 Moran's I of rasters that are random samplings of those values:
> M99 = sapply(1:99, function(i){v = r; v[]=sample(r[]);Moran(v)})
And let's see the distribution:
> hist(M99)
> range(M99)
[1] -0.02713834 0.02061854
And its clear that the Moran(r) is way outside the range of the simulations, so we can reject the null that the data are uncorrelated with respect to random sampling/rearrangement of the data.
To get the approximate pseudo p-value, see where the Moran stat for the data ranks amongst the simulations. Suppose the Moran stat for the data was 0.018 (using the 0.969 example from my code is a bit extreme), then compute the rank:
> rank(c(0.018, M99))[1]
[1] 98
Which shoes that 0.018 ranks 98 out of the 100 values (99 sims + itself). Hence reject H0 (no spatial autocorrelation under random arrangement hypothesis) with approximate p = 0.98.
If you do more simulations, then the general case is:
> R/(length(M99)+1)
[1] 0.98
for any number of simulations in M99.
This is what moran.mc from the spdep package does with data for polygons or other general neighbourhood structures.
Alternatively convert your raster to a grid of polygons and use the spdep functions with either 4-way or 8-way neighbours (or beyond...)
