Thomas Lumley tslumley

## selfharm.csv

          
            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

## stabwt.R
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))

## vaxsim.R
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)

## police-assault.R
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

## DHB-vax.R
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

## ma.R
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) %>%

## svygofchisq.R
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)){

## proflike.R


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

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

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

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

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

## svgam-example.do
use zip-merged.dta
svyset SDMVPSU [pw=WTINT2YR], stra(SDMVSTRA)
svy: zip malepartners RIDAGEYR i.RIDRETH1 DMDEDUC, infl(RIDAGEYR i.RIDRETH1 DMDEDUC)
	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
	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))
	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(102102,1,.9),102,102); vax[1:41,1:42]<-rbinom(4142,1,.5)
	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)^2365/times, y1=(sqrt(counts)+1)^2365/times,lty=3)

	## from https://fyi.org.nz/request/16238/response/61990/attach/3/Kelly%20P%20IR%2001%2021%2023344%20Final%20Response.pdf
	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
	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) %>%
	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)){



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

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

	lljoint<-function(a,d){
	k<-length(a)
	---
	title: "Rogue polls"
	author: "Thomas Lumley"
	date: "27/07/2020"
	output: word_document
	---

	```{r setup, include=FALSE}
	knitr::opts_chunk$set(echo = TRUE)
	```
	use zip-merged.dta
	svyset SDMVPSU [pw=WTINT2YR], stra(SDMVSTRA)
	svy: zip malepartners RIDAGEYR i.RIDRETH1 DMDEDUC, infl(RIDAGEYR i.RIDRETH1 DMDEDUC)