pairwise lmm with kinship

gistfile1.txt

Documentation
X: matrix of predictors (including the intercept)
y: matrix of outcomes
kinship: the kinship matrix (with 1s on the diagonal, not 2s)
pwt_mat: matrix of pairwise sampling weights (ie, reciprocal of pairwise sampling probabilities)

Output:

1. genetic variance divided by residual variance
2. residual variance
3 and 4: beta

pairwise.R

## pairwise lmm
   
pairwise<-function(X,y, kinship, pwt_mat){
    n<-nrow(X)
    ii<-rep(1:n,each=n)
    jj<-rep(1:n, n)
    keep<-ii!=jj
    ii<-ii[keep]
    jj<-jj[keep]
    pwt<-pwt_mat[cbind(ii,jj)]
    
    beta<-c(0,0)
    s2<-var(y)
    
    devfun<-function(sigma2ratio){

        Xi<-diag(n)+kinship*sigma2ratio
     
        ## v11 is a vector of (1,1) entries of the matrix var(Y)
        ## for each pair, similarly for the others
        v11<-Xi[cbind(ii,ii)]
        v12<-Xi[cbind(ii,jj)]
        v22<-Xi[cbind(jj,jj)]
        ## explicit 2x2 determinants
        det<-v11*v22-v12*v12
        ## explicit 2x2 inverses
        inv11<- v22/det
        inv22<- v11/det
        inv12<- -v12/det

        ## X matrices for first and second element of each pair
        Xii<-X[ii,,drop=FALSE]
        Xjj<-X[jj,,drop=FALSE]

        ## X^TWX
        xtwx<- crossprod(Xii,pwt*inv11*Xii)+
            crossprod(Xjj,pwt*inv22*Xjj)+
            crossprod(Xii,pwt*inv12*Xjj)+
            crossprod(Xjj,pwt*inv12*Xii)

        ## X^WY
        xtwy<-crossprod(Xii,pwt*inv11*y[ii])+
            crossprod(Xjj,pwt*inv22*y[jj])+
            crossprod(Xii,pwt*inv12*y[jj])+
            crossprod(Xjj,pwt*inv12*y[ii])

        ## betahat at the given variance parameter values
        beta<<-solve(xtwx,xtwy)
        Xbeta<-X%*%beta

        ## two residuals per pair
        r<-y-Xbeta
        r1<-r[ii]
        r2<-r[jj]

        ## -2 times Gaussian log profile pairwise likelihood
        qf<-crossprod(r1,pwt*inv11*r1)+
            crossprod(r2,pwt*inv22*r2)+
            crossprod(r1,pwt*inv12*r2)+
            crossprod(r2,pwt*inv12*r1)

        Nhat<-sum(pwt)*2
        s2<<-qf/Nhat
        
        ##sum(log(det)*pwt) + qf 
        sum(log(det)*pwt) + Nhat*log(qf*2*pi/Nhat)
        
    }
    fit<-optimise(Vectorize(devfun), lower=0, upper=10)
    c(fit$minimum, s2, beta)
    
}