tslumley

counts<-c(173,291,320,407,311,292,341,321,188)
times<-c(as.numeric(as.Date("2013-12-31")-as.Date("2013-4-1")),365,365,366,365,365,365,366,as.numeric(as.Date("2021-7-19")-as.Date("2021-1-1")))

barplot(counts*365/times,width=times,col="royalblue",names=2013:2021,ylab="Assaults, scaled to /year",bg="white")
plot(2013:2021,counts*365/times,ylim=c(200,450),pch=19,col="royalblue",xlab="year",ylab="Assaults, scaled to /year",bg="white")
segments(2013:2021,(sqrt(counts)-1)^2*365/times, y1=(sqrt(counts)+1)^2*365/times,lty=3)

## from https://fyi.org.nz/request/16238/response/61990/attach/3/Kelly%20P%20IR%2001%2021%2023344%20Final%20Response.pdf

tslumley

pop<-matrix(0,102,102)

## uniform
vax<-matrix(rbinom(102*102,1,.83),102,102)

## vertical stripes
vax<-matrix(rep(c(0,1,1,1,1,1,1,0,1,1,1,1),length=102*102),102,102)

## high-risk square
#vax<-matrix(rbinom(102*102,1,.9),102,102); vax[1:41,1:42]<-rbinom(41*42,1,.5)

tslumley

expit<-function(eta) exp(eta)/(1+exp(eta))
logit<-function(p) log(p/(1-p))

one.sim<-function(beta, N=1000){
pop<-data.frame(x=rnorm(50000))
pop$y<-with(pop,rbinom(50000,1,expit(beta*x-3-beta)))
	
rval<-replicate(N,{	
cases<-which(pop$y==1)
controls<-sample(which(pop$y==0),length(cases))

tslumley

tslumley

Based on the `rimu` package, the start of an implementation of 'mr' objects using Hadley's `vctrs` package, so they can be used in tibbles. 

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)]

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

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

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

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];
     }