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

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

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

library(rvest)
library(dplyr)
library(DHBins)

a<-read_html("https://www.health.govt.nz/our-work/diseases-and-conditions/covid-19-novel-coronavirus/covid-19-data-and-statistics/covid-19-vaccine-data")
a %>% 
 html_node("body") %>% 
 xml2::xml_find_all("//table") %>% nth(n=3) %>%html_text() ->b

gsub("\t"""strsplit(b"\n")[[1]]fixed=TRUE)->d 

tslumley

library(dplyr)
library(rlang)
ma<-read.csv("~/Downloads/weekly-covid-19-municipality-vaccination-report-5-20-2021 (1)/Age - municipality-Table 1.csv",skip=1)
ma %>% 
	mutate(pop = as.numeric(gsub(",","",Population)), vax=as.numeric(gsub(",","",Fully.vaccinated.individuals))) -> ma_all
	
ma_all %>% 
    filter(Age.Group %in% c("65-74 Years","75+ Years")) %>%
	select(County,Town,Age.Group,pop,vax) %>% 
	group_by(County,Town) %>%

tslumley

svygofchisq<-function(formula, p, design,...){
	means<-svytotal(formula, design,...)
	rval<-chisq.test(means,p=p)
	nm<-names(coef(means))
	ncat<-length(coef(means))
	means<-svycontrast(means, list(N=rep(1,ncat)), add=TRUE)
	pN<-split(cbind(p,0*diag(p)), paste0("p_",nm))
	names(pN)<-paste0("p_",nm)
	means<-svycontrast(means, pN,add=TRUE)	
	for(i in 1:length(nm)){

tslumley

llprofci<-function(tpos,fpos,tneg,fneg){

    ndis<-tpos+fneg
    nnot<-tneg+fpos

    lljoint<-function(a,d){
	k<-length(a)

tslumley

---
title: "Rogue polls"
author: "Thomas Lumley"
date: "27/07/2020"
output: word_document
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```

tslumley

use zip-merged.dta
svyset SDMVPSU [pw=WTINT2YR], stra(SDMVSTRA)
svy: zip malepartners RIDAGEYR i.RIDRETH1 DMDEDUC, infl(RIDAGEYR i.RIDRETH1 DMDEDUC)

tslumley

nnot<-100
ndis<-50
tpr<-0.9
fpr<-0.5



rr<-replicate(100,{
tpos <-rbinom(1,ndis,tpr)
fpos<-rbinom(1,nnot,fpr)