Based on Brant's test for the proportional odds assumption

brant_test.R

olr_brant_test<-function(formula, design,test=c("brant-original","omnidirectional-Wald")){
    test<-match.arg(test)
    m1<-svyolr(formula, design = design)
    K<-length(m1$lev)
    P<-length(m1$coef)
    
    get_infl<-function(k,formula,design){
        y<-formula[[2]]
        formula[[2]]<-bquote(I(as.numeric(.(y))>.(k)))
        mk<-svyglm(formula, design, family=quasibinomial, influence=TRUE)
        list(coef(mk), attr(mk,"influence"))
    }
    
    fits<-lapply(1:(K-1), function(k) get_infl(k, formula, design))
    infs<-do.call(cbind, lapply(fits, "[[",2))
    combined_V<- vcov(svytotal(infs/weights(design), design))
    coefs<-lapply(fits,"[[",1)
    
    
    B<-do.call(c,lapply(coefs, "[",-1))
    VarB<-(combined_V[-(1+(0:(K-2))*(P+1)),-(1+(0:(K-2))*(P+1))])
    VarBinv<-solve(VarB)
    
    if (test=="brant-original"){
        D<-diag(P)	
        for(k in 2:(K-1)){
            D<-rbind(D,diag(P))
        }
        D<-cbind(D,matrix(0,nrow=(P)*(K-1),ncol=(K-2)))
        for (k in 1:(K-2)){
            D[(P-1)*k+(1:(P)),P+k]<-coefs[[k+1]][-1]
        }       
        
        deltahat<- solve(t(D)%*%VarBinv%*%D, t(D)%*%VarBinv%*%B)
        VarDelta<-solve(t(D)%*%VarBinv%*%D)
        
        i<-P+(1:(K-2))
        brantTest<-deltahat[i]%*%solve(VarDelta[i,i])%*%deltahat[i]
        
        list(X2=brantTest, p=pchisq(brantTest, df=(K-2),lower.tail=FALSE), df=K-2, test=test)
    } else {
        D<-cbind(1,-diag(K-2))%x%diag(P)
        DB<-D%*%B
        bigTest<-t(DB)%*%solve(D%*%VarB%*%t(D))%*%DB
        list(X2=bigTest, p=pchisq(bigTest, df=(K-2)*P,lower.tail=FALSE), df=(K-2)*P, test=test)
        }


    }