😐

Thomas Lumley tslumley

😐

I do statistics and write software for it.

tslumley / truncate-platipus.R

Created September 25, 2025 05:46



	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 / eureka.R

Created November 14, 2024 00:37

"Pollen" data

tslumley / SE.R

Created November 3, 2024 00:27

Survey weighted SE

	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 / seven_of_nine.txt

Created November 12, 2023 05:54

Words of nine letters with seven distinct letters

tslumley / paired_test.R

Created September 1, 2023 09:52

	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.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{

tslumley / inprob.C

Created August 25, 2023 02:43

	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 / 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 / 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 / busdata.txt

Last active November 2, 2022 22:17

Buses