Simes method in R

Question

I have data that looks like this: 3 column SNPs their gene based on Annovar and a p-value for every SNP. What I would like is to aggregate the p values for every gene.

snps <- data.frame(
  snp_id = c("rs1", "rs2", "rs3", "rs4", "rs5", "rs6", "rs7", "rs8"),
  gene = c("gene1", "gene1", "gene1", "gene1", "gene2", "gene2", "gene2", "gene2"),
  pvalue = c(0.7703884, 0.9648540, 0.9648540, 0.9648540, 0.54, 0.03, 0.03, 0.8)
)

what I need to do is :

1. For each gene, sort the p-values in ascending order.
2. For each p-value in the sorted list, multiply it by the number of SNPs in that gene.
3. Divide the result from step 2 by the rank of the p-value in the sorted list (i.e., its index in the list).
4. Take the minimum value from step 3 as the new p-value for that gene.
5. Adjust the p-values using the false discovery rate (FDR) correction method.

This is the formula for each gene :combined p-value $\min(n p_{(i)}/i)$ where $p_{(1)}, …, p_{(n)}$ are the ordered p-values.

I tried this code with the usage of simes.test package in R:

library(mppa)

x=c(0.77,0.964,0.964,0.964)
simes.test(x)

The answer is 1.285. This is the source code of the simes.test function:

function (x, returnstat = FALSE) 
{
    r = rank(x)
    T = min(length(x) * x/r)
    if (returnstat) 
        c(T, T)
    else T
}

How is it possible for the pvalue to be larger then 1? What is the problem? I would be happy if you could share your knowledge in this situation.