tslumley

Shows the power of the Rao-Scott (Wald) test in a 2-d problem.  
Alternative is parametrised by displacement (fixed to give specified power for intrinsic test) and by angle theta.
Purple circle indicates power for intrinsic test (constant)
black points are corresponding power for RS-type test. 

Maybe facet this to get a three-d problem?

tslumley

stanfile<-write_stan_file(
"data {
  int<lower=2> K;
  int<lower=0> N;
  int<lower=1> D;
  array[N] int<lower=1,upper=K> ybaseline;
    array[N] int<lower=ybaseline, upper=K> y;
  array[N] row_vector[D] x;

tslumley

poll<-read.table("pollen.txt")
pairs(poll, col="#00000010",pch=".")
plot(poll[,1:2], col="#00000060",pch=".")

tslumley

library(survey)

population<-data.frame(y=1:1000,
	strata=rep(1:3,c(100,200,700)),
	popsizes=rep(c(100,200,700),c(100,200,700))
	)
with(population, by(y,strata,var))

sampsizes<-c("1"=10,"2"=30,"3"=50)	

tslumley

abilities
academics
accessing
accessory
activated
addiction
additions
admission
advantage
aerospace

tslumley

paired_test<-function(formula, data, subset, paired=NULL, ...){
	
	if (is.null(paired))
		return("do the current thing")
	
	## make sure paired is ~id
	
	if (!is.language(paired))
		stop("paired must be a formula")
	if(!(length(paired)==2 && paired[[1L]]=="~" && length(paired[[2L]])<=2))

tslumley

## New sample.int

sample_int<-function(n, size=NULL, replace=FALSE, prob=NULL,
                     useHash=(n >  1e+07 && !replace && is.null(prob) && (!is.null(size)) && size <= n/2),
                     method=c("sequential","marginal","poisson")){
    if (replace || is.null(prob)){
        if (is.null(size)){
            size<-n
        }
    } else{

tslumley

void tille_incprob(double a[], int *pn, int *plen){

     double a_sum=0;
     int i,n,len,l,l1;

     n=*pn;
     len=*plen;
     for(i=0;i<len; i++){
       a_sum+=a[i];
     }

tslumley

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

tslumley

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)