Skip to content

Instantly share code, notes, and snippets.

View tslumley's full-sized avatar
😐

Thomas Lumley tslumley

😐
View GitHub Profile
@tslumley
tslumley / seven_of_nine.txt
Created November 12, 2023 05:54
Words of nine letters with seven distinct letters
abilities
academics
accessing
accessory
activated
addiction
additions
admission
advantage
aerospace
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
tslumley / new_sample_int.R
Last active August 31, 2023 10:28
Modify sample.int with new options for PPS with replacement sampling that have the specified inclusion probabilities
## 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{
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
tslumley / pairwise.R
Last active March 31, 2023 00:28
pairwise lmm with kinship
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
tslumley / brant_test.R
Created February 11, 2023 04:31
Based on Brant's test for the proportional odds assumption
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
tslumley / busdata.txt
Last active November 2, 2022 22:17
Buses
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
tslumley / pdl.R
Last active October 13, 2022 06:45
Polynomial distributed lag glms
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
tslumley / README
Created July 23, 2022 08:37
vctrs class for multiple-response objects
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
tslumley / selfharm.csv
Created March 20, 2022 05:20
Hospital admissions for intentional self-harm (from OIA)
Calendar year and month Total 00 to 04 05 to 09 10 to 14 15 to 19 20 to 24 25 to 29 30 to 34 35 to 39 40 to 44 45 to 49 50 to 54 55 to 59 60 to 64 65 to 69 70 to 74 75 to 79 80 to 84 85+ year month date
1 201701 801 NA 2 25 169 150 92 59 48 68 49 44 40 19 15 8 3 5 5 2017 1 2017-01-28
2 201702 828 NA NA 44 203 139 89 53 54 56 74 46 19 18 11 11 2 1 8 2017 2 2017-02-28
3 201703 880 NA NA 38 216 153 81 63 64 70 62 45 37 22 8 7 9 2 3 2017 3 2017-03-28
4 201704 790 NA 2 40 182 162 96 69 33 45 38 51 24 22 9 8 1 4 4 2017 4 2017-04-28
5 201705 915 NA 1 48 222 182 105 56 57 58 65 43 22 24 10 6 9 3 4 2017 5 2017-05-28
6 201706 828 NA NA 46 212 151 85 62 54 55 54 38 26 15 13 6 3 1 7 2017 6 2017-06-28
7 201707 830 NA NA 30 189 148 99 60 63 53 65 45 33 18 11 7 4 1 4 2017 7 2017-07-28
8 201708 905 NA 1 64 233 151 103 64 55 51 67 35 29 25 14 4 2 3 4 2017 8 2017-08-28
9 201709 868 NA NA 35 231 143 98 80 55 56 54 48 22 15 13 10 4 1 3 2017 9 2017-09-28