Online Trick

Here's an example how to conveniently summarize data with the cast function (package reshape). By the way you see how this could be done "in-conveniently" by hand. You also see how a for-loop works and how a matrix is constructed and filled. In addition this serves as an illustrative example how flexible "indexing" in R works, as seen in the below loop! (download data) (this example is adapted from https://stat.ethz.ch/pipermail/r-sig-ecology/2011-May/002174.html)

# Set path were you downloaded data to:
load("C:\\...\\Ecoregions_Taxa.RData")

# see what's in there:
ls()

# investigate data3:
str(data3)

# i want to know how many taxa are within each ecoregion -
# more presicly i want to know how many orders, families, genera are there within each region:

require(reshape)

dfm<- melt(data3, id = "ECO_NAME")
dfc<- cast(dfm, ECO_NAME~variable, function(x) length(unique(x)))

# the same by hand -
# make new variable of ECO_NAME*variable combinations:
dfm$variable2 <- as.factor(paste(dfm$variable, dfm$ECO_NAME, sep = " - "))

# make vector to collect results for all unique ECO_NAME*variable combinations:
Ns <- data.frame(count = rep(NA, length(unique(dfm$variable2))),
                 row.names = unique(dfm$variable2))

# loop through all unique ECO_NAME*variable combinations
# and record length of unique values:
for (i in levels(dfm$variable2)){
                        subset = dfm$value[dfm$variable2 == i]
                Ns[i, "count"] <- length(unique(subset)) 
                }

# put counts in matrix/table, Ns$count is in order of taxonomy
# levels (= "variable"), so I have to fill by cols (byrow = F),
# as the matrix/table colums are chosen to be the taxonomy levels:
result <- matrix(Ns$count,
                 nrow = length(levels(dfm$ECO_NAME)),
                 ncol = length(levels(dfm$variable)),
                 dimnames = list(levels(dfm$ECO_NAME), levels(dfm$variable)), byrow = F)
print(result)

To cite reshape in publications, please use:

### two samples with abundance data:
require(vegan)

### two samples with abundance data:
### 1st sample of counts (numbers representing species 1-15)
dat1 <- c(1,2,0,3,0,1,4,20,0,2,21,3,15,23,30)

### 2nd sample of counts (numbers representing species 1-15)               
dat2 <- c(0,1,0,1,0,0,1,10,0,0,29,1,22,30,25)
 
(div1=diversity(t(dat1),"shannon"))
(div2=diversity(t(dat2),"shannon"))
(rich1=sum(dat1>0))
(rich2=sum(dat2>0))
 
(tr.diff.div=abs(div1-div2))          ### observed difference
(tr.diff.rich=abs(rich1-rich2))       ### ...
 
K=2000
 
pop.diff.div <- pop.diff.rich <- rep(NA,K)      ### dataframe for null population of differences
pop.diff.div[1]=tr.diff.div
pop.diff.rich[1]=tr.diff.rich
 
for(i in 2:K){                        ### loop to generate null pop.diff. of differences
 
  ind1<-sum(dat1)                     ### sum of individuals sample no.1
  ind2<-sum(dat2)                     ### sum of ... no.2
  pool<-c(rep(1:length(dat1),dat1),   ### pooled sample with numbers representing species
          rep(1:length(dat2),dat2))   ### replicated as often as no. of individuals
 
  temp1=sample(pool,ind1,replace=T)   ### resample no.1
  temp2=sample(pool,ind2,replace=T)   ### resample no.2
 
  ### calculate diversity:
  div1.temp=diversity(t(tabulate(temp1)),"shannon")
  div2.temp=diversity(t(tabulate(temp2)),"shannon")
 
  rich1.temp=sum(tabulate(temp1)>0)
  rich2.temp=sum(tabulate(temp2)>0)
 
  pop.diff.div[i]=abs(div1.temp-div2.temp)
  pop.diff.rich[i]=abs(rich1.temp-rich2.temp)
}
 
(p.div=sum(pop.diff.div>=abs(tr.diff.div))/K)
(p.rich=sum(pop.diff.rich>=abs(tr.diff.rich))/K)

### diagramms to show null-distributions with obs. differences
par(mfrow=c(2,1))
hist(pop.diff.div)
abline(v=tr.diff.div,lty=3,col=2,lwd=2)

### diagramms to show null-distributions with obs. differences
hist(pop.diff.rich)
abline(v=tr.diff.rich,lty=3,col=2,lwd=2)

Taking into account what Jari commented below (simplifying code & fixing row and column totals - still having in mind that the test may be anti-conservative!)

dat1 <- c(1,2,0,3,0,1,4,20,0,2,21,3,15,23,30)
dat2 <- c(0,1,0,1,0,0,1,10,0,0,29,1,22,30,25)
  
(div1=diversity(t(dat1)))
(div2=diversity(t(dat2)))
  
dat <- rbind(dat1, dat2)

B = 2000

sim <- r2dtable(B - 1, rowSums(dat), colSums(dat))
div <- sapply(sim, diversity)
pop.diff.div <- abs(div[1, ] - div[2, ])

(p <- sum(pop.diff.div >= abs(div1 - div2)) / B)
hist(pop.diff.div); abline(v = abs(div1 - div2))

Now with null model and vegan::oecosimu

### two samples with abundance data:
### 1st sample of counts (numbers representing species 1-15)
dat1 <- c(1,2,0,3,0,1,4,20,0,2,21,3,15,23,30)
 
### 2nd sample of counts (numbers representing species 1-15)              
dat2 <- c(0,1,0,1,0,0,1,10,0,0,29,1,22,30,25)

dat <- rbind(dat1, dat2)

divdiff <- function(x) abs(diversity(x[1, ]) - diversity(x[2, ]))
oecosimu(dat, divdiff, "r2dtable", nsimul = 1999)