Hi,

Assume a table as below:

    X = 
    
            col1    col2    col3
    row1    "A"      "0"     "1"
    row2    "B"      "2"     "NA"
    row3    "C"      "1"     "2"


I select combinations of two rows, using the code below:

    pair <- apply(X, 2, combn, m=2)

This returns a matrix of the form:

    pair = 
    
         [,1] [,2] [,3]
    [1,] "A"  "0"  "1" 
    [2,] "B"  "2"  NA  
    [3,] "A"  "0"  "1" 
    [4,] "C"  "1"  "2" 
    [5,] "B"  "2"  NA  
    [6,] "C"  "1"  "2"


I wish to iterate over pair, taking two rows at a time, i.e. first isolate `[1,]` and `[2,]`, then `[3,]` and `[4,]` and finally, `[5,]` and `[6,]`. These rows will then be passed as arguments to regression models, i.e. `lm(Y ~ row[i]*row[j])`.

I am dealing with a large dataset. Can anybody advise how to iterate over a matrix two rows at a time, assign those rows to variables and pass as arguments to a function?

Thanks,  
S ;-)

**Edit:** In response to the comments, I should specify that my problem concerns SNP and expression data where I aim to do a pairwise multiple regression analysis (first order regression) in order to assess any possible SNP-SNP interactions that may effect the expression phenotype.
A simple idiom for "all pairs of a matrix" M of dim(N, C)

    for( i in 1:(N-1) ){
      for( j in (i+1):N ){
        print(paste(i,j))
      }
    }

I've had to do this eQTL analysis. For any non-trivial number of genotypes or expression measurements, the naive approach (e.g. building a new lm for each comparison) is too slow. If PHENO is a matrix (not data.frame) of expression values and GENO the set of genotypes for a single locus at all of the samples in PHENO, you want something like

    for each GENO
      model = lm( t(PHENO)~GENO )

Extend the linear model to match whatever interactions you want to test. You can extract the p-values of coefficients of a solved linear model with

    calculate.lm.pval=function(linear.model){
        # Given a solved linear model, return the p-values of the coefficients
        n.obs = dim(linear.model$residuals)[1]
        n.models = dim(linear.model$coefficients)[2]
        cc = matrix(coef(linear.model), nrow=linear.model$rank, ncol=n.models)
        se = calculate.linear.model.se(linear.model)
        t.stat = cc/se
        rdf = n.obs - linear.model$rank
        dm=data.frame( 2*pt(abs(t.stat), rdf, lower.tail = FALSE)  )
        names(dm) = as.vector(dimnames(linear.model$coefficients)[[2]])
        rownames(dm) = dimnames(linear.model$coefficients)[[1]]
        dm
    }