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)

tslumley

mo day cancels
11 2 1806
11 1 1846
10 31 1887
10 30 1276
10 29 914
10 28 1828
10 27 1672
10 26 1686
10 25 1845

tslumley

lagged<-function(x,lag=1){
	if (lag==0) return(x)
        n<-length(x)
        c(rep(NA,lag),x[-( (n-lag+1):n)])
        }

pdlweights<-function(lag,degree,tiedown=c(F,F)){
	if (any(tiedown)) stop("Tiedown not working")
        contr.poly(lag,contrasts=F)[,tiedown[1]+(1:degree)]