groceryheist/ml_measurement_error_public
1
0
Fork 0
Code Activity

Compare commits

base: groceryheist:irr
Branches Tags
groceryheist:main groceryheist:master groceryheist:git-annex groceryheist:synced/git-annex groceryheist:osf groceryheist:synced/osf groceryheist:synced/master groceryheist:irr
...
compare: groceryheist:synced/master
Branches Tags
groceryheist:master groceryheist:main groceryheist:git-annex groceryheist:synced/git-annex groceryheist:osf groceryheist:synced/osf groceryheist:synced/master groceryheist:irr

10 Commits
irr ... synced/mas

Author SHA1 Message Date
ntq8312
c42b94110b git-annex in ntq8312@kibo:/data/ntq8312/ml_measurement_error_public 2022-11-09 16:18:48 -06:00
Nathan TeBlunthuis
5c931a7198 git-annex in nathante@n3246:/gscratch/comdata/users/nathante/ml_measurement_error_public 2022-11-02 17:46:04 -07:00
Nathan TeBlunthuis
e17a52e236 add stuff to get perspective scores from civil comments 2022-11-02 17:45:35 -07:00
Nathan TeBlunthuis
d0c5766bdf Merge branch 'master' of code:ml_measurement_error_public 2022-10-07 10:43:32 -07:00
Nathan TeBlunthuis
b52b4f7daa Update stuff. 2022-10-07 10:42:50 -07:00
Nathan TeBlunthuis
8c846037bf add TADA stuff 2022-10-05 16:45:32 -07:00
Nathan TeBlunthuis
979dc14b68 update simulation and mle code 2022-09-02 11:34:50 -07:00
Nathan TeBlunthuis
b04f18afb7 refactor example charts 2022-09-01 17:52:25 -07:00
Nathan TeBlunthuis
47e9367ed5 Update the core 4 simulations. 2022-08-30 13:50:42 -07:00
Nathan TeBlunthuis
2cd447c327 simplify simulation 02. 2022-08-11 17:34:56 -07:00
48 changed files with 2361 additions and 656 deletions
Expand all files Collapse all files

3
.gitmodules vendored Normal file
View File

@@ -0,0 +1,3 @@
[submodule "paper"]
path = paper
url = git@github.com:chainsawriot/measure.git

9
charts/.gitignore vendored Normal file
View File

@@ -0,0 +1,9 @@
*.aux
*.svg
*.pdf
*.png
*.aux
*.fdb_latexmk
*.log
*.fls
auto/*

4
charts/example_3_dag/Makefile → charts/Makefile
View File

@@ -6,8 +6,7 @@ all: $(patsubst %.tex,%.svg,$(wildcard *.tex)) $(patsubst %.tex,%.png,$(wildcard
convert -density 300 -transparent white $< $@
%.svg: %.pdf
/usr/bin/inkscape $< --export-plain-svg=$@
/usr/bin/inkscape $< --export-plain-svg --export-type=svg --export-filename=$@
%.pdf: %.tex
latexmk -f -pdf $<
@@ -18,6 +17,7 @@ clean:
rm -f vc
rm *.svg
viewpdf: all
evince *.pdf

46
charts/example_1_dag.tex Normal file
View File

@@ -0,0 +1,46 @@
\documentclass[12pt]{standalone}
\usepackage{ucs}
\usepackage[utf8x]{inputenc}
\usepackage[T1]{fontenc}
\usepackage{textcomp}
\renewcommand{\rmdefault}{ugm}
\renewcommand{\sfdefault}{phv}
\usepackage[garamond]{mathdesign}
\usepackage{tikz}
\usetikzlibrary{positioning, shapes, arrows, shadows, }
\begin{document}
\tikzset{
observed/.style={circle, draw},
partly observed/.style 2 args={draw, fill=#2, path picture={
\fill[#1, sharp corners] (path picture bounding box.south west) -|
(path picture bounding box.north east) -- cycle;},
circle},
unobserved/.style={draw, circle, fill=gray!40},
residual/.style={draw, rectangle}
}
\tikzset{>=latex}
\begin{tikzpicture}
\node[observed] (y) {$Y$};
\node[unobserved, above=of y] (x) {$X$};
\node[observed, left=of x] (w) {$W$};
% \node[unobserved, above=of w] (k) {$K$};
\node[observed,right=of x] (z) {$Z$};
% \node[residual,below=of y] (e) {$\varepsilon$};
% \node[residual,below=of w] (xi) {$\xi$};
\draw[->] (z) to (y);
\draw[->] (z) -- (x);
\draw[->] (x) -- (y);
\draw[->] (x) -- (w);
% \draw[->] (y) -- (w);
% \draw[->] (x) -- (xi);
% \draw[->] (w) -- (xi);
\end{tikzpicture}
\end{document}

12
charts/example_1_dag/.gitignore vendored
View File

@@ -1,12 +0,0 @@
/auto/
/vc
/refs-processed.bib
/*.aux
/*.bbl
/*.blg
/*.log
/*.out
/*.pdf
/*.fdb_latexmk
/*.fls
/*.ttt

28
charts/example_1_dag/Makefile
View File

@@ -1,28 +0,0 @@
#!/usr/bin/make
all: $(patsubst %.tex,%.svg,$(wildcard *.tex)) $(patsubst %.tex,%.png,$(wildcard *.tex))
%.png: %.pdf
convert -density 300 -transparent white $< $@
%.svg: %.pdf
/usr/bin/inkscape $< --export-plain-svg=$@
%.pdf: %.tex
latexmk -f -pdf $<
clean:
latexmk -C *.tex
rm -f *.tmp
rm -f vc
rm *.svg
viewpdf: all
evince *.pdf
vc:
vc-git
pdf: all
.PHONY: clean all

46
charts/example_1_dag/example_1_dag.tex
View File

@@ -1,46 +0,0 @@
\documentclass[12pt]{standalone}
\usepackage{ucs}
\usepackage[utf8x]{inputenc}
\usepackage[T1]{fontenc}
\usepackage{textcomp}
\renewcommand{\rmdefault}{ugm}
\renewcommand{\sfdefault}{phv}
\usepackage[garamond]{mathdesign}
\usepackage{tikz}
\usetikzlibrary{positioning, shapes, arrows, shadows, }
\begin{document}
\tikzset{>=latex}
\begin{tikzpicture}[
observed/.style={circle, draw},
unobserved/.style={draw, circle, fill=gray!40},
legend/.style={rectangle, draw},
partly observed/.style 2 args={draw, fill=#2, path picture={
\fill[#1, sharp corners] (path picture bounding box.south west) -|
(path picture bounding box.north east) -- cycle;},
circle}
]
\node[observed] (y) {Y};
% \node[observed,above=of y, xshift=1cm] (r) {R};
\node[partly observed={white}{gray!40}, above=of y] (x) {X};
\node[observed, left=of x,xshift=-1] (w) {W};
\node[unobserved, above=of w, xshift=1cm] (k) {K};
\matrix [draw, below, yshift=-0.2cm, font=\small, align=center, column sep=2\pgflinewidth, inner sep=0.6em, outer sep=0em, nodes={align=center, anchor=center}] at (current bounding box.south){
\node[observed,label=right:observed] {}; \\
\node[unobserved,label=right:unobserved]{}; \\
};
\draw[->] (x) -- (y);
\draw[->] (x) -- (w);
\draw[-] (k) -- (x);
\draw[->] (k) -- (w);
\end{tikzpicture}
\end{document}

47
charts/example_2_dag.tex Normal file
View File

@@ -0,0 +1,47 @@
\documentclass[12pt]{standalone}
\usepackage{ucs}
\usepackage[utf8x]{inputenc}
\usepackage[T1]{fontenc}
\usepackage{textcomp}
\renewcommand{\rmdefault}{ugm}
\renewcommand{\sfdefault}{phv}
\usepackage[garamond]{mathdesign}
\usepackage{tikz}
\usetikzlibrary{positioning, shapes, arrows, shadows}
\begin{document}
\tikzset{
observed/.style={circle, draw},
partly observed/.style 2 args={draw, fill=#2, path picture={
\fill[#1, sharp corners] (path picture bounding box.south west) -|
(path picture bounding box.north east) -- cycle;},
circle},
unobserved/.style={draw, circle, fill=gray!40},
residual/.style={draw, rectangle}
}
\tikzset{>=latex}
\begin{tikzpicture}
\node[observed] (y) {$Y$};
\node[unobserved, above=of y] (x) {$X$};
\node[observed, left=of x] (w) {$W$};
\node[observed,right=of x] (z) {$Z$};
\draw[->] (z) to (y);
\draw[->] (z) -- (x);
\draw[->] (x) -- (y);
\draw[->] (x) -- (w);
\draw[->] (x) to (y);
\draw[->] (y) -- (w);
\draw[->] (y) -- (w);
\end{tikzpicture}
\end{document}

12
charts/example_2_dag/.gitignore vendored
View File

@@ -1,12 +0,0 @@
/auto/
/vc
/refs-processed.bib
/*.aux
/*.bbl
/*.blg
/*.log
/*.out
/*.pdf
/*.fdb_latexmk
/*.fls
/*.ttt

28
charts/example_2_dag/Makefile
View File

@@ -1,28 +0,0 @@
#!/usr/bin/make
all: $(patsubst %.tex,%.svg,$(wildcard *.tex)) $(patsubst %.tex,%.png,$(wildcard *.tex))
%.png: %.pdf
convert -density 300 -transparent white $< $@
%.svg: %.pdf
/usr/bin/inkscape $< --export-plain-svg=$@
%.pdf: %.tex
latexmk -f -pdf $<
clean:
latexmk -C *.tex
rm -f *.tmp
rm -f vc
rm *.svg
viewpdf: all
evince *.pdf
vc:
vc-git
pdf: all
.PHONY: clean all

47
charts/example_2_dag/example_2_dag.tex
View File

@@ -1,47 +0,0 @@
\documentclass[12pt]{standalone}
\usepackage{ucs}
\usepackage[utf8x]{inputenc}
\usepackage[T1]{fontenc}
\usepackage{textcomp}
\renewcommand{\rmdefault}{ugm}
\renewcommand{\sfdefault}{phv}
\usepackage[garamond]{mathdesign}
\usepackage{tikz}
\usetikzlibrary{positioning, shapes, arrows, shadows}
\begin{document}
\tikzset{>=latex}
\begin{tikzpicture}[
observed/.style={circle, draw},
partly observed/.style 2 args={draw, fill=#2, path picture={
\fill[#1, sharp corners] (path picture bounding box.south west) -|
(path picture bounding box.north east) -- cycle;},
circle},
unobserved/.style={draw, circle, fill=gray!40}
]
\node[observed] (y) {Y};
\node[observed,above=of y, xshift=1cm] (r) {R};
\node[partly observed={white}{gray!40}, left=of r] (x) {X};
\node[observed, left=of x] (w) {W};
\node[unobserved, above=of w, xshift=1cm] (k) {K};
\matrix [draw, below, yshift=-0.2cm, font=\small, align=center, column sep=2\pgflinewidth, inner sep=0.6em, outer sep=0em, nodes={align=center, anchor=center}] at (current bounding box.south){
\node[observed,label=right:observed] {}; \\
\node[unobserved,label=right:unobserved]{}; \\
};
\draw[->] (r) -- (y);
\draw[->] (r) -- (x);
\draw[->] (x) -- (y);
\draw[->] (x) -- (w);
\draw[->] (k) -- (x);
\draw[->] (k) -- (w);
\end{tikzpicture}
\end{document}

42
charts/example_3_dag.tex Normal file
View File

@@ -0,0 +1,42 @@
\documentclass[12pt]{standalone}
\usepackage{ucs}
\usepackage[utf8x]{inputenc}
\usepackage[T1]{fontenc}
\usepackage{textcomp}
\renewcommand{\rmdefault}{ugm}
\renewcommand{\sfdefault}{phv}
\usepackage[garamond]{mathdesign}
\usepackage{tikz}
\usetikzlibrary{positioning, shapes, arrows, shadows}
\begin{document}
\tikzset{
observed/.style={circle, draw},
partly observed/.style 2 args={draw, fill=#2, path picture={
\fill[#1, sharp corners] (path picture bounding box.south west) -|
(path picture bounding box.north east) -- cycle;},
circle},
unobserved/.style={draw, circle, fill=gray!40},
residual/.style={draw, rectangle}
}
\tikzset{>=latex}
\begin{tikzpicture}
\node[unobserved] (y) {$Y$};
\node[observed, above=of y] (x) {$X$};
\node[observed, left=of x] (w) {$W$};
% \node[unobserved, above=of w] (k) {$K$};
\node[observed,right=of x] (z) {$Z$};
% \node[residual,below=of y] (e) {$\varepsilon$};
% \node[residual,below=of w] (xi) {$\xi$};
\draw[->] (z) to (y);
\draw[->] (x) -- (y);
\draw[->] (y) -- (w);
\end{tikzpicture}
\end{document}

12
charts/example_3_dag/.gitignore vendored
View File

@@ -1,12 +0,0 @@
/auto/
/vc
/refs-processed.bib
/*.aux
/*.bbl
/*.blg
/*.log
/*.out
/*.pdf
/*.fdb_latexmk
/*.fls
/*.ttt

51
charts/example_3_dag/example_3_dag.tex
View File

@@ -1,51 +0,0 @@
\documentclass[12pt]{standalone}
\usepackage{ucs}
\usepackage[utf8x]{inputenc}
\usepackage[T1]{fontenc}
\usepackage{textcomp}
\renewcommand{\rmdefault}{ugm}
\renewcommand{\sfdefault}{phv}
\usepackage[garamond]{mathdesign}
\usepackage{tikz}
\usetikzlibrary{positioning, shapes, arrows, shadows}
\begin{document}
\tikzset{>=latex}
\begin{tikzpicture}[
observed/.style={circle, draw},
partly observed/.style 2 args={draw, fill=#2, path picture={
\fill[#1, sharp corners] (path picture bounding box.south west) -|
(path picture bounding box.north east) -- cycle;},
circle},
unobserved/.style={draw, circle, fill=gray!40}
]
\node[observed] (y) {Y};
% \node[observed,above=of y, xshift=1cm] (r) {R};
\node[partly observed={white}{gray!40}, above = of y] (x) {X};
\node[observed, left=of x] (w) {W};
\node[unobserved, above=of w, xshift=1cm] (k) {K};
\node[unobserved, right=of k] (u) {U};
% \draw[->] (r) -- (y);
% \draw[->] (r) -- (x);
\draw[->] (x) -- (y);
\draw[->] (x) -- (w);
\draw[->] (x) -- (k);
\draw[->] (k) -- (w);
\draw[->] (u) to [out=270,in=30] (y);
\draw[->] (u) -- (k);
\matrix [draw, below, yshift=-0.2cm, inner sep=0.6em, outer sep=0em, nodes={align=center, anchor=center}] at (current bounding box.south){
\node[observed,label=right:observed] {}; \\
\node[unobserved,label=right:unobserved]{}; \\
};
\end{tikzpicture}
\end{document}

47
charts/example_4_dag.tex Normal file
View File

@@ -0,0 +1,47 @@
\documentclass[12pt]{standalone}
\usepackage{ucs}
\usepackage[utf8x]{inputenc}
\usepackage[T1]{fontenc}
\usepackage{textcomp}
\renewcommand{\rmdefault}{ugm}
\renewcommand{\sfdefault}{phv}
\usepackage[garamond]{mathdesign}
\usepackage{tikz}
\usetikzlibrary{positioning, shapes, arrows, shadows}
\begin{document}
\tikzset{
observed/.style={circle, draw},
partly observed/.style 2 args={draw, fill=#2, path picture={
\fill[#1, sharp corners] (path picture bounding box.south west) -|
(path picture bounding box.north east) -- cycle;},
circle},
unobserved/.style={draw, circle, fill=gray!40},
residual/.style={draw, rectangle}
}
\tikzset{>=latex}
\begin{tikzpicture}
\node[unobserved] (y) {$Y$};
\node[observed={white}{gray!40}, above=of y] (x) {$X$};
\node[observed, left=of x] (w) {$W$};
% \node[unobserved, above=of w] (k) {$K$};
\node[observed,right=of x] (z) {$Z$};
% \node[residual,below=of y] (e) {$\varepsilon$};
% \node[residual,below=of w] (xi) {$\xi$};
\draw[->] (x) -- (y);
% \draw[->] (x) -- (w);
\draw[->] (y) -- (w);
% \draw[->] (k) -- (w);
\draw[->] (z) -- (y);
% \draw[->] (z) -- (k);
% \draw[->] (y) -- (xi);
% \draw[->] (w) -- (xi);
\end{tikzpicture}
\end{document}

35
charts/example_legend.tex Normal file
View File

@@ -0,0 +1,35 @@
\documentclass[12pt]{standalone}
\usepackage{ucs}
\usepackage[utf8x]{inputenc}
\usepackage[T1]{fontenc}
\usepackage{textcomp}
\renewcommand{\rmdefault}{ugm}
\renewcommand{\sfdefault}{phv}
\usepackage[garamond]{mathdesign}
\usepackage{tikz}
\usetikzlibrary{positioning, shapes, arrows, shadows}
\begin{document}
\tikzset{
observed/.style={circle, draw},
partly observed/.style 2 args={draw, fill=#2, path picture={
\fill[#1, sharp corners] (path picture bounding box.south west) -|
(path picture bounding box.north east) -- cycle;},
circle},
unobserved/.style={draw, circle, fill=gray!40},
residual/.style={draw, rectangle}
}
\tikzset{>=latex}
\begin{tikzpicture}
\matrix [draw, below, font=\small, align=center, column sep=2\pgflinewidth, inner sep=0.4em, outer sep=0em, nodes={align=center, anchor=center}] at (current bounding box.south){
\node[observed,label=right:observed] {}; \\
\node[unobserved,label=right:automatically classified]{}; \\
\node[residual,label=right:error term]{}; \\
};
\end{tikzpicture}
\end{document}

1
civil_comments/all_data.csv Symbolic link
View File

@@ -0,0 +1 @@
../.git/annex/objects/6v/fJ/SHA256E-s916052376--a85b5ba7e9a8cda38b91ea6e3957a4f2bfff17bb52f22c935595cbe47cc54d94.csv/SHA256E-s916052376--a85b5ba7e9a8cda38b91ea6e3957a4f2bfff17bb52f22c935595cbe47cc54d94.csv
1 ../.git/annex/objects/6v/fJ/SHA256E-s916052376--a85b5ba7e9a8cda38b91ea6e3957a4f2bfff17bb52f22c935595cbe47cc54d94.csv/SHA256E-s916052376--a85b5ba7e9a8cda38b91ea6e3957a4f2bfff17bb52f22c935595cbe47cc54d94.csv

39
civil_comments/get_perspective_scores.py Normal file
View File

@@ -0,0 +1,39 @@
from googleapiclient import discovery
import json
import csv
from pathlib import Path
from time import sleep
from itertools import islice
API_KEY = open('perspective_api_key').read()
client = discovery.build("commentanalyzer","v1alpha",developerKey=API_KEY,discoveryServiceUrl="https://commentanalyzer.googleapis.com/$discovery/rest?version=v1alpha1",static_discovery=False,)
csvreader = csv.DictReader(open("all_data.csv"), dialect='unix')
outfile = Path("perspective_results.json")
already_scored = set()
if outfile.exists():
already_scored = set([json.loads(l)['id'] for l in open(str(outfile),'r')])
print(f"loaded {len(already_scored)} scored comments")
with open("perspective_results.json",'a') as of:
for line in csvreader:
if line['id'] not in already_scored and len(line.get('comment_text','')) > 0:
analyze_request = {'comment':{'text':line['comment_text']},
'languages':['en'],
'requestedAttributes':{'TOXICITY':{},
"SEVERE_TOXICITY":{},
"IDENTITY_ATTACK":{},
"INSULT":{},
"PROFANITY":{},
"THREAT":{}}}
response = client.comments().analyze(body=analyze_request).execute()
response['id'] = line['id']
result = json.dumps(response)
of.write(result + '\n')
of.flush()
sleep(0.10)

1
civil_comments/identity_individual_annotations.csv Symbolic link
View File

@@ -0,0 +1 @@
../.git/annex/objects/qP/Xw/SHA256E-s106388260--7b8e9f21c5110d32e337137f8b4fe50987ec1b59fdbfd56a4717cdc13e509ec3.csv/SHA256E-s106388260--7b8e9f21c5110d32e337137f8b4fe50987ec1b59fdbfd56a4717cdc13e509ec3.csv
1 ../.git/annex/objects/qP/Xw/SHA256E-s106388260--7b8e9f21c5110d32e337137f8b4fe50987ec1b59fdbfd56a4717cdc13e509ec3.csv/SHA256E-s106388260--7b8e9f21c5110d32e337137f8b4fe50987ec1b59fdbfd56a4717cdc13e509ec3.csv

BIN
civil_comments/perspective_api_key.gpg Normal file
View File

Binary file not shown.

1
civil_comments/perspective_results.json Symbolic link
View File

@@ -0,0 +1 @@
../.git/annex/objects/ZV/z8/SHA256E-s2293825121--6cdc8f8fb64fad2e51027e2564928e8938bf5fc6ca0cd6c31cb2e67aafe0a203.json/SHA256E-s2293825121--6cdc8f8fb64fad2e51027e2564928e8938bf5fc6ca0cd6c31cb2e67aafe0a203.json

1
civil_comments/toxicity_individual_annotations.csv Symbolic link
View File

@@ -0,0 +1 @@
../.git/annex/objects/FF/WZ/SHA256E-s417648663--c85bda15b964a24869ae11f76092bde6f4b18236dd1cbe17539526b3b5b736cf.csv/SHA256E-s417648663--c85bda15b964a24869ae11f76092bde6f4b18236dd1cbe17539526b3b5b736cf.csv
1 ../.git/annex/objects/FF/WZ/SHA256E-s417648663--c85bda15b964a24869ae11f76092bde6f4b18236dd1cbe17539526b3b5b736cf.csv/SHA256E-s417648663--c85bda15b964a24869ae11f76092bde6f4b18236dd1cbe17539526b3b5b736cf.csv

1
paper Submodule

Submodule paper added at b135cac19e

BIN
presentations/TADA_2022/Tada Slides.pptx.pdf Normal file
View File

Binary file not shown.

BIN
presentations/TADA_2022/tada_2020.mkv Normal file
View File

Binary file not shown.

BIN
presentations/TADA_2022/tada_2020.mp4 Normal file
View File

Binary file not shown.

91
simulations/.Rhistory
View File

@@ -63,3 +63,94 @@ list.files()
install.packages("filelock")
q()
n
df
df
outcome_formula <- y ~ x + z
outcome_family=gaussian()
proxy_formula <- w_pred ~ x
truth_formula <- x ~ z
params <- start
ll.y.obs.x0
ll.y.obs.x1
rater_formula <- x.obs ~ x
rater_formula
rater.modle.matrix.obs.x0
rater.model.matrix.obs.x0
names(rater.model.matrix.obs.x0)
head(rater.model.matrix.obs.x0)
df.obs
ll.x.obs.0
rater.params
rater.params %*% t(rater.model.matrix.x.obs.0[df.obs$xobs.0==1])
df.obs$xobs.0==1
df.obs$x.obs.0==1
ll.x.obs.0[df.obs$x.obs.0==1]
rater.model.matrix.x.obs.0[df.obs$x.obs.0==1,]
df.obs$x.obs.0==1
n.rater.model.covars <- dim(rater.model.matrix.x.obs.0)[2]
rater.params <- params[param.idx:n.rater.model.covars]
rater.params
ll.x.obs.0[df.obs$x.obs.0==1] <- plogis(rater.params %*% t(rater.model.matrix.x.obs.0[df.obs$x.obs.0==1,]), log=TRUE)
t(rater.model.matrix.x.obs.0[df.obs$x.obs.0==1,]
)
dimt(rater.model.matrix.x.obs.0[df.obs$x.obs.0==1,])
dim(t(rater.model.matrix.x.obs.0[df.obs$x.obs.0==1,]))
dim(ll.x.obs.0[df.obs$x.obs.0==1])
rater.params
rater.params
rater.params
rater_formula
rater.params
)
1+1
q()
n
outcome_formula <- y ~ x + z
proxy_formula <- w_pred ~ x + z + y
truth_formula <- x ~ z
proxy_formula
eyboardio Model 01 - Kaleidoscope locally built
df <- df.triple.proxy.mle
outcome_family='gaussian'
outcome_family=gaussian()
proxy_formulas=list(proxy_formula,x.obs.0~x, x.obs.1~x)
proxy_formulas
proxy_familites <- rep(binomial(link='logit'),3)
proxy_families = rep(binomial(link='logit'),3)
proxy_families
proxy_families = list(binomial(link='logit'),binomial(link='logit'),binomial(link='logit'))
proxy_families
proxy_families[[1]]
proxy.params
i
proxy_params
proxy.params
params
params <- start
df.triple.proxy.mle
df
coder.formulas <- c(x.obs.0 ~ x, x.obs.1 ~x)
outcome.formula
outcome_formula
depvar(outcome_formula
)
outcome_formula$terms
terms(outcome_formula)
q()
n
df.triple.proxy.mle
triple.proxy.mle
df
df <- df.triple.proxy
outcome_family <- binomial(link='logit')
outcome_formula <- y ~x+z
proxy_formula <- w_pred ~ y
coder_formulas=list(y.obs.1~y,y.obs.2~y); proxy_formula=w_pred~y; proxy_family=binomial(link='logit'))
coder_formulas=list(y.obs.1~y,y.obs.2~y); proxy_formula=w_pred~y; proxy_family=binomial(link='logit')
coder_formulas=list(y.obs.0~y,y.obs.1~y)
traceback()
df
df
outcome.model.matrix
q()
n

46
simulations/01_two_covariates.R
View File

@@ -32,7 +32,7 @@ source("simulation_base.R")
simulate_data <- function(N, m, B0=0, Bxy=0.2, Bzy=-0.2, Bzx=0.2, y_explained_variance=0.025, prediction_accuracy=0.73, seed=1){
set.seed(seed)
z <- rbinom(N, 1, 0.5)
z <- rnorm(N,sd=0.5)
# x.var.epsilon <- var(Bzx *z) * ((1-zx_explained_variance)/zx_explained_variance)
xprime <- Bzx * z #+ x.var.epsilon
x <- rbinom(N,1,plogis(xprime))
@@ -71,31 +71,35 @@ parser <- add_argument(parser, "--outfile", help='output file', default='example
parser <- add_argument(parser, "--y_explained_variance", help='what proportion of the variance of y can be explained?', default=0.05)
# parser <- add_argument(parser, "--zx_explained_variance", help='what proportion of the variance of x can be explained by z?', default=0.3)
parser <- add_argument(parser, "--prediction_accuracy", help='how accurate is the predictive model?', default=0.73)
parser <- add_argument(parser, "--Bzx", help='coefficient of z on x?', default=1)
args <- parse_args(parser)
parser <- add_argument(parser, "--outcome_formula", help='formula for the outcome variable', default="y~x+z")
parser <- add_argument(parser, "--proxy_formula", help='formula for the proxy variable', default="w_pred~x")
parser <- add_argument(parser, "--truth_formula", help='formula for the true variable', default="x~z")
parser <- add_argument(parser, "--Bzx", help='Effect of z on x', default=0.3)
parser <- add_argument(parser, "--Bzy", help='Effect of z on y', default=-0.3)
parser <- add_argument(parser, "--Bxy", help='Effect of x on y', default=0.3)
args <- parse_args(parser)
B0 <- 0
Bxy <- 0.3
Bzy <- -0.3
Bxy <- args$Bxy
Bzy <- args$Bzy
Bzx <- args$Bzx
if (args$m < args$N){
df <- simulate_data(args$N, args$m, B0, Bxy, Bzy, Bzx, seed=args$seed + 500, y_explained_variance = args$y_explained_variance, prediction_accuracy=args$prediction_accuracy)
df <- simulate_data(args$N, args$m, B0, Bxy, Bzy, Bzx, seed=args$seed + 500, y_explained_variance = args$y_explained_variance, prediction_accuracy=args$prediction_accuracy)
result <- list('N'=args$N,'m'=args$m,'B0'=B0,'Bxy'=Bxy, Bzx=Bzx, 'Bzy'=Bzy, 'seed'=args$seed, 'y_explained_variance'=args$y_explained_variance, 'prediction_accuracy'=args$prediction_accuracy, 'accuracy_imbalance_difference'=args$accuracy_imbalance_difference, 'outcome_formula'=args$outcome_formula, 'proxy_formula'=args$proxy_formula,truth_formula=args$truth_formula, error='')
result <- list('N'=args$N,'m'=args$m,'B0'=B0,'Bxy'=Bxy,'Bzy'=Bzy, 'Bzx'=Bzx, 'seed'=args$seed, 'y_explained_variance' = args$y_explained_variance, 'zx_explained_variance' = args$zx_explained_variance, "prediction_accuracy"=args$prediction_accuracy, "error"="")
outline <- run_simulation(df, result, outcome_formula=as.formula(args$outcome_formula), proxy_formula=as.formula(args$proxy_formula), truth_formula=as.formula(args$truth_formula))
outline <- run_simulation(df, result)
outfile_lock <- lock(paste0(args$outfile, '_lock'),exclusive=TRUE)
if(file.exists(args$outfile)){
logdata <- read_feather(args$outfile)
logdata <- rbind(logdata,as.data.table(outline),fill=TRUE)
} else {
logdata <- as.data.table(outline)
}
print(outline)
write_feather(logdata, args$outfile)
unlock(outfile_lock)
outfile_lock <- lock(paste0(args$outfile, '_lock'),exclusive=TRUE)
if(file.exists(args$outfile)){
logdata <- read_feather(args$outfile)
logdata <- rbind(logdata,as.data.table(outline),fill=TRUE)
} else {
logdata <- as.data.table(outline)
}
print(outline)
write_feather(logdata, args$outfile)
unlock(outfile_lock)

146
simulations/02_indep_differential.R
View File

@@ -31,11 +31,11 @@ source("simulation_base.R")
## one way to do it is by adding correlation to x.obs and y that isn't in w.
## in other words, the model is missing an important feature of x.obs that's related to y.
simulate_data <- function(N, m, B0, Bxy, Bzx, Bzy, seed, y_explained_variance=0.025, prediction_accuracy=0.73, accuracy_imbalance_difference=0.3){
simulate_data <- function(N, m, B0, Bxy, Bzx, Bzy, seed, y_explained_variance=0.025, prediction_accuracy=0.73, y_bias=-0.8,z_bias=0,accuracy_imbalance_difference=0.3){
set.seed(seed)
# make w and y dependent
z <- rbinom(N, 1, 0.5)
x <- rbinom(N, 1, Bzx * z + 0.5)
z <- rnorm(N,sd=0.5)
x <- rbinom(N, 1, plogis(Bzx * z))
y.var.epsilon <- (var(Bzy * z) + var(Bxy *x) + 2*cov(Bzy*z,Bxy*x)) * ((1-y_explained_variance)/y_explained_variance)
y.epsilon <- rnorm(N, sd = sqrt(y.var.epsilon))
@@ -49,109 +49,119 @@ simulate_data <- function(N, m, B0, Bxy, Bzx, Bzy, seed, y_explained_variance=0.
df <- df[, x.obs := x]
}
## px <- mean(x)
## probablity of an error is correlated with y
## pz <- mean(z)
## accuracy_imbalance_ratio <- (prediction_accuracy + accuracy_imbalance_difference/2) / (prediction_accuracy - accuracy_imbalance_difference/2)
## # this works because of conditional probability
## accuracy_x0 <- prediction_accuracy / (px*(accuracy_imbalance_ratio) + (1-px))
## accuracy_x1 <- accuracy_imbalance_ratio * accuracy_x0
## accuracy_z0 <- prediction_accuracy / (pz*(accuracy_imbalance_ratio) + (1-pz))
## accuracy_z1 <- accuracy_imbalance_ratio * accuracy_z0
## x0 <- df[x==0]$x
## x1 <- df[x==1]$x
## nx1 <- nrow(df[x==1])
## nx0 <- nrow(df[x==0])
## z0x0 <- df[(z==0) & (x==0)]$x
## z0x1 <- df[(z==0) & (x==1)]$x
## z1x0 <- df[(z==1) & (x==0)]$x
## z1x1 <- df[(z==1) & (x==1)]$x
## yx0 <- df[x==0]$y
## yx1 <- df[x==1]$y
## yz0x0 <- df[(z==0) & (x==0)]$y
## yz0x1 <- df[(z==0) & (x==1)]$y
## yz1x0 <- df[(z==1) & (x==0)]$y
## yz1x1 <- df[(z==1) & (x==1)]$y
# tranform yz0.1 into a logistic distribution with mean accuracy_z0
## acc.x0 <- plogis(0.5*scale(yx0) + qlogis(accuracy_x0))
## acc.x1 <- plogis(1.5*scale(yx1) + qlogis(accuracy_x1))
## nz0x0 <- nrow(df[(z==0) & (x==0)])
## nz0x1 <- nrow(df[(z==0) & (x==1)])
## nz1x0 <- nrow(df[(z==1) & (x==0)])
## nz1x1 <- nrow(df[(z==1) & (x==1)])
## w0x0 <- (1-x0)**2 + (-1)**(1-x0) * acc.x0
## w0x1 <- (1-x1)**2 + (-1)**(1-x1) * acc.x1
pz <- mean(z)
accuracy_imbalance_ratio <- (prediction_accuracy + accuracy_imbalance_difference/2) / (prediction_accuracy - accuracy_imbalance_difference/2)
## yz1 <- df[z==1]$y
## yz1 <- df[z==1]$y
# this works because of conditional probability
accuracy_z0 <- prediction_accuracy / (pz*(accuracy_imbalance_ratio) + (1-pz))
accuracy_z1 <- accuracy_imbalance_ratio * accuracy_z0
## # tranform yz0.1 into a logistic distribution with mean accuracy_z0
## acc.z0x0 <- plogis(0.5*scale(yz0x0) + qlogis(accuracy_z0))
## acc.z0x1 <- plogis(0.5*scale(yz0x1) + qlogis(accuracy_z0))
## acc.z1x0 <- plogis(1.5*scale(yz1x0) + qlogis(accuracy_z1))
## acc.z1x1 <- plogis(1.5*scale(yz1x1) + qlogis(accuracy_z1))
z0x0 <- df[(z==0) & (x==0)]$x
z0x1 <- df[(z==0) & (x==1)]$x
z1x0 <- df[(z==1) & (x==0)]$x
z1x1 <- df[(z==1) & (x==1)]$x
## w0z0x0 <- (1-z0x0)**2 + (-1)**(1-z0x0) * acc.z0x0
## w0z0x1 <- (1-z0x1)**2 + (-1)**(1-z0x1) * acc.z0x1
## w0z1x0 <- (1-z1x0)**2 + (-1)**(1-z1x0) * acc.z1x0
## w0z1x1 <- (1-z1x1)**2 + (-1)**(1-z1x1) * acc.z1x1
yz0x0 <- df[(z==0) & (x==0)]$y
yz0x1 <- df[(z==0) & (x==1)]$y
yz1x0 <- df[(z==1) & (x==0)]$y
yz1x1 <- df[(z==1) & (x==1)]$y
## ##perrorz0 <- w0z0*(pyz0)
## ##perrorz1 <- w0z1*(pyz1)
nz0x0 <- nrow(df[(z==0) & (x==0)])
nz0x1 <- nrow(df[(z==0) & (x==1)])
nz1x0 <- nrow(df[(z==1) & (x==0)])
nz1x1 <- nrow(df[(z==1) & (x==1)])
## w0z0x0.noisy.odds <- rlogis(nz0x0,qlogis(w0z0x0))
## w0z0x1.noisy.odds <- rlogis(nz0x1,qlogis(w0z0x1))
## w0z1x0.noisy.odds <- rlogis(nz1x0,qlogis(w0z1x0))
## w0z1x1.noisy.odds <- rlogis(nz1x1,qlogis(w0z1x1))
yz1 <- df[z==1]$y
yz1 <- df[z==1]$y
## df[(z==0)&(x==0),w:=plogis(w0z0x0.noisy.odds)]
## df[(z==0)&(x==1),w:=plogis(w0z0x1.noisy.odds)]
## df[(z==1)&(x==0),w:=plogis(w0z1x0.noisy.odds)]
## df[(z==1)&(x==1),w:=plogis(w0z1x1.noisy.odds)]
# tranform yz0.1 into a logistic distribution with mean accuracy_z0
acc.z0x0 <- plogis(0.5*scale(yz0x0) + qlogis(accuracy_z0))
acc.z0x1 <- plogis(0.5*scale(yz0x1) + qlogis(accuracy_z0))
acc.z1x0 <- plogis(1.5*scale(yz1x0) + qlogis(accuracy_z1))
acc.z1x1 <- plogis(1.5*scale(yz1x1) + qlogis(accuracy_z1))
## df[,w_pred:=as.integer(w > 0.5)]
## print(mean(df[z==0]$x == df[z==0]$w_pred))
## print(mean(df[z==1]$x == df[z==1]$w_pred))
## print(mean(df$w_pred == df$x))
w0z0x0 <- (1-z0x0)**2 + (-1)**(1-z0x0) * acc.z0x0
w0z0x1 <- (1-z0x1)**2 + (-1)**(1-z0x1) * acc.z0x1
w0z1x0 <- (1-z1x0)**2 + (-1)**(1-z1x0) * acc.z1x0
w0z1x1 <- (1-z1x1)**2 + (-1)**(1-z1x1) * acc.z1x1
resids <- resid(lm(y~x + z))
odds.x1 <- qlogis(prediction_accuracy) + y_bias*qlogis(pnorm(resids[x==1])) + z_bias * qlogis(pnorm(z,sd(z)))
odds.x0 <- qlogis(prediction_accuracy,lower.tail=F) + y_bias*qlogis(pnorm(resids[x==0])) + z_bias * qlogis(pnorm(z,sd(z)))
##perrorz0 <- w0z0*(pyz0)
##perrorz1 <- w0z1*(pyz1)
## acc.x0 <- p.correct[df[,x==0]]
## acc.x1 <- p.correct[df[,x==1]]
w0z0x0.noisy.odds <- rlogis(nz0x0,qlogis(w0z0x0))
w0z0x1.noisy.odds <- rlogis(nz0x1,qlogis(w0z0x1))
w0z1x0.noisy.odds <- rlogis(nz1x0,qlogis(w0z1x0))
w0z1x1.noisy.odds <- rlogis(nz1x1,qlogis(w0z1x1))
df[x==0,w:=plogis(rlogis(.N,odds.x0))]
df[x==1,w:=plogis(rlogis(.N,odds.x1))]
df[,w_pred := as.integer(w > 0.5)]
df[(z==0)&(x==0),w:=plogis(w0z0x0.noisy.odds)]
df[(z==0)&(x==1),w:=plogis(w0z0x1.noisy.odds)]
df[(z==1)&(x==0),w:=plogis(w0z1x0.noisy.odds)]
df[(z==1)&(x==1),w:=plogis(w0z1x1.noisy.odds)]
df[,w_pred:=as.integer(w > 0.5)]
print(mean(df[z==0]$x == df[z==0]$w_pred))
print(mean(df[z==1]$x == df[z==1]$w_pred))
print(mean(df$w_pred == df$x))
print(mean(df[y>=0]$w_pred == df[y>=0]$x))
print(mean(df[y<=0]$w_pred == df[y<=0]$x))
return(df)
}
parser <- arg_parser("Simulate data and fit corrected models")
parser <- add_argument(parser, "--N", default=1400, help="number of observations of w")
parser <- add_argument(parser, "--N", default=5000, help="number of observations of w")
parser <- add_argument(parser, "--m", default=500, help="m the number of ground truth observations")
parser <- add_argument(parser, "--seed", default=50, help='seed for the rng')
parser <- add_argument(parser, "--seed", default=51, help='seed for the rng')
parser <- add_argument(parser, "--outfile", help='output file', default='example_2.feather')
parser <- add_argument(parser, "--y_explained_variance", help='what proportion of the variance of y can be explained?', default=0.01)
parser <- add_argument(parser, "--prediction_accuracy", help='how accurate is the predictive model?', default=0.73)
parser <- add_argument(parser, "--y_explained_variance", help='what proportion of the variance of y can be explained?', default=0.1)
parser <- add_argument(parser, "--prediction_accuracy", help='how accurate is the predictive model?', default=0.75)
parser <- add_argument(parser, "--accuracy_imbalance_difference", help='how much more accurate is the predictive model for one class than the other?', default=0.3)
parser <- add_argument(parser, "--Bzx", help='Effect of z on x', default=0.3)
parser <- add_argument(parser, "--Bzy", help='Effect of z on y', default=-0.3)
parser <- add_argument(parser, "--Bxy", help='Effect of z on y', default=0.3)
parser <- add_argument(parser, "--outcome_formula", help='formula for the outcome variable', default="y~x+z")
parser <- add_argument(parser, "--proxy_formula", help='formula for the proxy variable', default="w_pred~y*z*x")
parser <- add_argument(parser, "--y_bias", help='coefficient of y on the probability a classification is correct', default=-0.5)
parser <- add_argument(parser, "--z_bias", help='coefficient of z on the probability a classification is correct', default=0)
parser <- add_argument(parser, "--truth_formula", help='formula for the true variable', default="x~z")
args <- parse_args(parser)
B0 <- 0
Bxy <- 0.3
Bxy <- args$Bxy
Bzy <- args$Bzy
Bzx <- args$Bzx
if(args$m < args$N){
df <- simulate_data(args$N, args$m, B0, Bxy, args$Bzx, Bzy, args$seed, args$y_explained_variance, args$prediction_accuracy, args$accuracy_imbalance_difference)
result <- list('N'=args$N,'m'=args$m,'B0'=B0,'Bxy'=Bxy, Bzx=args$Bzx, 'Bzy'=Bzy, 'seed'=args$seed, 'y_explained_variance'=args$y_explained_variance, 'prediction_accuracy'=args$prediction_accuracy, 'accuracy_imbalance_difference'=args$accuracy_imbalance_difference, error='')
df <- simulate_data(args$N, args$m, B0, Bxy, Bzx, Bzy, args$seed, args$y_explained_variance, args$prediction_accuracy, y_bias=args$y_bias)
outline <- run_simulation(df, result, outcome_formula=y~x+z, proxy_formula=w_pred~x+z+y+x:y, truth_formula=x~z)
## df.pc <- df[,.(x,y,z,w_pred,w)]
## # df.pc <- df.pc[,err:=x-w_pred]
## pc.df <- pc(suffStat=list(C=cor(df.pc),n=nrow(df.pc)),indepTest=gaussCItest,labels=names(df.pc),alpha=0.05)
## plot(pc.df)
outfile_lock <- lock(paste0(args$outfile, '_lock'),exclusive=TRUE)
result <- list('N'=args$N,'m'=args$m,'B0'=B0,'Bxy'=Bxy, Bzx=args$Bzx, 'Bzy'=Bzy, 'seed'=args$seed, 'y_explained_variance'=args$y_explained_variance, 'prediction_accuracy'=args$prediction_accuracy, 'accuracy_imbalance_difference'=args$accuracy_imbalance_difference, 'y_bias'=args$y_bias,'outcome_formula'=args$outcome_formula, 'proxy_formula'=args$proxy_formula,truth_formula=args$truth_formula, error='')
outline <- run_simulation(df, result, outcome_formula=as.formula(args$outcome_formula), proxy_formula=as.formula(args$proxy_formula), truth_formula=as.formula(args$truth_formula))
outfile_lock <- lock(paste0(args$outfile, '_lock'),exclusive=TRUE)
if(file.exists(args$outfile)){
logdata <- read_feather(args$outfile)
logdata <- rbind(logdata,as.data.table(outline), fill=TRUE)

109
simulations/03_depvar.R Normal file
View File

@@ -0,0 +1,109 @@
### EXAMPLE 1: demonstrates how measurement error can lead to a type sign error in a covariate
### What kind of data invalidates fong + tyler?
### Even when you have a good predictor, if it's biased against a covariate you can get the wrong sign.
### Even when you include the proxy variable in the regression.
### But with some ground truth and multiple imputation, you can fix it.
library(argparser)
library(mecor)
library(ggplot2)
library(data.table)
library(filelock)
library(arrow)
library(Amelia)
library(Zelig)
library(predictionError)
options(amelia.parallel="no",
amelia.ncpus=1)
setDTthreads(40)
source("simulation_base.R")
## SETUP:
### we want to estimate x -> y; x is MAR
### we have x -> k; k -> w; x -> w is used to predict x via the model w.
### A realistic scenario is that we have an NLP model predicting something like "racial harassment" in social media comments
### The labels x are binary, but the model provides a continuous predictor
### simulation:
#### how much power do we get from the model in the first place? (sweeping N and m)
####
## one way to do it is by adding correlation to x.obs and y that isn't in w.
## in other words, the model is missing an important feature of x.obs that's related to y.
simulate_data <- function(N, m, B0, Bxy, Bzy, Bzx, seed, prediction_accuracy=0.73, log.likelihood.gain = 0.1){
set.seed(seed)
set.seed(seed)
# make w and y dependent
z <- rnorm(N, sd=0.5)
x <- rbinom(N, 1, plogis(Bzx*z))
ystar <- Bzy * z + Bxy * x + B0
y <- rbinom(N,1,plogis(ystar))
df <- data.table(x=x,y=y,ystar=ystar,z=z)
if(m < N){
df <- df[sample(nrow(df), m), y.obs := y]
} else {
df <- df[, y.obs := y]
}
odds.y1 <- qlogis(prediction_accuracy)
odds.y0 <- qlogis(prediction_accuracy,lower.tail=F)
df[y==0,w:=plogis(rlogis(.N,odds.y0))]
df[y==1,w:=plogis(rlogis(.N,odds.y1))]
df[,w_pred := as.integer(w > 0.5)]
print(mean(df[x==0]$y == df[x==0]$w_pred))
print(mean(df[x==1]$y == df[x==1]$w_pred))
print(mean(df$w_pred == df$y))
return(df)
}
parser <- arg_parser("Simulate data and fit corrected models")
parser <- add_argument(parser, "--N", default=10000, help="number of observations of w")
parser <- add_argument(parser, "--m", default=500, help="m the number of ground truth observations")
parser <- add_argument(parser, "--seed", default=17, help='seed for the rng')
parser <- add_argument(parser, "--outfile", help='output file', default='example_2.feather')
parser <- add_argument(parser, "--y_explained_variance", help='what proportion of the variance of y can be explained?', default=0.1)
parser <- add_argument(parser, "--prediction_accuracy", help='how accurate is the predictive model?', default=0.72)
## parser <- add_argument(parser, "--x_bias_y1", help='how is the classifier biased when y = 1?', default=-0.75)
## parser <- add_argument(parser, "--x_bias_y0", help='how is the classifier biased when y = 0 ?', default=0.75)
parser <- add_argument(parser, "--Bxy", help='coefficient of x on y', default=0.01)
parser <- add_argument(parser, "--Bzy", help='coeffficient of z on y', default=-0.01)
parser <- add_argument(parser, "--Bzx", help='coeffficient of z on x', default=-0.5)
parser <- add_argument(parser, "--outcome_formula", help='formula for the outcome variable', default="y~x+z")
parser <- add_argument(parser, "--proxy_formula", help='formula for the proxy variable', default="w_pred~y")
args <- parse_args(parser)
B0 <- 0
Bxy <- args$Bxy
Bzy <- args$Bzy
Bzx <- args$Bzx
if(args$m < args$N){
df <- simulate_data(args$N, args$m, B0, Bxy, Bzy, Bzx, args$seed, args$prediction_accuracy)
# result <- list('N'=args$N,'m'=args$m,'B0'=B0,'Bxy'=Bxy,'Bzy'=Bzy, 'seed'=args$seed, 'y_explained_variance'=args$y_explained_variance, 'prediction_accuracy'=args$prediction_accuracy, 'x_bias_y0'=args$x_bias_y0,'x_bias_y1'=args$x_bias_y1,'outcome_formula' = args$outcome_formula, 'proxy_formula' = args$proxy_formula)
result <- list('N'=args$N,'m'=args$m,'B0'=B0,'Bxy'=Bxy,'Bzy'=Bzy, 'Bzx'=Bzx,'seed'=args$seed, 'y_explained_variance'=args$y_explained_variance, 'prediction_accuracy'=args$prediction_accuracy, 'outcome_formula' = args$outcome_formula, 'proxy_formula' = args$proxy_formula)
outline <- run_simulation_depvar(df, result, outcome_formula = as.formula(args$outcome_formula), proxy_formula = as.formula(args$proxy_formula))
outfile_lock <- lock(paste0(args$outfile, '_lock'),exclusive=TRUE)
if(file.exists(args$outfile)){
logdata <- read_feather(args$outfile)
logdata <- rbind(logdata,as.data.table(outline),fill=TRUE)
} else {
logdata <- as.data.table(outline)
}
print(outline)
write_feather(logdata, args$outfile)
unlock(outfile_lock)
}

63
simulations/03_depvar_differential.R
View File

@@ -31,13 +31,14 @@ source("simulation_base.R")
## one way to do it is by adding correlation to x.obs and y that isn't in w.
## in other words, the model is missing an important feature of x.obs that's related to y.
simulate_data <- function(N, m, B0, Bxy, Bzy, seed, prediction_accuracy=0.73, accuracy_imbalance_difference=0.3){
simulate_data <- function(N, m, B0, Bxy, Bzy, seed, prediction_accuracy=0.73, x_bias=-0.75){
set.seed(seed)
# make w and y dependent
z <- rbinom(N, 1, 0.5)
x <- rbinom(N, 1, 0.5)
ystar <- Bzy * z + Bxy * x
ystar <- Bzy * z + Bxy * x + B0
y <- rbinom(N,1,plogis(ystar))
# glm(y ~ x + z, family="binomial")
@@ -50,39 +51,17 @@ simulate_data <- function(N, m, B0, Bxy, Bzy, seed, prediction_accuracy=0.73, ac
df <- df[, y.obs := y]
}
df <- df[,w_pred:=y]
odds.y1 <- qlogis(prediction_accuracy) + x_bias*df[y==1]$x
odds.y0 <- qlogis(prediction_accuracy,lower.tail=F) + x_bias*df[y==0]$x
pz <- mean(z)
df[y==0,w:=plogis(rlogis(.N,odds.y0))]
df[y==1,w:=plogis(rlogis(.N,odds.y1))]
accuracy_imbalance_ratio <- (prediction_accuracy + accuracy_imbalance_difference/2) / (prediction_accuracy - accuracy_imbalance_difference/2)
df[,w_pred := as.integer(w > 0.5)]
# this works because of conditional probability
accuracy_z0 <- prediction_accuracy / (pz*(accuracy_imbalance_ratio) + (1-pz))
accuracy_z1 <- accuracy_imbalance_ratio * accuracy_z0
yz0 <- df[z==0]$y
yz1 <- df[z==1]$y
nz1 <- nrow(df[z==1])
nz0 <- nrow(df[z==0])
acc_z0 <- plogis(0.7*scale(yz0) + qlogis(accuracy_z0))
acc_z1 <- plogis(1.3*scale(yz1) + qlogis(accuracy_z1))
w0z0 <- (1-yz0)**2 + (-1)**(1-yz0) * acc_z0
w0z1 <- (1-yz1)**2 + (-1)**(1-yz1) * acc_z1
w0z0.noisy.odds <- rlogis(nz0,qlogis(w0z0))
w0z1.noisy.odds <- rlogis(nz1,qlogis(w0z1))
df[z==0,w:=plogis(w0z0.noisy.odds)]
df[z==1,w:=plogis(w0z1.noisy.odds)]
df[,w_pred:=as.integer(w > 0.5)]
print(mean(df[y==0]$y == df[y==0]$w_pred))
print(mean(df[y==1]$y == df[y==1]$w_pred))
print(mean(df[x==0]$y == df[x==0]$w_pred))
print(mean(df[x==1]$y == df[x==1]$w_pred))
print(mean(df$w_pred == df$y))
return(df)
}
@@ -92,21 +71,29 @@ parser <- add_argument(parser, "--m", default=500, help="m the number of ground
parser <- add_argument(parser, "--seed", default=17, help='seed for the rng')
parser <- add_argument(parser, "--outfile", help='output file', default='example_2.feather')
parser <- add_argument(parser, "--y_explained_variance", help='what proportion of the variance of y can be explained?', default=0.005)
parser <- add_argument(parser, "--prediction_accuracy", help='how accurate is the predictive model?', default=0.73)
parser <- add_argument(parser, "--accuracy_imbalance_difference", help='how much more accurate is the predictive model for one class than the other?', default=0.3)
parser <- add_argument(parser, "--prediction_accuracy", help='how accurate is the predictive model?', default=0.8)
## parser <- add_argument(parser, "--x_bias_y1", help='how is the classifier biased when y = 1?', default=-0.75)
## parser <- add_argument(parser, "--x_bias_y0", help='how is the classifier biased when y = 0 ?', default=0.75)
parser <- add_argument(parser, "--x_bias", help='how is the classifier biased?', default=0.75)
parser <- add_argument(parser, "--Bxy", help='coefficient of x on y', default=0.3)
parser <- add_argument(parser, "--Bzy", help='coeffficient of z on y', default=-0.3)
parser <- add_argument(parser, "--outcome_formula", help='formula for the outcome variable', default="y~x+z")
parser <- add_argument(parser, "--proxy_formula", help='formula for the proxy variable', default="w_pred~y*x")
args <- parse_args(parser)
B0 <- 0
Bxy <- 0.7
Bzy <- -0.7
Bxy <- args$Bxy
Bzy <- args$Bzy
if(args$m < args$N){
df <- simulate_data(args$N, args$m, B0, Bxy, Bzy, args$seed, args$prediction_accuracy, args$accuracy_imbalance_difference)
df <- simulate_data(args$N, args$m, B0, Bxy, Bzy, args$seed, args$prediction_accuracy, args$x_bias_y0, args$x_bias_y1)
result <- list('N'=args$N,'m'=args$m,'B0'=B0,'Bxy'=Bxy,'Bzy'=Bzy, 'seed'=args$seed, 'y_explained_variance'=args$y_explained_variance, 'prediction_accuracy'=args$prediction_accuracy, 'accuracy_imbalance_difference'=args$accuracy_imbalance_difference)
# result <- list('N'=args$N,'m'=args$m,'B0'=B0,'Bxy'=Bxy,'Bzy'=Bzy, 'seed'=args$seed, 'y_explained_variance'=args$y_explained_variance, 'prediction_accuracy'=args$prediction_accuracy, 'x_bias_y0'=args$x_bias_y0,'x_bias_y1'=args$x_bias_y1,'outcome_formula' = args$outcome_formula, 'proxy_formula' = args$proxy_formula)
result <- list('N'=args$N,'m'=args$m,'B0'=B0,'Bxy'=Bxy,'Bzy'=Bzy, 'seed'=args$seed, 'y_explained_variance'=args$y_explained_variance, 'prediction_accuracy'=args$prediction_accuracy, 'x_bias'=args$x_bias,'x_bias'=args$x_bias,'outcome_formula' = args$outcome_formula, 'proxy_formula' = args$proxy_formula)
outline <- run_simulation_depvar(df, result, outcome_formula = y ~ x + z, proxy_formula = w_pred ~ y*x + y*z + z*x)
outline <- run_simulation_depvar(df, result, outcome_formula = as.formula(args$outcome_formula), proxy_formula = as.formula(args$proxy_formula))
outfile_lock <- lock(paste0(args$outfile, '_lock'),exclusive=TRUE)

109
simulations/04_depvar_differential.R Normal file
View File

@@ -0,0 +1,109 @@
### EXAMPLE 1: demonstrates how measurement error can lead to a type sign error in a covariate
### What kind of data invalidates fong + tyler?
### Even when you have a good predictor, if it's biased against a covariate you can get the wrong sign.
### Even when you include the proxy variable in the regression.
### But with some ground truth and multiple imputation, you can fix it.
library(argparser)
library(mecor)
library(ggplot2)
library(data.table)
library(filelock)
library(arrow)
library(Amelia)
library(Zelig)
library(predictionError)
options(amelia.parallel="no",
amelia.ncpus=1)
setDTthreads(40)
source("simulation_base.R")
## SETUP:
### we want to estimate x -> y; x is MAR
### we have x -> k; k -> w; x -> w is used to predict x via the model w.
### A realistic scenario is that we have an NLP model predicting something like "racial harassment" in social media comments
### The labels x are binary, but the model provides a continuous predictor
### simulation:
#### how much power do we get from the model in the first place? (sweeping N and m)
####
## one way to do it is by adding correlation to x.obs and y that isn't in w.
## in other words, the model is missing an important feature of x.obs that's related to y.
simulate_data <- function(N, m, B0, Bxy, Bzy, seed, prediction_accuracy=0.73, z_bias=-0.75){
set.seed(seed)
# make w and y dependent
z <- rnorm(N,sd=0.5)
x <- rbinom(N,1,0.5)
ystar <- Bzy * z + Bxy * x + B0
y <- rbinom(N,1,plogis(ystar))
# glm(y ~ x + z, family="binomial")
df <- data.table(x=x,y=y,ystar=ystar,z=z)
if(m < N){
df <- df[sample(nrow(df), m), y.obs := y]
} else {
df <- df[, y.obs := y]
}
odds.y1 <- qlogis(prediction_accuracy) + z_bias*df[y==1]$z
odds.y0 <- qlogis(prediction_accuracy,lower.tail=F) + z_bias*df[y==0]$z
df[y==0,w:=plogis(rlogis(.N,odds.y0))]
df[y==1,w:=plogis(rlogis(.N,odds.y1))]
df[,w_pred := as.integer(w > 0.5)]
print(mean(df[x==0]$y == df[x==0]$w_pred))
print(mean(df[x==1]$y == df[x==1]$w_pred))
print(mean(df$w_pred == df$y))
return(df)
}
parser <- arg_parser("Simulate data and fit corrected models")
parser <- add_argument(parser, "--N", default=1000, help="number of observations of w")
parser <- add_argument(parser, "--m", default=500, help="m the number of ground truth observations")
parser <- add_argument(parser, "--seed", default=17, help='seed for the rng')
parser <- add_argument(parser, "--outfile", help='output file', default='example_4.feather')
parser <- add_argument(parser, "--prediction_accuracy", help='how accurate is the predictive model?', default=0.79)
## parser <- add_argument(parser, "--z_bias_y1", help='how is the classifier biased when y = 1?', default=-0.75)
## parser <- add_argument(parser, "--z_bias_y0", help='how is the classifier biased when y = 0 ?', default=0.75)
parser <- add_argument(parser, "--z_bias", help='how is the classifier biased?', default=1.5)
parser <- add_argument(parser, "--Bxy", help='coefficient of x on y', default=0.1)
parser <- add_argument(parser, "--Bzy", help='coeffficient of z on y', default=-0.1)
parser <- add_argument(parser, "--outcome_formula", help='formula for the outcome variable', default="y~x+z")
parser <- add_argument(parser, "--proxy_formula", help='formula for the proxy variable', default="w_pred~y+z")
args <- parse_args(parser)
B0 <- 0
Bxy <- args$Bxy
Bzy <- args$Bzy
if(args$m < args$N){
df <- simulate_data(args$N, args$m, B0, Bxy, Bzy, args$seed, args$prediction_accuracy, args$z_bias)
# result <- list('N'=args$N,'m'=args$m,'B0'=B0,'Bxy'=Bxy,'Bzy'=Bzy, 'seed'=args$seed, 'y_explained_variance'=args$y_explained_variance, 'prediction_accuracy'=args$prediction_accuracy, 'z_bias_y0'=args$z_bias_y0,'z_bias_y1'=args$z_bias_y1,'outcome_formula' = args$outcome_formula, 'proxy_formula' = args$proxy_formula)
result <- list('N'=args$N,'m'=args$m,'B0'=B0,'Bxy'=Bxy,'Bzy'=Bzy, 'seed'=args$seed, 'y_explained_variance'=args$y_explained_variance, 'prediction_accuracy'=args$prediction_accuracy, 'z_bias'=args$z_bias,'outcome_formula' = args$outcome_formula, 'proxy_formula' = args$proxy_formula)
outline <- run_simulation_depvar(df, result, outcome_formula = as.formula(args$outcome_formula), proxy_formula = as.formula(args$proxy_formula))
outfile_lock <- lock(paste0(args$outfile, '_lock'),exclusive=TRUE)
if(file.exists(args$outfile)){
logdata <- read_feather(args$outfile)
logdata <- rbind(logdata,as.data.table(outline),fill=TRUE)
} else {
logdata <- as.data.table(outline)
}
print(outline)
write_feather(logdata, args$outfile)
unlock(outfile_lock)
}

116
simulations/05_irr_indep.R Normal file
View File

@@ -0,0 +1,116 @@
### EXAMPLE 2_b: demonstrates how measurement error can lead to a type
### sign error in a covariate This is the same as example 2, only
### instead of x->k we have k->x. Even when you have a good
### predictor, if it's biased against a covariate you can get the
### wrong sign. Even when you include the proxy variable in the
### regression. But with some ground truth and multiple imputation,
### you can fix it.
library(argparser)
library(mecor)
library(ggplot2)
library(data.table)
library(filelock)
library(arrow)
library(Amelia)
library(Zelig)
library(predictionError)
options(amelia.parallel="no", amelia.ncpus=1)
source("irr_simulation_base.R")
## SETUP:
### we want to estimate x -> y; x is MAR
### we have x -> k; k -> w; x -> w is used to predict x via the model w.
### A realistic scenario is that we have an NLP model predicting something like "racial harassment" in social media comments
### The labels x are binary, but the model provides a continuous predictor
### simulation:
#### how much power do we get from the model in the first place? (sweeping N and m)
####
simulate_data <- function(N, m, B0=0, Bxy=0.2, Bzy=-0.2, Bzx=0.2, y_explained_variance=0.025, prediction_accuracy=0.73, coder_accuracy=0.9, seed=1){
set.seed(seed)
z <- rbinom(N, 1, 0.5)
# x.var.epsilon <- var(Bzx *z) * ((1-zx_explained_variance)/zx_explained_variance)
xprime <- Bzx * z #+ x.var.epsilon
x <- rbinom(N,1,plogis(xprime))
y.var.epsilon <- (var(Bzy * z) + var(Bxy *x) + 2*cov(Bxy*x,Bzy*z)) * ((1-y_explained_variance)/y_explained_variance)
y.epsilon <- rnorm(N, sd = sqrt(y.var.epsilon))
y <- Bzy * z + Bxy * x + y.epsilon + B0
df <- data.table(x=x,y=y,z=z)
if(m < N){
df <- df[sample(nrow(df), m), x.obs := x]
} else {
df <- df[, x.obs := x]
}
coder.0.correct <- rbinom(m, 1, coder_accuracy)
coder.1.correct <- rbinom(m, 1, coder_accuracy)
df[!is.na(x.obs),x.obs.0 := as.numeric((x.obs & coder.0.correct) | (!x.obs & !coder.0.correct))]
df[!is.na(x.obs),x.obs.1 := as.numeric((x.obs & coder.1.correct) | (!x.obs & !coder.1.correct))]
## how can you make a model with a specific accuracy?
w0 =(1-x)**2 + (-1)**(1-x) * prediction_accuracy
## how can you make a model with a specific accuracy, with a continuous latent variable.
# now it makes the same amount of mistake to each point, probably
# add mean0 noise to the odds.
w.noisey.odds = rlogis(N,qlogis(w0))
df[,w := plogis(w.noisey.odds)]
df[,w_pred:=as.integer(w > 0.5)]
(mean(df$x==df$w_pred))
return(df)
}
parser <- arg_parser("Simulate data and fit corrected models")
parser <- add_argument(parser, "--N", default=1000, help="number of observations of w")
parser <- add_argument(parser, "--m", default=150, help="m the number of ground truth observations")
parser <- add_argument(parser, "--seed", default=1, help='seed for the rng')
parser <- add_argument(parser, "--outfile", help='output file', default='example_1.feather')
parser <- add_argument(parser, "--y_explained_variance", help='what proportion of the variance of y can be explained?', default=0.1)
# parser <- add_argument(parser, "--zx_explained_variance", help='what proportion of the variance of x can be explained by z?', default=0.3)
parser <- add_argument(parser, "--prediction_accuracy", help='how accurate is the predictive model?', default=0.73)
parser <- add_argument(parser, "--coder_accuracy", help='how accurate are the human coders?', default=0.85)
parser <- add_argument(parser, "--outcome_formula", help='formula for the outcome variable', default="y~x+z")
parser <- add_argument(parser, "--proxy_formula", help='formula for the proxy variable', default="w_pred~x")
# parser <- add_argument(parser, "--rater_formula", help='formula for the true variable', default="x.obs~x")
parser <- add_argument(parser, "--truth_formula", help='formula for the true variable', default="x~z")
parser <- add_argument(parser, "--Bzx", help='Effect of z on x', default=-0.3)
parser <- add_argument(parser, "--Bzy", help='Effect of z on y', default=0.27)
parser <- add_argument(parser, "--Bxy", help='Effect of x on y', default=-0.33)
args <- parse_args(parser)
B0 <- 0
Bxy <- args$Bxy
Bzy <- args$Bzy
Bzx <- args$Bzx
if (args$m < args$N){
df <- simulate_data(args$N, args$m, B0, Bxy, Bzy, Bzx, seed=args$seed, y_explained_variance = args$y_explained_variance, prediction_accuracy=args$prediction_accuracy, coder_accuracy=args$coder_accuracy)
result <- list('N'=args$N,'m'=args$m,'B0'=B0,'Bxy'=Bxy, Bzx=Bzx, 'Bzy'=Bzy, 'seed'=args$seed, 'y_explained_variance'=args$y_explained_variance, 'prediction_accuracy'=args$prediction_accuracy, 'accuracy_imbalance_difference'=args$accuracy_imbalance_difference, 'outcome_formula'=args$outcome_formula, 'truth_formula'=args$truth_formula, 'proxy_formula'=args$proxy_formula,truth_formula=args$truth_formula, 'coder_accuracy'=args$coder_accuracy, error='')
outline <- run_simulation(df, result, outcome_formula=as.formula(args$outcome_formula), proxy_formula=as.formula(args$proxy_formula), truth_formula=as.formula(args$truth_formula))
outfile_lock <- lock(paste0(args$outfile, '_lock'),exclusive=TRUE)
if(file.exists(args$outfile)){
logdata <- read_feather(args$outfile)
logdata <- rbind(logdata,as.data.table(outline),fill=TRUE)
} else {
logdata <- as.data.table(outline)
}
print(outline)
write_feather(logdata, args$outfile)
unlock(outfile_lock)
}

101
simulations/06_irr_dv.R Normal file
View File

@@ -0,0 +1,101 @@
library(argparser)
library(mecor)
library(ggplot2)
library(data.table)
library(filelock)
library(arrow)
library(Amelia)
library(Zelig)
library(predictionError)
options(amelia.parallel="no",
amelia.ncpus=1)
setDTthreads(40)
source("irr_dv_simulation_base.R")
## one way to do it is by adding correlation to x.obs and y that isn't in w.
## in other words, the model is missing an important feature of x.obs that's related to y.
simulate_data <- function(N, m, B0, Bxy, Bzy, seed, prediction_accuracy=0.73, coder_accuracy=0.8){
set.seed(seed)
# make w and y dependent
z <- rbinom(N, 1, 0.5)
x <- rbinom(N, 1, 0.5)
ystar <- Bzy * z + Bxy * x + B0
y <- rbinom(N,1,plogis(ystar))
# glm(y ~ x + z, family="binomial")
df <- data.table(x=x,y=y,ystar=ystar,z=z)
df <- df[sample(nrow(df), m), y.obs := y]
coder.0.correct <- rbinom(m, 1, coder_accuracy)
coder.1.correct <- rbinom(m, 1, coder_accuracy)
df[!is.na(y.obs),y.obs.0 := as.numeric((.SD$y.obs & coder.0.correct) | (!.SD$y.obs & !coder.0.correct))]
df[!is.na(y.obs),y.obs.1 := as.numeric((.SD$y.obs & coder.1.correct) | (!.SD$y.obs & !coder.1.correct))]
odds.y1 <- qlogis(prediction_accuracy)
odds.y0 <- qlogis(prediction_accuracy,lower.tail=F)
df[y==0,w:=plogis(rlogis(.N,odds.y0))]
df[y==1,w:=plogis(rlogis(.N,odds.y1))]
df[,w_pred := as.integer(w > 0.5)]
print(mean(df$y == df$y.obs.0,na.rm=T))
print(mean(df$y == df$y.obs.1,na.rm=T))
print(mean(df[x==0]$y == df[x==0]$w_pred))
print(mean(df[x==1]$y == df[x==1]$w_pred))
print(mean(df$w_pred == df$y))
return(df)
}
parser <- arg_parser("Simulate data and fit corrected models")
parser <- add_argument(parser, "--N", default=5000, help="number of observations of w")
parser <- add_argument(parser, "--m", default=500, help="m the number of ground truth observations")
parser <- add_argument(parser, "--seed", default=16, help='seed for the rng')
parser <- add_argument(parser, "--outfile", help='output file', default='example_2.feather')
parser <- add_argument(parser, "--y_explained_variance", help='what proportion of the variance of y can be explained?', default=0.1)
parser <- add_argument(parser, "--prediction_accuracy", help='how accurate is the predictive model?', default=0.73)
## parser <- add_argument(parser, "--x_bias_y1", help='how is the classifier biased when y = 1?', default=-0.75)
## parser <- add_argument(parser, "--x_bias_y0", help='how is the classifier biased when y = 0 ?', default=0.75)
parser <- add_argument(parser, "--Bxy", help='coefficient of x on y', default=0.3)
parser <- add_argument(parser, "--Bzy", help='coeffficient of z on y', default=-0.3)
parser <- add_argument(parser, "--outcome_formula", help='formula for the outcome variable', default="y~x+z")
parser <- add_argument(parser, "--proxy_formula", help='formula for the proxy variable', default="w_pred~y+y.obs.1+y.obs.0")
parser <- add_argument(parser, "--coder_accuracy", help='How accurate are the coders?', default=0.8)
args <- parse_args(parser)
B0 <- 0
Bxy <- args$Bxy
Bzy <- args$Bzy
df <- simulate_data(args$N, args$m, B0, Bxy, Bzy, args$seed, args$prediction_accuracy, args$coder_accuracy)
# result <- list('N'=args$N,'m'=args$m,'B0'=B0,'Bxy'=Bxy,'Bzy'=Bzy, 'seed'=args$seed, 'y_explained_variance'=args$y_explained_variance, 'prediction_accuracy'=args$prediction_accuracy, 'x_bias_y0'=args$x_bias_y0,'x_bias_y1'=args$x_bias_y1,'outcome_formula' = args$outcome_formula, 'proxy_formula' = args$proxy_formula)
result <- list('N'=args$N,'m'=args$m,'B0'=B0,'Bxy'=Bxy,'Bzy'=Bzy, 'seed'=args$seed, 'y_explained_variance'=args$y_explained_variance, 'prediction_accuracy'=args$prediction_accuracy, 'outcome_formula' = args$outcome_formula, 'proxy_formula' = args$proxy_formula)
outline <- run_simulation_depvar(df, result, outcome_formula = as.formula(args$outcome_formula), proxy_formula = as.formula(args$proxy_formula))
outfile_lock <- lock(paste0(args$outfile, '_lock'),exclusive=TRUE)
if(file.exists(args$outfile)){
logdata <- read_feather(args$outfile)
logdata <- rbind(logdata,as.data.table(outline),fill=TRUE)
} else {
logdata <- as.data.table(outline)
}
print(outline)
write_feather(logdata, args$outfile)
unlock(outfile_lock)
warnings()

127
simulations/Makefile
View File

@@ -1,20 +1,21 @@
SHELL=bash
Ns=[1000,3600,14400]
ms=[75,150,300]
seeds=[$(shell seq -s, 1 250)]
Ns=[1000, 5000, 10000]
ms=[100, 200, 400]
seeds=[$(shell seq -s, 1 500)]
explained_variances=[0.1]
all:remembr.RDS
all:remembr.RDS remember_irr.RDS
supplement: remember_robustness_misspec.RDS
srun=srun -A comdata -p compute-bigmem --time=6:00:00 --mem 4G -c 1
srun=sbatch --wait --verbose run_job.sbatch
joblists:example_1_jobs example_2_jobs example_3_jobs
# test_true_z_jobs: test_true_z.R simulation_base.R
# grid_sweep.py --command "Rscript test_true_z.R" --arg_dict '{"N":${Ns},"m":${ms}, "seed":${seeds}, "outfile":["test_true_z.feather"], "y_explained_variancevari":${explained_variances}, "Bzx":${Bzx}}' --outfile test_true_z_jobsb
# sbatch --wait --verbose run_job.sbatch grid_sweep.py --command "Rscript test_true_z.R" --arg_dict '{"N":${Ns},"m":${ms}, "seed":${seeds}, "outfile":["test_true_z.feather"], "y_explained_variancevari":${explained_variances}, "Bzx":${Bzx}}' --outfile test_true_z_jobsb
# test_true_z.feather: test_true_z_jobs
# rm -f test_true_z.feather
@@ -22,42 +23,127 @@ joblists:example_1_jobs example_2_jobs example_3_jobs
# sbatch --wait --verbose --array=3001-6001 run_simulation.sbatch 0 test_true_z_jobs
example_1_jobs: 01_two_covariates.R simulation_base.R
grid_sweep.py --command "Rscript 01_two_covariates.R" --arg_dict '{"N":${Ns},"m":${ms}, "seed":${seeds}, "outfile":["example_1.feather"], "y_explained_variance":${explained_variances}, "Bzx":[0.1]}' --outfile example_1_jobs
example_1_jobs: 01_two_covariates.R simulation_base.R grid_sweep.py pl_methods.R
sbatch --wait --verbose run_job.sbatch grid_sweep.py --command "Rscript 01_two_covariates.R" --arg_dict '{"N":${Ns},"m":${ms}, "seed":${seeds}, "outfile":["example_1.feather"], "y_explained_variance":${explained_variances}, "Bzx":[1]}' --outfile example_1_jobs
example_1.feather: example_1_jobs
rm -f example_1.feather
sbatch --wait --verbose --array=1-$(shell cat example_1_jobs | wc -l) run_simulation.sbatch 0 example_1_jobs
# sbatch --wait --verbose --array=3001-6001 run_simulation.sbatch 0 example_1_jobs
sbatch --wait --verbose --array=1-1000 run_simulation.sbatch 0 example_1_jobs
sbatch --wait --verbose --array=1001-2000 run_simulation.sbatch 0 example_1_jobs
sbatch --wait --verbose --array=2001-3000 run_simulation.sbatch 0 example_1_jobs
sbatch --wait --verbose --array=3001-4000 run_simulation.sbatch 0 example_1_jobs
sbatch --wait --verbose --array=4001-$(shell cat example_1_jobs | wc -l) run_simulation.sbatch 0 example_1_jobs
example_2_jobs: 02_indep_differential.R simulation_base.R
grid_sweep.py --command "Rscript 02_indep_differential.R" --arg_dict '{"N":${Ns},"m":${ms}, "seed":${seeds}, "outfile":["example_2.feather"],"y_explained_variance":${explained_variances}, "accuracy_imbalance_difference":[0.3], "Bzy":[0.3]}' --outfile example_2_jobs
example_2_jobs: 02_indep_differential.R simulation_base.R grid_sweep.py pl_methods.R
sbatch --wait --verbose run_job.sbatch grid_sweep.py --command "Rscript 02_indep_differential.R" --arg_dict '{"N":${Ns},"m":${ms}, "seed":${seeds}, "outfile":["example_2.feather"],"y_explained_variance":${explained_variances}, "Bzy":[-0.3],"Bxy":[0.3],"Bzx":[1], "outcome_formula":["y~x+z"], "proxy_formula":["w_pred~y*z*x"]}' --outfile example_2_jobs
example_2.feather: example_2_jobs
rm -f example_2.feather
sbatch --wait --verbose --array=1-$(shell cat example_2_jobs | wc -l) run_simulation.sbatch 0 example_2_jobs
# sbatch --wait --verbose --array=3001-6001 run_simulation.sbatch 0 example_2_jobs
sbatch --wait --verbose --array=1-1000 run_simulation.sbatch 0 example_2_jobs
sbatch --wait --verbose --array=1001-2000 run_simulation.sbatch 0 example_2_jobs
sbatch --wait --verbose --array=2001-3000 run_simulation.sbatch 0 example_2_jobs
sbatch --wait --verbose --array=3001-4000 run_simulation.sbatch 0 example_2_jobs
sbatch --wait --verbose --array=4001-$(shell cat example_2_jobs | wc -l) run_simulation.sbatch 0 example_2_jobs
# example_2_B_jobs: example_2_B.R
# grid_sweep.py --command "Rscript example_2_B.R" --arg_dict '{"N":${Ns},"m":${ms}, "seed":${seeds}, "outfile":["example_2_B.feather"]}' --outfile example_2_B_jobs
# sbatch --wait --verbose run_job.sbatch grid_sweep.py --command "Rscript example_2_B.R" --arg_dict '{"N":${Ns},"m":${ms}, "seed":${seeds}, "outfile":["example_2_B.feather"]}' --outfile example_2_B_jobs
# example_2_B.feather: example_2_B_jobs
# rm -f example_2_B.feather
# sbatch --wait --verbose --array=1-3000 run_simulation.sbatch 0 example_2_B_jobs
example_3_jobs: 03_depvar_differential.R simulation_base.R
grid_sweep.py --command "Rscript 03_depvar_differential.R" --arg_dict '{"N":${Ns},"m":${ms}, "seed":${seeds}, "outfile":["example_3.feather"], "y_explained_variance":${explained_variances}}' --outfile example_3_jobs
example_3_jobs: 03_depvar.R simulation_base.R grid_sweep.py pl_methods.R
sbatch --wait --verbose run_job.sbatch grid_sweep.py --command "Rscript 03_depvar.R" --arg_dict '{"N":${Ns},"m":${ms}, "Bxy":[0.7],"Bzy":[-0.7],"seed":${seeds}, "outfile":["example_3.feather"], "y_explained_variance":${explained_variances}}' --outfile example_3_jobs
example_3.feather: example_3_jobs
rm -f example_3.feather
sbatch --wait --verbose --array=1-$(shell cat example_3_jobs | wc -l) run_simulation.sbatch 0 example_3_jobs
sbatch --wait --verbose --array=1-1000 run_simulation.sbatch 0 example_3_jobs
sbatch --wait --verbose --array=1001-2000 run_simulation.sbatch 0 example_3_jobs
sbatch --wait --verbose --array=2001-3000 run_simulation.sbatch 0 example_3_jobs
sbatch --wait --verbose --array=3001-4000 run_simulation.sbatch 0 example_3_jobs
sbatch --wait --verbose --array=4001-$(shell cat example_3_jobs | wc -l) run_simulation.sbatch 0 example_3_jobs
example_4_jobs: 04_depvar_differential.R simulation_base.R grid_sweep.py pl_methods.R
sbatch --wait --verbose run_job.sbatch grid_sweep.py --command "Rscript 04_depvar_differential.R" --arg_dict '{"N":${Ns},"Bxy":[0.7],"Bzy":[-0.7],"m":${ms}, "seed":${seeds}, "outfile":["example_4.feather"], "z_bias":[0.5]}' --outfile example_4_jobs
example_4.feather: example_4_jobs
rm -f example_4.feather
sbatch --wait --verbose --array=1-1000 run_simulation.sbatch 0 example_4_jobs
sbatch --wait --verbose --array=1001-2000 run_simulation.sbatch 0 example_4_jobs
sbatch --wait --verbose --array=2001-3000 run_simulation.sbatch 0 example_4_jobs
sbatch --wait --verbose --array=3001-4000 run_simulation.sbatch 0 example_4_jobs
sbatch --wait --verbose --array=4001-$(shell cat example_4_jobs | wc -l) run_simulation.sbatch 0 example_4_jobs
remembr.RDS:example_1.feather example_2.feather example_3.feather plot_example.R plot_dv_example.R
remembr.RDS:example_1.feather example_2.feather example_3.feather example_4.feather plot_example.R plot_dv_example.R summarize_estimator.R
rm -f remembr.RDS
${srun} Rscript plot_example.R --infile example_1.feather --name "plot.df.example.1"
${srun} Rscript plot_example.R --infile example_2.feather --name "plot.df.example.2"
${srun} Rscript plot_dv_example.R --infile example_3.feather --name "plot.df.example.3"
${srun} Rscript plot_dv_example.R --infile example_4.feather --name "plot.df.example.4"
irr_Ns = [1000]
irr_ms = [150,300,600]
irr_seeds=${seeds}
irr_explained_variances=${explained_variances}
irr_coder_accuracy=[0.80]
example_5_jobs: 05_irr_indep.R irr_simulation_base.R grid_sweep.py pl_methods.R measerr_methods.R
sbatch --wait --verbose run_job.sbatch grid_sweep.py --command "Rscript 05_irr_indep.R" --arg_dict '{"N":${irr_Ns},"m":${irr_ms}, "seed":${irr_seeds}, "outfile":["example_5.feather"], "y_explained_variance":${irr_explained_variances}, "coder_accuracy":${irr_coder_accuracy}}' --outfile example_5_jobs
example_5.feather:example_5_jobs
rm -f example_5.feather
sbatch --wait --verbose --array=1-1000 run_simulation.sbatch 0 example_5_jobs
sbatch --wait --verbose --array=1001-$(shell cat example_5_jobs | wc -l) run_simulation.sbatch 1000 example_5_jobs
# sbatch --wait --verbose --array=2001-3000 run_simulation.sbatch 2000 example_5_jobs
# sbatch --wait --verbose --array=3001-4000 run_simulation.sbatch 3000 example_5_jobs
# sbatch --wait --verbose --array=2001-$(shell cat example_5_jobs | wc -l) run_simulation.sbatch 4000 example_5_jobs
# example_6_jobs: 06_irr_dv.R irr_dv_simulation_base.R grid_sweep.py pl_methods.R
# sbatch --wait --verbose run_job.sbatch grid_sweep.py --command "Rscript 06_irr_dv.R" --arg_dict '{"N":${irr_Ns},"m":${irr_ms}, "seed":${irr_seeds}, "outfile":["example_6.feather"], "y_explained_variance":${irr_explained_variances},"coder_accuracy":${irr_coder_accuracy}}' --outfile example_6_jobs
# example_6.feather:example_6_jobs
# rm -f example_6.feather
# sbatch --wait --verbose --array=1-1000 run_simulation.sbatch 0 example_6_jobs
# sbatch --wait --verbose --array=1001-2000 run_simulation.sbatch 1000 example_6_jobs
# sbatch --wait --verbose --array=2001-$(shell cat example_6_jobs | wc -l) run_simulation.sbatch 2000 example_6_jobs
remember_irr.RDS: example_5.feather plot_irr_example.R plot_irr_dv_example.R summarize_estimator.R
rm -f remember_irr.RDS
sbatch --wait --verbose run_job.sbatch Rscript plot_irr_example.R --infile example_5.feather --name "plot.df.example.5"
# sbatch --wait --verbose run_job.sbatch Rscript plot_irr_dv_example.R --infile example_6.feather --name "plot.df.example.6"
robustness_1_jobs: 02_indep_differential.R simulation_base.R grid_sweep.py
sbatch --wait --verbose run_job.sbatch grid_sweep.py --command "Rscript 02_indep_differential.R" --arg_dict '{"N":${Ns},"m":${ms}, "seed":${seeds}, "outfile":["robustness_1.feather"],"y_explained_variance":${explained_variances}, "Bzy":[-0.3],"Bxy":[0.3],"Bzx":[0.3], "outcome_formula":["y~x+z"], "proxy_formula":["w_pred~y+x"], "truth_formula":["x~1"]}' --outfile robustness_1_jobs
robustness_1.feather: robustness_1_jobs
rm -f robustness_1.feather
sbatch --wait --verbose --array=1-$(shell cat robustness_1_jobs | wc -l) run_simulation.sbatch 0 robustness_1_jobs
robustness_1_dv_jobs: simulation_base.R 04_depvar_differential.R grid_sweep.py
${srun} bash -c "source ~/.bashrc && grid_sweep.py --command 'Rscript 04_depvar_differential.R' --arg_dict \"{'N':${Ns},'m':${ms}, 'seed':${seeds}, 'outfile':['robustness_1_dv.feather'], 'y_explained_variance':${explained_variances}, 'proxy_formula':['w_pred~y']}\" --outfile robustness_1_dv_jobs"
robustness_1_dv.feather: robustness_1_dv_jobs
rm -f robustness_1_dv.feather
sbatch --wait --verbose --array=1-$(shell cat example_3_jobs | wc -l) run_simulation.sbatch 0 robustness_1_dv_jobs
remember_robustness_misspec.RDS: robustness_1.feather robustness_1_dv.feather
rm -f remember_robustness_misspec.RDS
sbatch --wait --verbose run_job.sbatch Rscript plot_example.R --infile robustness_1.feather --name "plot.df.robustness.1" --remember-file "remember_robustness_misspec.RDS"
sbatch --wait --verbose run_job.sbatch Rscript plot_dv_example.R --infile robustness_1_dv.feather --name "plot.df.robustness.1.dv" --remember-file "remember_robustness_mispec.RDS"
clean:
rm *.feather
@@ -66,7 +152,7 @@ clean:
# sbatch --wait --verbose --array=3001-6001 run_simulation.sbatch 0 example_2_B_jobs
# example_2_B_mecor_jobs:
# grid_sweep.py --command "Rscript example_2_B_mecor.R" --arg_dict '{"N":${Ns},"m":${ms}, "seed":${seeds}, "outfile":["example_2_B_mecor.feather"]}' --outfile example_2_B_mecor_jobs
# sbatch --wait --verbose run_job.sbatch grid_sweep.py --command "Rscript example_2_B_mecor.R" --arg_dict '{"N":${Ns},"m":${ms}, "seed":${seeds}, "outfile":["example_2_B_mecor.feather"]}' --outfile example_2_B_mecor_jobs
# example_2_B_mecor.feather:example_2_B_mecor.R example_2_B_mecor_jobs
# rm -f example_2_B_mecor.feather
@@ -75,3 +161,4 @@ clean:
.PHONY: supplement

BIN
simulations/example_2_B.feather
View File

Binary file not shown.

BIN
simulations/example_3.feather
View File

Binary file not shown.

7
simulations/grid_sweep.py
View File

@@ -2,8 +2,13 @@
import fire
from itertools import product
import pyRemembeR
def main(command, arg_dict, outfile):
def main(command, arg_dict, outfile, remember_file='remember_grid_sweep.RDS'):
remember = pyRemembeR.remember.Remember()
remember.set_file(remember_file)
remember[outfile] = arg_dict
remember.save_to_r()
keys = []
values = []

212
simulations/irr_dv_simulation_base.R Normal file
View File

@@ -0,0 +1,212 @@
library(matrixStats) # for numerically stable logsumexps
options(amelia.parallel="no",
amelia.ncpus=1)
library(Amelia)
source("pl_methods.R")
source("measerr_methods_2.R") ## for my more generic function.
run_simulation_depvar <- function(df, result, outcome_formula = y ~ x + z, coder_formulas = list(y.obs.0 ~ 1, y.obs.1 ~ 1), proxy_formula = w_pred ~ y.obs.1+y.obs.0+y){
(accuracy <- df[,mean(w_pred==y)])
result <- append(result, list(accuracy=accuracy))
(error.cor.z <- cor(df$x, df$w_pred - df$z))
(error.cor.x <- cor(df$x, df$w_pred - df$y))
(error.cor.y <- cor(df$y, df$y - df$w_pred))
result <- append(result, list(error.cor.x = error.cor.x,
error.cor.z = error.cor.z,
error.cor.y = error.cor.y))
model.null <- glm(y~1, data=df,family=binomial(link='logit'))
(model.true <- glm(y ~ x + z, data=df,family=binomial(link='logit')))
(lik.ratio <- exp(logLik(model.true) - logLik(model.null)))
true.ci.Bxy <- confint(model.true)['x',]
true.ci.Bzy <- confint(model.true)['z',]
result <- append(result, list(lik.ratio=lik.ratio))
result <- append(result, list(Bxy.est.true=coef(model.true)['x'],
Bzy.est.true=coef(model.true)['z'],
Bxy.ci.upper.true = true.ci.Bxy[2],
Bxy.ci.lower.true = true.ci.Bxy[1],
Bzy.ci.upper.true = true.ci.Bzy[2],
Bzy.ci.lower.true = true.ci.Bzy[1]))
(model.naive <- lm(y~w_pred+z, data=df))
naive.ci.Bxy <- confint(model.naive)['w_pred',]
naive.ci.Bzy <- confint(model.naive)['z',]
result <- append(result, list(Bxy.est.naive=coef(model.naive)['w_pred'],
Bzy.est.naive=coef(model.naive)['z'],
Bxy.ci.upper.naive = naive.ci.Bxy[2],
Bxy.ci.lower.naive = naive.ci.Bxy[1],
Bzy.ci.upper.naive = naive.ci.Bzy[2],
Bzy.ci.lower.naive = naive.ci.Bzy[1]))
loa0.feasible <- glm(y.obs.0 ~ x + z, data = df[!(is.na(y.obs.0))], family=binomial(link='logit'))
loa0.ci.Bxy <- confint(loa0.feasible)['x',]
loa0.ci.Bzy <- confint(loa0.feasible)['z',]
result <- append(result, list(Bxy.est.loa0.feasible=coef(loa0.feasible)['x'],
Bzy.est.loa0.feasible=coef(loa0.feasible)['z'],
Bxy.ci.upper.loa0.feasible = loa0.ci.Bxy[2],
Bxy.ci.lower.loa0.feasible = loa0.ci.Bxy[1],
Bzy.ci.upper.loa0.feasible = loa0.ci.Bzy[2],
Bzy.ci.lower.loa0.feasible = loa0.ci.Bzy[1]))
## df.loa0.mle <- copy(df)
## df.loa0.mle[,y:=y.obs.0]
## loa0.mle <- measerr_mle_dv(df.loa0.mle, outcome_formula=outcome_formula, proxy_formula=proxy_formula)
## fisher.info <- solve(loa0.mle$hessian)
## coef <- loa0.mle$par
## ci.upper <- coef + sqrt(diag(fisher.info)) * 1.96
## ci.lower <- coef - sqrt(diag(fisher.info)) * 1.96
## result <- append(result, list(Bxy.est.loa0.mle=coef['x'],
## Bzy.est.loa0.mle=coef['z'],
## Bxy.ci.upper.loa0.mle = ci.upper['x'],
## Bxy.ci.lower.loa0.mle = ci.lower['x'],
## Bzy.ci.upper.loa0.mle = ci.upper['z'],
## Bzy.ci.lower.loa0.mle = ci.upper['z']))
loco.feasible <- glm(y.obs.0 ~ x + z, data = df[(!is.na(y.obs.0)) & (y.obs.1 == y.obs.0)], family=binomial(link='logit'))
loco.feasible.ci.Bxy <- confint(loco.feasible)['x',]
loco.feasible.ci.Bzy <- confint(loco.feasible)['z',]
result <- append(result, list(Bxy.est.loco.feasible=coef(loco.feasible)['x'],
Bzy.est.loco.feasible=coef(loco.feasible)['z'],
Bxy.ci.upper.loco.feasible = loco.feasible.ci.Bxy[2],
Bxy.ci.lower.loco.feasible = loco.feasible.ci.Bxy[1],
Bzy.ci.upper.loco.feasible = loco.feasible.ci.Bzy[2],
Bzy.ci.lower.loco.feasible = loco.feasible.ci.Bzy[1]))
## df.double.proxy <- copy(df)
## df.double.proxy <- df.double.proxy[,y.obs:=NA]
## df.double.proxy <- df.double.proxy[,y:=NA]
## double.proxy.mle <- measerr_irr_mle_dv(df.double.proxy, outcome_formula=y~x+z, outcome_family=binomial(link='logit'), coder_formulas=list(y.obs.0 ~ y), proxy_formula=w_pred ~ y.obs.0 + y, proxy_family=binomial(link='logit'))
## print(double.proxy.mle$hessian)
## fisher.info <- solve(double.proxy.mle$hessian)
## coef <- double.proxy.mle$par
## ci.upper <- coef + sqrt(diag(fisher.info)) * 1.96
## ci.lower <- coef - sqrt(diag(fisher.info)) * 1.96
## result <- append(result, list(Bxy.est.double.proxy=coef['x'],
## Bzy.est.double.proxy=coef['z'],
## Bxy.ci.upper.double.proxy = ci.upper['x'],
## Bxy.ci.lower.double.proxy = ci.lower['x'],
## Bzy.ci.upper.double.proxy = ci.upper['z'],
## Bzy.ci.lower.double.proxy = ci.lower['z']))
df.triple.proxy <- copy(df)
df.triple.proxy <- df.triple.proxy[,y.obs:=NA]
df.triple.proxy <- df.triple.proxy[,y:=NA]
triple.proxy.mle <- measerr_irr_mle_dv(df.triple.proxy, outcome_formula=outcome_formula, outcome_family=binomial(link='logit'), coder_formulas=coder_formulas, proxy_formula=proxy_formula, proxy_family=binomial(link='logit'))
print(triple.proxy.mle$hessian)
fisher.info <- solve(triple.proxy.mle$hessian)
print(fisher.info)
coef <- triple.proxy.mle$par
ci.upper <- coef + sqrt(diag(fisher.info)) * 1.96
ci.lower <- coef - sqrt(diag(fisher.info)) * 1.96
result <- append(result, list(Bxy.est.triple.proxy=coef['x'],
Bzy.est.triple.proxy=coef['z'],
Bxy.ci.upper.triple.proxy = ci.upper['x'],
Bxy.ci.lower.triple.proxy = ci.lower['x'],
Bzy.ci.upper.triple.proxy = ci.upper['z'],
Bzy.ci.lower.triple.proxy = ci.lower['z']))
## df.loco.mle <- copy(df)
## df.loco.mle[,y.obs:=NA]
## df.loco.mle[(y.obs.0)==(y.obs.1),y.obs:=y.obs.0]
## df.loco.mle[,y.true:=y]
## df.loco.mle[,y:=y.obs]
## print(df.loco.mle[!is.na(y.obs.1),mean(y.true==y,na.rm=TRUE)])
## loco.mle <- measerr_mle_dv(df.loco.mle, outcome_formula=outcome_formula, proxy_formula=proxy_formula)
## fisher.info <- solve(loco.mle$hessian)
## coef <- loco.mle$par
## ci.upper <- coef + sqrt(diag(fisher.info)) * 1.96
## ci.lower <- coef - sqrt(diag(fisher.info)) * 1.96
## result <- append(result, list(Bxy.est.loco.mle=coef['x'],
## Bzy.est.loco.mle=coef['z'],
## Bxy.ci.upper.loco.mle = ci.upper['x'],
## Bxy.ci.lower.loco.mle = ci.lower['x'],
## Bzy.ci.upper.loco.mle = ci.upper['z'],
## Bzy.ci.lower.loco.mle = ci.lower['z']))
## my implementatoin of liklihood based correction
mod.zhang <- zhang.mle.dv(df.loco.mle)
coef <- coef(mod.zhang)
ci <- confint(mod.zhang,method='quad')
result <- append(result,
list(Bxy.est.zhang = coef['Bxy'],
Bxy.ci.upper.zhang = ci['Bxy','97.5 %'],
Bxy.ci.lower.zhang = ci['Bxy','2.5 %'],
Bzy.est.zhang = coef['Bzy'],
Bzy.ci.upper.zhang = ci['Bzy','97.5 %'],
Bzy.ci.lower.zhang = ci['Bzy','2.5 %']))
print(df.loco.mle)
# amelia says use normal distribution for binary variables.
tryCatch({
amelia.out.k <- amelia(df.loco.mle, m=200, p2s=0, idvars=c('y','ystar','w','y.obs.1','y.obs.0','y.true'))
mod.amelia.k <- zelig(y.obs~x+z, model='ls', data=amelia.out.k$imputations, cite=FALSE)
(coefse <- combine_coef_se(mod.amelia.k, messages=FALSE))
est.x.mi <- coefse['x','Estimate']
est.x.se <- coefse['x','Std.Error']
result <- append(result,
list(Bxy.est.amelia.full = est.x.mi,
Bxy.ci.upper.amelia.full = est.x.mi + 1.96 * est.x.se,
Bxy.ci.lower.amelia.full = est.x.mi - 1.96 * est.x.se
))
est.z.mi <- coefse['z','Estimate']
est.z.se <- coefse['z','Std.Error']
result <- append(result,
list(Bzy.est.amelia.full = est.z.mi,
Bzy.ci.upper.amelia.full = est.z.mi + 1.96 * est.z.se,
Bzy.ci.lower.amelia.full = est.z.mi - 1.96 * est.z.se
))
},
error = function(e){
message("An error occurred:\n",e)
result$error <- paste0(result$error,'\n', e)
})
## mle.irr <- measerr_irr_mle( df, outcome_formula = outcome_formula, rater_formula = rater_formula, proxy_formula=proxy_formula, truth_formula=truth_formula)
## fisher.info <- solve(mle.irr$hessian)
## coef <- mle.irr$par
## ci.upper <- coef + sqrt(diag(fisher.info)) * 1.96
## ci.lower <- coef - sqrt(diag(fisher.info)) * 1.96
## result <- append(result,
## list(Bxy.est.mle = coef['x'],
## Bxy.ci.upper.mle = ci.upper['x'],
## Bxy.ci.lower.mle = ci.lower['x'],
## Bzy.est.mle = coef['z'],
## Bzy.ci.upper.mle = ci.upper['z'],
## Bzy.ci.lower.mle = ci.lower['z']))
return(result)
}

222
simulations/irr_simulation_base.R Normal file
View File

@@ -0,0 +1,222 @@
library(matrixStats) # for numerically stable logsumexps
options(amelia.parallel="no",
amelia.ncpus=1)
library(Amelia)
source("measerr_methods.R")
source("pl_methods.R")
run_simulation <- function(df, result, outcome_formula = y ~ x + z, proxy_formula = w_pred ~ x, coder_formulas=list(x.obs.1~x, x.obs.0~x), truth_formula = x ~ z){
accuracy <- df[,mean(w_pred==x)]
result <- append(result, list(accuracy=accuracy))
(model.true <- lm(y ~ x + z, data=df))
true.ci.Bxy <- confint(model.true)['x',]
true.ci.Bzy <- confint(model.true)['z',]
result <- append(result, list(Bxy.est.true=coef(model.true)['x'],
Bzy.est.true=coef(model.true)['z'],
Bxy.ci.upper.true = true.ci.Bxy[2],
Bxy.ci.lower.true = true.ci.Bxy[1],
Bzy.ci.upper.true = true.ci.Bzy[2],
Bzy.ci.lower.true = true.ci.Bzy[1]))
loa0.feasible <- lm(y ~ x.obs.0 + z, data = df[!(is.na(x.obs.1))])
loa0.ci.Bxy <- confint(loa0.feasible)['x.obs.0',]
loa0.ci.Bzy <- confint(loa0.feasible)['z',]
result <- append(result, list(Bxy.est.loa0.feasible=coef(loa0.feasible)['x.obs.0'],
Bzy.est.loa0.feasible=coef(loa0.feasible)['z'],
Bxy.ci.upper.loa0.feasible = loa0.ci.Bxy[2],
Bxy.ci.lower.loa0.feasible = loa0.ci.Bxy[1],
Bzy.ci.upper.loa0.feasible = loa0.ci.Bzy[2],
Bzy.ci.lower.loa0.feasible = loa0.ci.Bzy[1]))
print("fitting loa0 model")
df.loa0.mle <- copy(df)
df.loa0.mle[,x:=x.obs.0]
loa0.mle <- measerr_mle(df.loa0.mle, outcome_formula=outcome_formula, proxy_formula=proxy_formula, truth_formula=truth_formula)
fisher.info <- solve(loa0.mle$hessian)
coef <- loa0.mle$par
ci.upper <- coef + sqrt(diag(fisher.info)) * 1.96
ci.lower <- coef - sqrt(diag(fisher.info)) * 1.96
result <- append(result, list(Bxy.est.loa0.mle=coef['x'],
Bzy.est.loa0.mle=coef['z'],
Bxy.ci.upper.loa0.mle = ci.upper['x'],
Bxy.ci.lower.loa0.mle = ci.lower['x'],
Bzy.ci.upper.loa0.mle = ci.upper['z'],
Bzy.ci.lower.loa0.mle = ci.upper['z']))
loco.feasible <- lm(y ~ x.obs.1 + z, data = df[(!is.na(x.obs.1)) & (x.obs.1 == x.obs.0)])
loco.feasible.ci.Bxy <- confint(loco.feasible)['x.obs.1',]
loco.feasible.ci.Bzy <- confint(loco.feasible)['z',]
result <- append(result, list(Bxy.est.loco.feasible=coef(loco.feasible)['x.obs.1'],
Bzy.est.loco.feasible=coef(loco.feasible)['z'],
Bxy.ci.upper.loco.feasible = loco.feasible.ci.Bxy[2],
Bxy.ci.lower.loco.feasible = loco.feasible.ci.Bxy[1],
Bzy.ci.upper.loco.feasible = loco.feasible.ci.Bzy[2],
Bzy.ci.lower.loco.feasible = loco.feasible.ci.Bzy[1]))
(model.naive <- lm(y~w_pred+z, data=df))
naive.ci.Bxy <- confint(model.naive)['w_pred',]
naive.ci.Bzy <- confint(model.naive)['z',]
result <- append(result, list(Bxy.est.naive=coef(model.naive)['w_pred'],
Bzy.est.naive=coef(model.naive)['z'],
Bxy.ci.upper.naive = naive.ci.Bxy[2],
Bxy.ci.lower.naive = naive.ci.Bxy[1],
Bzy.ci.upper.naive = naive.ci.Bzy[2],
Bzy.ci.lower.naive = naive.ci.Bzy[1]))
print("fitting loco model")
df.loco.mle <- copy(df)
df.loco.mle[,x.obs:=NA]
df.loco.mle[(x.obs.0)==(x.obs.1),x.obs:=x.obs.0]
df.loco.mle[,x.true:=x]
df.loco.mle[,x:=x.obs]
print(df.loco.mle[!is.na(x.obs.1),mean(x.true==x,na.rm=TRUE)])
loco.accuracy <- df.loco.mle[(!is.na(x.obs.1)) & (x.obs.1 == x.obs.0),mean(x.obs.1 == x.true)]
loco.mle <- measerr_mle(df.loco.mle, outcome_formula=outcome_formula, proxy_formula=proxy_formula, truth_formula=truth_formula)
fisher.info <- solve(loco.mle$hessian)
coef <- loco.mle$par
ci.upper <- coef + sqrt(diag(fisher.info)) * 1.96
ci.lower <- coef - sqrt(diag(fisher.info)) * 1.96
result <- append(result, list(loco.accuracy=loco.accuracy,
Bxy.est.loco.mle=coef['x'],
Bzy.est.loco.mle=coef['z'],
Bxy.ci.upper.loco.mle = ci.upper['x'],
Bxy.ci.lower.loco.mle = ci.lower['x'],
Bzy.ci.upper.loco.mle = ci.upper['z'],
Bzy.ci.lower.loco.mle = ci.lower['z']))
df.double.proxy.mle <- copy(df)
df.double.proxy.mle[,x.obs:=NA]
print("fitting double proxy model")
double.proxy.mle <- measerr_irr_mle(df.double.proxy.mle, outcome_formula=outcome_formula, proxy_formula=proxy_formula, coder_formulas=coder_formulas[1], truth_formula=truth_formula)
fisher.info <- solve(double.proxy.mle$hessian)
coef <- double.proxy.mle$par
ci.upper <- coef + sqrt(diag(fisher.info)) * 1.96
ci.lower <- coef - sqrt(diag(fisher.info)) * 1.96
result <- append(result, list(
Bxy.est.double.proxy=coef['x'],
Bzy.est.double.proxy=coef['z'],
Bxy.ci.upper.double.proxy = ci.upper['x'],
Bxy.ci.lower.double.proxy = ci.lower['x'],
Bzy.ci.upper.double.proxy = ci.upper['z'],
Bzy.ci.lower.double.proxy = ci.lower['z']))
df.triple.proxy.mle <- copy(df)
df.triple.proxy.mle[,x.obs:=NA]
print("fitting triple proxy model")
triple.proxy.mle <- measerr_irr_mle(df.triple.proxy.mle, outcome_formula=outcome_formula, proxy_formula=proxy_formula, coder_formulas=coder_formulas, truth_formula=truth_formula)
fisher.info <- solve(triple.proxy.mle$hessian)
coef <- triple.proxy.mle$par
ci.upper <- coef + sqrt(diag(fisher.info)) * 1.96
ci.lower <- coef - sqrt(diag(fisher.info)) * 1.96
result <- append(result, list(
Bxy.est.triple.proxy=coef['x'],
Bzy.est.triple.proxy=coef['z'],
Bxy.ci.upper.triple.proxy = ci.upper['x'],
Bxy.ci.lower.triple.proxy = ci.lower['x'],
Bzy.ci.upper.triple.proxy = ci.upper['z'],
Bzy.ci.lower.triple.proxy = ci.lower['z']))
tryCatch({
amelia.out.k <- amelia(df.loco.mle, m=200, p2s=0, idvars=c('x.true','w','x.obs.1','x.obs.0','x'))
mod.amelia.k <- zelig(y~x.obs+z, model='ls', data=amelia.out.k$imputations, cite=FALSE)
(coefse <- combine_coef_se(mod.amelia.k, messages=FALSE))
est.x.mi <- coefse['x.obs','Estimate']
est.x.se <- coefse['x.obs','Std.Error']
result <- append(result,
list(Bxy.est.amelia.full = est.x.mi,
Bxy.ci.upper.amelia.full = est.x.mi + 1.96 * est.x.se,
Bxy.ci.lower.amelia.full = est.x.mi - 1.96 * est.x.se
))
est.z.mi <- coefse['z','Estimate']
est.z.se <- coefse['z','Std.Error']
result <- append(result,
list(Bzy.est.amelia.full = est.z.mi,
Bzy.ci.upper.amelia.full = est.z.mi + 1.96 * est.z.se,
Bzy.ci.lower.amelia.full = est.z.mi - 1.96 * est.z.se
))
},
error = function(e){
message("An error occurred:\n",e)
result$error <-paste0(result$error,'\n', e)
}
)
tryCatch({
mod.zhang.lik <- zhang.mle.iv(df.loco.mle)
coef <- coef(mod.zhang.lik)
ci <- confint(mod.zhang.lik,method='quad')
result <- append(result,
list(Bxy.est.zhang = coef['Bxy'],
Bxy.ci.upper.zhang = ci['Bxy','97.5 %'],
Bxy.ci.lower.zhang = ci['Bxy','2.5 %'],
Bzy.est.zhang = coef['Bzy'],
Bzy.ci.upper.zhang = ci['Bzy','97.5 %'],
Bzy.ci.lower.zhang = ci['Bzy','2.5 %']))
},
error = function(e){
message("An error occurred:\n",e)
result$error <- paste0(result$error,'\n', e)
})
df <- df.loco.mle
N <- nrow(df)
m <- nrow(df[!is.na(x.obs)])
p <- v <- train <- rep(0,N)
M <- m
p[(M+1):(N)] <- 1
v[1:(M)] <- 1
df <- df[order(x.obs)]
y <- df[,y]
x <- df[,x.obs]
z <- df[,z]
w <- df[,w_pred]
# gmm gets pretty close
(gmm.res <- predicted_covariates(y, x, z, w, v, train, p, max_iter=100, verbose=TRUE))
result <- append(result,
list(Bxy.est.gmm = gmm.res$beta[1,1],
Bxy.ci.upper.gmm = gmm.res$confint[1,2],
Bxy.ci.lower.gmm = gmm.res$confint[1,1],
gmm.ER_pval = gmm.res$ER_pval
))
result <- append(result,
list(Bzy.est.gmm = gmm.res$beta[2,1],
Bzy.ci.upper.gmm = gmm.res$confint[2,2],
Bzy.ci.lower.gmm = gmm.res$confint[2,1]))
return(result)
}

409
simulations/measerr_methods.R
View File

@@ -1,6 +1,7 @@
library(formula.tools)
library(matrixStats)
library(optimx)
library(bbmle)
## df: dataframe to model
## outcome_formula: formula for y | x, z
## outcome_family: family for y | x, z
@@ -17,7 +18,7 @@ library(matrixStats)
## outcome_formula <- y ~ x + z; proxy_formula <- w_pred ~ y + x + z + x:z + x:y + z:y
measerr_mle_dv <- function(df, outcome_formula, outcome_family=binomial(link='logit'), proxy_formula, proxy_family=binomial(link='logit')){
measerr_mle_dv <- function(df, outcome_formula, outcome_family=binomial(link='logit'), proxy_formula, proxy_family=binomial(link='logit'),method='optim'){
nll <- function(params){
df.obs <- model.frame(outcome_formula, df)
@@ -57,7 +58,7 @@ measerr_mle_dv <- function(df, outcome_formula, outcome_family=binomial(link='lo
df.unobs.y1 <- copy(df.unobs)
df.unobs.y1[[response.var]] <- 1
df.unobs.y0 <- copy(df.unobs)
df.unobs.y0[[response.var]] <- 1
df.unobs.y0[[response.var]] <- 0
## integrate out y
outcome.model.matrix.y1 <- model.matrix(outcome_formula, df.unobs.y1)
@@ -98,23 +99,33 @@ measerr_mle_dv <- function(df, outcome_formula, outcome_family=binomial(link='lo
start <- rep(0.1,length(params))
names(start) <- params
fit <- optim(start, fn = nll, lower=lower, method='L-BFGS-B', hessian=TRUE, control=list(maxit=1e6))
if(method=='optim'){
fit <- optim(start, fn = nll, lower=lower, method='L-BFGS-B', hessian=TRUE, control=list(maxit=1e6))
} else {
quoted.names <- gsub("[\\(\\)]",'',names(start))
print(quoted.names)
text <- paste("function(", paste0(quoted.names,'=',start,collapse=','),"){params<-c(",paste0(quoted.names,collapse=','),");return(nll(params))}")
measerr_mle_nll <- eval(parse(text=text))
names(start) <- quoted.names
names(lower) <- quoted.names
fit <- mle2(minuslogl=measerr_mle_nll, start=start, lower=lower, parnames=params,control=list(maxit=1e6),method='L-BFGS-B')
}
return(fit)
}
measerr_mle <- function(df, outcome_formula, outcome_family=gaussian(), proxy_formula, proxy_family=binomial(link='logit'), truth_formula, truth_family=binomial(link='logit')){
measrr_mle_nll <- function(params){
df.obs <- model.frame(outcome_formula, df)
measerr_mle <- function(df, outcome_formula, outcome_family=gaussian(), proxy_formula, proxy_family=binomial(link='logit'), truth_formula, truth_family=binomial(link='logit'),method='optim'){
proxy.variable <- all.vars(proxy_formula)[1]
proxy.model.matrix <- model.matrix(proxy_formula, df)
response.var <- all.vars(outcome_formula)[1]
y.obs <- with(df.obs,eval(parse(text=response.var)))
outcome.model.matrix <- model.matrix(outcome_formula, df)
df.obs <- model.frame(outcome_formula, df)
response.var <- all.vars(outcome_formula)[1]
proxy.variable <- all.vars(proxy_formula)[1]
truth.variable <- all.vars(truth_formula)[1]
outcome.model.matrix <- model.matrix(outcome_formula, df)
proxy.model.matrix <- model.matrix(proxy_formula, df)
y.obs <- with(df.obs,eval(parse(text=response.var)))
measerr_mle_nll <- function(params){
param.idx <- 1
n.outcome.model.covars <- dim(outcome.model.matrix)[2]
outcome.params <- params[param.idx:n.outcome.model.covars]
@@ -124,6 +135,7 @@ measerr_mle <- function(df, outcome_formula, outcome_family=gaussian(), proxy_fo
if(outcome_family$family == "gaussian"){
sigma.y <- params[param.idx]
param.idx <- param.idx + 1
# outcome_formula likelihood using linear regression
ll.y.obs <- dnorm(y.obs, outcome.params %*% t(outcome.model.matrix),sd=sigma.y, log=TRUE)
}
@@ -135,12 +147,14 @@ measerr_mle <- function(df, outcome_formula, outcome_family=gaussian(), proxy_fo
if( (proxy_family$family=="binomial") & (proxy_family$link=='logit')){
ll.w.obs <- vector(mode='numeric',length=dim(proxy.model.matrix)[1])
# proxy_formula likelihood using logistic regression
ll.w.obs[proxy.obs==1] <- plogis(proxy.params %*% t(proxy.model.matrix[proxy.obs==1,]),log=TRUE)
ll.w.obs[proxy.obs==0] <- plogis(proxy.params %*% t(proxy.model.matrix[proxy.obs==0,]),log=TRUE, lower.tail=FALSE)
}
df.obs <- model.frame(truth_formula, df)
truth.variable <- all.vars(truth_formula)[1]
truth.obs <- with(df.obs, eval(parse(text=truth.variable)))
truth.model.matrix <- model.matrix(truth_formula,df)
n.truth.model.covars <- dim(truth.model.matrix)[2]
@@ -149,10 +163,13 @@ measerr_mle <- function(df, outcome_formula, outcome_family=gaussian(), proxy_fo
if( (truth_family$family=="binomial") & (truth_family$link=='logit')){
ll.x.obs <- vector(mode='numeric',length=dim(truth.model.matrix)[1])
# truth_formula likelihood using logistic regression
ll.x.obs[truth.obs==1] <- plogis(truth.params %*% t(truth.model.matrix[truth.obs==1,]),log=TRUE)
ll.x.obs[truth.obs==0] <- plogis(truth.params %*% t(truth.model.matrix[truth.obs==0,]),log=TRUE, lower.tail=FALSE)
}
# add the three likelihoods
ll.obs <- sum(ll.y.obs + ll.w.obs + ll.x.obs)
## likelihood for the predicted data
@@ -169,8 +186,10 @@ measerr_mle <- function(df, outcome_formula, outcome_family=gaussian(), proxy_fo
outcome.model.matrix.x0 <- model.matrix(outcome_formula, df.unobs.x0)
outcome.model.matrix.x1 <- model.matrix(outcome_formula, df.unobs.x1)
if(outcome_family$family=="gaussian"){
ll.y.x0 <- dnorm(outcome.unobs, outcome.params %*% t(outcome.model.matrix.x0), sd=sigma.y, log=TRUE)
ll.y.x1 <- dnorm(outcome.unobs, outcome.params %*% t(outcome.model.matrix.x1), sd=sigma.y, log=TRUE)
# likelihood of outcome
ll.y.x0 <- dnorm(outcome.unobs, outcome.params %*% t(outcome.model.matrix.x0), sd=sigma.y, log=TRUE)
ll.y.x1 <- dnorm(outcome.unobs, outcome.params %*% t(outcome.model.matrix.x1), sd=sigma.y, log=TRUE)
}
if( (proxy_family$family=='binomial') & (proxy_family$link=='logit')){
@@ -181,6 +200,7 @@ measerr_mle <- function(df, outcome_formula, outcome_family=gaussian(), proxy_fo
ll.w.x0 <- vector(mode='numeric', length=dim(df.unobs)[1])
ll.w.x1 <- vector(mode='numeric', length=dim(df.unobs)[1])
# likelihood of proxy
ll.w.x0[proxy.unobs==1] <- plogis(proxy.params %*% t(proxy.model.matrix.x0[proxy.unobs==1,]), log=TRUE)
ll.w.x1[proxy.unobs==1] <- plogis(proxy.params %*% t(proxy.model.matrix.x1[proxy.unobs==1,]), log=TRUE)
@@ -190,8 +210,9 @@ measerr_mle <- function(df, outcome_formula, outcome_family=gaussian(), proxy_fo
if(truth_family$link=='logit'){
truth.model.matrix <- model.matrix(truth_formula, df.unobs.x0)
ll.x.x0 <- plogis(truth.params %*% t(truth.model.matrix), log=TRUE)
ll.x.x1 <- plogis(truth.params %*% t(truth.model.matrix), log=TRUE, lower.tail=FALSE)
# likelihood of truth
ll.x.x1 <- plogis(truth.params %*% t(truth.model.matrix), log=TRUE)
ll.x.x0 <- plogis(truth.params %*% t(truth.model.matrix), log=TRUE, lower.tail=FALSE)
}
}
@@ -221,7 +242,355 @@ measerr_mle <- function(df, outcome_formula, outcome_family=gaussian(), proxy_fo
start <- rep(0.1,length(params))
names(start) <- params
fit <- optim(start, fn = measrr_mle_nll, lower=lower, method='L-BFGS-B', hessian=TRUE, control=list(maxit=1e6))
if(method=='optim'){
fit <- optim(start, fn = measerr_mle_nll, lower=lower, method='L-BFGS-B', hessian=TRUE, control=list(maxit=1e6))
} else { # method='mle2'
quoted.names <- gsub("[\\(\\)]",'',names(start))
text <- paste("function(", paste0(quoted.names,'=',start,collapse=','),"){params<-c(",paste0(quoted.names,collapse=','),");return(measerr_mle_nll(params))}")
measerr_mle_nll_mle <- eval(parse(text=text))
names(start) <- quoted.names
names(lower) <- quoted.names
fit <- mle2(minuslogl=measerr_mle_nll_mle, start=start, lower=lower, parnames=params,control=list(maxit=1e6),method='L-BFGS-B')
}
return(fit)
}
## Experimental, but probably works.
measerr_irr_mle <- function(df, outcome_formula, outcome_family=gaussian(), coder_formulas=list(x.obs.0~x, x.obs.1~x), proxy_formula, proxy_family=binomial(link='logit'), truth_formula, truth_family=binomial(link='logit'),method='optim'){
### in this scenario, the ground truth also has measurement error, but we have repeated measures for it.
# this time we never get to observe the true X
outcome.model.matrix <- model.matrix(outcome_formula, df)
proxy.model.matrix <- model.matrix(proxy_formula, df)
response.var <- all.vars(outcome_formula)[1]
proxy.var <- all.vars(proxy_formula)[1]
param.var <- all.vars(truth_formula)[1]
truth.var<- all.vars(truth_formula)[1]
y <- with(df,eval(parse(text=response.var)))
nll <- function(params){
param.idx <- 1
n.outcome.model.covars <- dim(outcome.model.matrix)[2]
outcome.params <- params[param.idx:n.outcome.model.covars]
param.idx <- param.idx + n.outcome.model.covars
if(outcome_family$family == "gaussian"){
sigma.y <- params[param.idx]
param.idx <- param.idx + 1
}
n.proxy.model.covars <- dim(proxy.model.matrix)[2]
proxy.params <- params[param.idx:(n.proxy.model.covars+param.idx-1)]
param.idx <- param.idx + n.proxy.model.covars
df.temp <- copy(df)
if((truth_family$family == "binomial")
& (truth_family$link=='logit')){
integrate.grid <- expand.grid(replicate(1 + length(coder_formulas), c(0,1), simplify=FALSE))
ll.parts <- matrix(nrow=nrow(df),ncol=nrow(integrate.grid))
for(i in 1:nrow(integrate.grid)){
# setup the dataframe for this row
row <- integrate.grid[i,]
df.temp[[param.var]] <- row[[1]]
ci <- 2
for(coder_formula in coder_formulas){
coder.var <- all.vars(coder_formula)[1]
df.temp[[coder.var]] <- row[[ci]]
ci <- ci + 1
}
outcome.model.matrix <- model.matrix(outcome_formula, df.temp)
if(outcome_family$family == "gaussian"){
ll.y <- dnorm(y, outcome.params %*% t(outcome.model.matrix), sd=sigma.y, log=TRUE)
}
if(proxy_family$family=="binomial" & (proxy_family$link=='logit')){
proxy.model.matrix <- model.matrix(proxy_formula, df.temp)
ll.w <- vector(mode='numeric', length=dim(proxy.model.matrix)[1])
proxyvar <- with(df.temp,eval(parse(text=proxy.var)))
ll.w[proxyvar==1] <- plogis(proxy.params %*% t(proxy.model.matrix[proxyvar==1,]),log=TRUE)
ll.w[proxyvar==0] <- plogis(proxy.params %*% t(proxy.model.matrix[proxyvar==0,]),log=TRUE,lower.tail=FALSE)
}
## probability of the coded variables
coder.lls <- matrix(nrow=nrow(df.temp),ncol=length(coder_formulas))
ci <- 1
for(coder_formula in coder_formulas){
coder.model.matrix <- model.matrix(coder_formula, df.temp)
n.coder.model.covars <- dim(coder.model.matrix)[2]
coder.params <- params[param.idx:(n.coder.model.covars + param.idx - 1)]
param.idx <- param.idx + n.coder.model.covars
coder.var <- all.vars(coder_formula)[1]
x.obs <- with(df.temp, eval(parse(text=coder.var)))
true.codervar <- df[[all.vars(coder_formula)[1]]]
ll.coder <- vector(mode='numeric', length=dim(coder.model.matrix)[1])
ll.coder[x.obs==1] <- plogis(coder.params %*% t(coder.model.matrix[x.obs==1,]),log=TRUE)
ll.coder[x.obs==0] <- plogis(coder.params %*% t(coder.model.matrix[x.obs==0,]),log=TRUE,lower.tail=FALSE)
# don't count when we know the observed value, unless we're accounting for observed value
ll.coder[(!is.na(true.codervar)) & (true.codervar != x.obs)] <- NA
coder.lls[,ci] <- ll.coder
ci <- ci + 1
}
truth.model.matrix <- model.matrix(truth_formula, df.temp)
n.truth.model.covars <- dim(truth.model.matrix)[2]
truth.params <- params[param.idx:(n.truth.model.covars + param.idx - 1)]
for(coder_formula in coder_formulas){
coder.model.matrix <- model.matrix(coder_formula, df.temp)
n.coder.model.covars <- dim(coder.model.matrix)[2]
param.idx <- param.idx - n.coder.model.covars
}
x <- with(df.temp, eval(parse(text=truth.var)))
ll.truth <- vector(mode='numeric', length=dim(truth.model.matrix)[1])
ll.truth[x==1] <- plogis(truth.params %*% t(truth.model.matrix[x==1,]), log=TRUE)
ll.truth[x==0] <- plogis(truth.params %*% t(truth.model.matrix[x==0,]), log=TRUE, lower.tail=FALSE)
true.truthvar <- df[[all.vars(truth_formula)[1]]]
if(!is.null(true.truthvar)){
# ll.truth[(!is.na(true.truthvar)) & (true.truthvar != truthvar)] <- -Inf
# ll.truth[(!is.na(true.truthvar)) & (true.truthvar == truthvar)] <- 0
}
ll.parts[,i] <- ll.y + ll.w + apply(coder.lls,1,sum) + ll.truth
}
lls <- rowLogSumExps(ll.parts,na.rm=TRUE)
## likelihood of observed data
target <- -1 * sum(lls)
return(target)
}
}
outcome.params <- colnames(model.matrix(outcome_formula,df))
lower <- rep(-Inf, length(outcome.params))
if(outcome_family$family=='gaussian'){
params <- c(outcome.params, 'sigma_y')
lower <- c(lower, 0.00001)
} else {
params <- outcome.params
}
proxy.params <- colnames(model.matrix(proxy_formula, df))
params <- c(params, paste0('proxy_',proxy.params))
positive.params <- paste0('proxy_',truth.var)
lower <- c(lower, rep(-Inf, length(proxy.params)))
names(lower) <- params
lower[positive.params] <- 0.01
ci <- 0
for(coder_formula in coder_formulas){
coder.params <- colnames(model.matrix(coder_formula,df))
params <- c(params, paste0('coder_',ci,coder.params))
positive.params <- paste0('coder_', ci, truth.var)
ci <- ci + 1
lower <- c(lower, rep(-Inf, length(coder.params)))
names(lower) <- params
lower[positive.params] <- 0.01
}
truth.params <- colnames(model.matrix(truth_formula, df))
params <- c(params, paste0('truth_', truth.params))
lower <- c(lower, rep(-Inf, length(truth.params)))
start <- rep(0.1,length(params))
names(start) <- params
names(lower) <- params
if(method=='optim'){
print(start)
fit <- optim(start, fn = nll, lower=lower, method='L-BFGS-B', hessian=TRUE, control=list(maxit=1e6))
} else {
quoted.names <- gsub("[\\(\\)]",'',names(start))
print(quoted.names)
text <- paste("function(", paste0(quoted.names,'=',start,collapse=','),"){params<-c(",paste0(quoted.names,collapse=','),");return(nll(params))}")
measerr_mle_nll <- eval(parse(text=text))
names(start) <- quoted.names
names(lower) <- quoted.names
fit <- mle2(minuslogl=measerr_mle_nll, start=start, lower=lower, method='L-BFGS-B',control=list(maxit=1e6))
}
return(fit)
}
## Experimental, and does not work.
measerr_irr_mle_dv <- function(df, outcome_formula, outcome_family=binomial(link='logit'), coder_formulas=list(y.obs.0~y+w_pred+y.obs.1,y.obs.1~y+w_pred+y.obs.0), proxy_formula=w_pred~y, proxy_family=binomial(link='logit'),method='optim'){
integrate.grid <- expand.grid(replicate(1 + length(coder_formulas), c(0,1), simplify=FALSE))
print(integrate.grid)
outcome.model.matrix <- model.matrix(outcome_formula, df)
n.outcome.model.covars <- dim(outcome.model.matrix)[2]
### in this scenario, the ground truth also has measurement error, but we have repeated measures for it.
# this time we never get to observe the true X
nll <- function(params){
param.idx <- 1
outcome.params <- params[param.idx:n.outcome.model.covars]
param.idx <- param.idx + n.outcome.model.covars
proxy.model.matrix <- model.matrix(proxy_formula, df)
n.proxy.model.covars <- dim(proxy.model.matrix)[2]
response.var <- all.vars(outcome_formula)[1]
if(outcome_family$family == "gaussian"){
sigma.y <- params[param.idx]
param.idx <- param.idx + 1
}
proxy.params <- params[param.idx:(n.proxy.model.covars+param.idx-1)]
param.idx <- param.idx + n.proxy.model.covars
df.temp <- copy(df)
if((outcome_family$family == "binomial")
& (outcome_family$link=='logit')){
ll.parts <- matrix(nrow=nrow(df),ncol=nrow(integrate.grid))
for(i in 1:nrow(integrate.grid)){
# setup the dataframe for this row
row <- integrate.grid[i,]
df.temp[[response.var]] <- row[[1]]
ci <- 2
for(coder_formula in coder_formulas){
codervar <- all.vars(coder_formula)[1]
df.temp[[codervar]] <- row[[ci]]
ci <- ci + 1
}
outcome.model.matrix <- model.matrix(outcome_formula, df.temp)
if(outcome_family$family == "gaussian"){
ll.y <- dnorm(df.temp[[response.var]], outcome.params %*% t(outcome.model.matrix), sd=sigma.y, log=T)
}
if(outcome_family$family == "binomial" & (outcome_family$link=='logit')){
ll.y <- vector(mode='numeric',length=nrow(df.temp))
ll.y[df.temp[[response.var]]==1] <- plogis( outcome.params %*% t(outcome.model.matrix), log=TRUE)
ll.y[df.temp[[response.var]]==0] <- plogis( outcome.params %*% t(outcome.model.matrix), log=TRUE,lower.tail=FALSE)
}
if(proxy_family$family=="binomial" & (proxy_family$link=='logit')){
proxy.model.matrix <- model.matrix(proxy_formula, df.temp)
ll.w <- vector(mode='numeric', length=dim(proxy.model.matrix)[1])
proxyvar <- with(df.temp,eval(parse(text=all.vars(proxy_formula)[1])))
ll.w[proxyvar==1] <- plogis(proxy.params %*% t(proxy.model.matrix[proxyvar==1,]),log=TRUE)
ll.w[proxyvar==0] <- plogis(proxy.params %*% t(proxy.model.matrix[proxyvar==0,]),log=TRUE,lower.tail=FALSE)
}
## probability of the coded variables
coder.lls <- matrix(nrow=nrow(df.temp),ncol=length(coder_formulas))
ci <- 1
for(coder_formula in coder_formulas){
coder.model.matrix <- model.matrix(coder_formula, df.temp)
n.coder.model.covars <- dim(coder.model.matrix)[2]
coder.params <- params[param.idx:(n.coder.model.covars + param.idx - 1)]
param.idx <- param.idx + n.coder.model.covars
codervar <- with(df.temp, eval(parse(text=all.vars(coder_formula)[1])))
true.codervar <- df[[all.vars(coder_formula)[1]]]
ll.coder <- vector(mode='numeric', length=dim(coder.model.matrix)[1])
ll.coder[codervar==1] <- plogis(coder.params %*% t(coder.model.matrix[codervar==1,]),log=TRUE)
ll.coder[codervar==0] <- plogis(coder.params %*% t(coder.model.matrix[codervar==0,]),log=TRUE,lower.tail=FALSE)
# don't count when we know the observed value, unless we're accounting for observed value
ll.coder[(!is.na(true.codervar)) & (true.codervar != codervar)] <- NA
coder.lls[,ci] <- ll.coder
ci <- ci + 1
}
for(coder_formula in coder_formulas){
coder.model.matrix <- model.matrix(coder_formula, df.temp)
n.coder.model.covars <- dim(coder.model.matrix)[2]
param.idx <- param.idx - n.coder.model.covars
}
ll.parts[,i] <- ll.y + ll.w + apply(coder.lls,1,function(x) sum(x))
}
lls <- rowLogSumExps(ll.parts,na.rm=TRUE)
## likelihood of observed data
target <- -1 * sum(lls)
print(target)
print(params)
return(target)
}
}
outcome.params <- colnames(model.matrix(outcome_formula,df))
response.var <- all.vars(outcome_formula)[1]
lower <- rep(-Inf, length(outcome.params))
if(outcome_family$family=='gaussian'){
params <- c(outcome.params, 'sigma_y')
lower <- c(lower, 0.00001)
} else {
params <- outcome.params
}
## constrain the model of the coder and proxy vars
## this is to ensure identifiability
## it is a safe assumption because the coders aren't hostile (wrong more often than right)
## so we can assume that y ~Bw, B is positive
proxy.params <- colnames(model.matrix(proxy_formula, df))
positive.params <- paste0('proxy_',response.var)
params <- c(params, paste0('proxy_',proxy.params))
lower <- c(lower, rep(-Inf, length(proxy.params)))
names(lower) <- params
lower[positive.params] <- 0.001
ci <- 0
for(coder_formula in coder_formulas){
coder.params <- colnames(model.matrix(coder_formula,df))
latent.coder.params <- coder.params %in% response.var
params <- c(params, paste0('coder_',ci,coder.params))
positive.params <- paste0('coder_',ci,response.var)
ci <- ci + 1
lower <- c(lower, rep(-Inf, length(coder.params)))
names(lower) <-params
lower[positive.params] <- 0.001
}
## init by using the "loco model"
temp.df <- copy(df)
temp.df <- temp.df[y.obs.1 == y.obs.0, y:=y.obs.1]
loco.model <- glm(outcome_formula, temp.df, family=outcome_family)
start <- rep(1,length(params))
names(start) <- params
start[names(coef(loco.model))] <- coef(loco.model)
names(lower) <- params
if(method=='optim'){
print(lower)
fit <- optim(start, fn = nll, lower=lower, method='L-BFGS-B', hessian=TRUE,control=list(maxit=1e6))
} else {
quoted.names <- gsub("[\\(\\)]",'',names(start))
print(quoted.names)
text <- paste("function(", paste0(quoted.names,'=',start,collapse=','),"){params<-c(",paste0(quoted.names,collapse=','),");return(nll(params))}")
measerr_mle_nll <- eval(parse(text=text))
names(start) <- quoted.names
names(lower) <- quoted.names
fit <- mle2(minuslogl=measerr_mle_nll, start=start, lower=lower, parnames=params,control=list(maxit=1e6),method='L-BFGS-B')
}
return(fit)
}

94
simulations/plot_dv_example.R
View File

@@ -6,53 +6,52 @@ library(filelock)
library(argparser)
parser <- arg_parser("Simulate data and fit corrected models.")
parser <- add_argument(parser, "--infile", default="", help="name of the file to read.")
parser <- add_argument(parser, "--infile", default="example_4.feather", help="name of the file to read.")
parser <- add_argument(parser, "--remember-file", default="remembr.RDS", help="name of the remember file.")
parser <- add_argument(parser, "--name", default="", help="The name to safe the data to in the remember file.")
args <- parse_args(parser)
## summarize.estimator <- function(df, suffix='naive', coefname='x'){
## part <- df[,c('N',
## 'm',
## 'Bxy',
## paste0('B',coefname,'y.est.',suffix),
## paste0('B',coefname,'y.ci.lower.',suffix),
## paste0('B',coefname,'y.ci.upper.',suffix),
## 'y_explained_variance',
## 'Bzy'
## ),
## with=FALSE]
summarize.estimator <- function(df, suffix='naive', coefname='x'){
## true.in.ci <- as.integer((part$Bxy >= part[[paste0('B',coefname,'y.ci.lower.',suffix)]]) & (part$Bxy <= part[[paste0('B',coefname,'y.ci.upper.',suffix)]]))
## zero.in.ci <- as.integer(0 >= part[[paste0('B',coefname,'y.ci.lower.',suffix)]]) & (0 <= part[[paste0('B',coefname,'y.ci.upper.',suffix)]])
## bias <- part$Bxy - part[[paste0('B',coefname,'y.est.',suffix)]]
## sign.correct <- as.integer(sign(part$Bxy) == sign(part[[paste0('B',coefname,'y.est.',suffix)]]))
part <- df[,c('N',
'm',
'Bxy',
paste0('B',coefname,'y.est.',suffix),
paste0('B',coefname,'y.ci.lower.',suffix),
paste0('B',coefname,'y.ci.upper.',suffix),
'y_explained_variance',
'Bzy',
'accuracy_imbalance_difference'
),
with=FALSE]
## part <- part[,':='(true.in.ci = true.in.ci,
## zero.in.ci = zero.in.ci,
## bias=bias,
## sign.correct =sign.correct)]
true.in.ci <- as.integer((part$Bxy >= part[[paste0('B',coefname,'y.ci.lower.',suffix)]]) & (part$Bxy <= part[[paste0('B',coefname,'y.ci.upper.',suffix)]]))
zero.in.ci <- as.integer(0 >= part[[paste0('B',coefname,'y.ci.lower.',suffix)]]) & (0 <= part[[paste0('B',coefname,'y.ci.upper.',suffix)]])
bias <- part$Bxy - part[[paste0('B',coefname,'y.est.',suffix)]]
sign.correct <- as.integer(sign(part$Bxy) == sign(part[[paste0('B',coefname,'y.est.',suffix)]]))
## part.plot <- part[, .(p.true.in.ci = mean(true.in.ci),
## mean.bias = mean(bias),
## mean.est = mean(.SD[[paste0('B',coefname,'y.est.',suffix)]]),
## var.est = var(.SD[[paste0('B',coefname,'y.est.',suffix)]]),
## est.upper.95 = quantile(.SD[[paste0('B',coefname,'y.est.',suffix)]],0.95),
## est.lower.95 = quantile(.SD[[paste0('B',coefname,'y.est.',suffix)]],0.05),
## N.sims = .N,
## p.sign.correct = mean(as.integer(sign.correct & (! zero.in.ci))),
## variable=coefname,
## method=suffix
## ),
## by=c("N","m",'Bzy','y_explained_variance')
## ]
part <- part[,':='(true.in.ci = true.in.ci,
zero.in.ci = zero.in.ci,
bias=bias,
sign.correct =sign.correct)]
part.plot <- part[, .(p.true.in.ci = mean(true.in.ci),
mean.bias = mean(bias),
mean.est = mean(.SD[[paste0('B',coefname,'y.est.',suffix)]]),
var.est = var(.SD[[paste0('B',coefname,'y.est.',suffix)]]),
est.upper.95 = quantile(.SD[[paste0('B',coefname,'y.est.',suffix)]],0.95),
est.lower.95 = quantile(.SD[[paste0('B',coefname,'y.est.',suffix)]],0.05),
N.sims = .N,
p.sign.correct = mean(as.integer(sign.correct & (! zero.in.ci))),
variable=coefname,
method=suffix
),
by=c("N","m",'Bzy','accuracy_imbalance_difference','y_explained_variance')
]
return(part.plot)
}
## return(part.plot)
## }
source("summarize_estimator.R")
build_plot_dataset <- function(df){
@@ -85,12 +84,25 @@ build_plot_dataset <- function(df){
return(plot.df)
}
df <- read_feather(args$infile)
plot.df <- build_plot_dataset(df)
change.remember.file(args$remember_file, clear=TRUE)
sims.df <- read_feather(args$infile)
sims.df[,Bzx:=NA]
sims.df[,y_explained_variance:=NA]
sims.df[,accuracy_imbalance_difference:=NA]
plot.df <- build_plot_dataset(sims.df)
remember(plot.df,args$name)
set.remember.prefix(gsub("plot.df.","",args$name))
remember(median(sims.df$cor.xz),'med.cor.xz')
remember(median(sims.df$accuracy),'med.accuracy')
remember(median(sims.df$error.cor.x),'med.error.cor.x')
remember(median(sims.df$error.cor.z),'med.error.cor.z')
remember(median(sims.df$lik.ratio),'med.lik.ratio')
## df[gmm.ER_pval<0.05]
## plot.df.test <- plot.df[,':='(method=factor(method,levels=c("Naive","Multiple imputation", "Multiple imputation (Classifier features unobserved)","Regression Calibration","2SLS+gmm","Bespoke MLE", "Feasible"),ordered=T),

122
simulations/plot_example.R
View File

@@ -5,52 +5,58 @@ library(ggplot2)
library(filelock)
library(argparser)
source("summarize_estimator.R")
parser <- arg_parser("Simulate data and fit corrected models.")
parser <- add_argument(parser, "--infile", default="", help="name of the file to read.")
parser <- add_argument(parser, "--infile", default="example_2.feather", help="name of the file to read.")
parser <- add_argument(parser, "--remember-file", default="remembr.RDS", help="name of the remember file.")
parser <- add_argument(parser, "--name", default="", help="The name to safe the data to in the remember file.")
args <- parse_args(parser)
summarize.estimator <- function(df, suffix='naive', coefname='x'){
part <- df[,c('N',
'm',
'Bxy',
paste0('B',coefname,'y.est.',suffix),
paste0('B',coefname,'y.ci.lower.',suffix),
paste0('B',coefname,'y.ci.upper.',suffix),
'y_explained_variance',
'Bzx',
'Bzy',
'accuracy_imbalance_difference'
),
with=FALSE]
true.in.ci <- as.integer((part$Bxy >= part[[paste0('B',coefname,'y.ci.lower.',suffix)]]) & (part$Bxy <= part[[paste0('B',coefname,'y.ci.upper.',suffix)]]))
zero.in.ci <- as.integer(0 >= part[[paste0('B',coefname,'y.ci.lower.',suffix)]]) & (0 <= part[[paste0('B',coefname,'y.ci.upper.',suffix)]])
bias <- part$Bxy - part[[paste0('B',coefname,'y.est.',suffix)]]
sign.correct <- as.integer(sign(part$Bxy) == sign(part[[paste0('B',coefname,'y.est.',suffix)]]))
## summarize.estimator <- function(df, suffix='naive', coefname='x'){
part <- part[,':='(true.in.ci = true.in.ci,
zero.in.ci = zero.in.ci,
bias=bias,
sign.correct =sign.correct)]
## part <- df[,c('N',
## 'm',
## 'Bxy',
## paste0('B',coefname,'y.est.',suffix),
## paste0('B',coefname,'y.ci.lower.',suffix),
## paste0('B',coefname,'y.ci.upper.',suffix),
## 'y_explained_variance',
## 'Bzx',
## 'Bzy',
## 'accuracy_imbalance_difference'
## ),
## with=FALSE]
part.plot <- part[, .(p.true.in.ci = mean(true.in.ci),
mean.bias = mean(bias),
mean.est = mean(.SD[[paste0('B',coefname,'y.est.',suffix)]]),
var.est = var(.SD[[paste0('B',coefname,'y.est.',suffix)]]),
est.upper.95 = quantile(.SD[[paste0('B',coefname,'y.est.',suffix)]],0.95,na.rm=T),
est.lower.95 = quantile(.SD[[paste0('B',coefname,'y.est.',suffix)]],0.05,na.rm=T),
N.sims = .N,
p.sign.correct = mean(as.integer(sign.correct & (! zero.in.ci))),
variable=coefname,
method=suffix
),
by=c("N","m",'y_explained_variance','Bzx', 'Bzy', 'accuracy_imbalance_difference')
]
## true.in.ci <- as.integer((part$Bxy >= part[[paste0('B',coefname,'y.ci.lower.',suffix)]]) & (part$Bxy <= part[[paste0('B',coefname,'y.ci.upper.',suffix)]]))
## zero.in.ci <- as.integer(0 >= part[[paste0('B',coefname,'y.ci.lower.',suffix)]]) & (0 <= part[[paste0('B',coefname,'y.ci.upper.',suffix)]])
## bias <- part$Bxy - part[[paste0('B',coefname,'y.est.',suffix)]]
## sign.correct <- as.integer(sign(part$Bxy) == sign(part[[paste0('B',coefname,'y.est.',suffix)]]))
return(part.plot)
}
## part <- part[,':='(true.in.ci = true.in.ci,
## zero.in.ci = zero.in.ci,
## bias=bias,
## sign.correct =sign.correct)]
## part.plot <- part[, .(p.true.in.ci = mean(true.in.ci),
## mean.bias = mean(bias),
## mean.est = mean(.SD[[paste0('B',coefname,'y.est.',suffix)]]),
## var.est = var(.SD[[paste0('B',coefname,'y.est.',suffix)]]),
## est.upper.95 = quantile(.SD[[paste0('B',coefname,'y.est.',suffix)]],0.95,na.rm=T),
## est.lower.95 = quantile(.SD[[paste0('B',coefname,'y.est.',suffix)]],0.05,na.rm=T),
## N.sims = .N,
## p.sign.correct = mean(as.integer(sign.correct & (! zero.in.ci))),
## variable=coefname,
## method=suffix
## ),
## by=c("N","m",'y_explained_variance','Bzx', 'Bzy', 'accuracy_imbalance_difference')
## ]
## return(part.plot)
## }
build_plot_dataset <- function(df){
@@ -70,13 +76,13 @@ build_plot_dataset <- function(df){
z.amelia.full <- summarize.estimator(df, 'amelia.full', 'z')
x.mecor <- summarize.estimator(df, 'mecor', 'x')
## x.mecor <- summarize.estimator(df, 'mecor', 'x')
z.mecor <- summarize.estimator(df, 'mecor', 'z')
## z.mecor <- summarize.estimator(df, 'mecor', 'z')
x.mecor <- summarize.estimator(df, 'mecor', 'x')
## x.mecor <- summarize.estimator(df, 'mecor', 'x')
z.mecor <- summarize.estimator(df, 'mecor', 'z')
## z.mecor <- summarize.estimator(df, 'mecor', 'z')
x.mle <- summarize.estimator(df, 'mle', 'x')
@@ -91,30 +97,48 @@ build_plot_dataset <- function(df){
z.gmm <- summarize.estimator(df, 'gmm', 'z')
accuracy <- df[,mean(accuracy)]
plot.df <- rbindlist(list(x.true,z.true,x.naive,z.naive,x.amelia.full,z.amelia.full,x.mecor, z.mecor, x.gmm, z.gmm, x.feasible, z.feasible,z.mle, x.mle, x.zhang, z.zhang, x.gmm, z.gmm),use.names=T)
plot.df <- rbindlist(list(x.true,z.true,x.naive,z.naive,x.amelia.full,z.amelia.full,x.gmm, z.gmm, x.feasible, z.feasible,z.mle, x.mle, x.zhang, z.zhang, x.gmm, z.gmm),use.names=T)
plot.df[,accuracy := accuracy]
plot.df <- plot.df[,":="(sd.est=sqrt(var.est)/N.sims)]
return(plot.df)
}
plot.df <- read_feather(args$infile)
sims.df <- read_feather(args$infile)
unique(sims.df[,.N,by=.(N,m)])
print(unique(sims.df$N))
# df <- df[apply(df,1,function(x) !any(is.na(x)))]
if(!('Bzx' %in% names(plot.df)))
plot.df[,Bzx:=NA]
if(!('Bzx' %in% names(sims.df)))
sims.df[,Bzx:=NA]
if(!('accuracy_imbalance_difference' %in% names(plot.df)))
plot.df[,accuracy_imbalance_difference:=NA]
if(!('accuracy_imbalance_difference' %in% names(sims.df)))
sims.df[,accuracy_imbalance_difference:=NA]
unique(plot.df[,'accuracy_imbalance_difference'])
unique(sims.df[,'accuracy_imbalance_difference'])
change.remember.file(args$remember_file, clear=TRUE)
#plot.df <- build_plot_dataset(df[accuracy_imbalance_difference==0.1][N==700])
plot.df <- build_plot_dataset(plot.df)
plot.df <- build_plot_dataset(sims.df)
remember(plot.df,args$name)
set.remember.prefix(gsub("plot.df.","",args$name))
remember(median(sims.df$cor.xz),'med.cor.xz')
remember(median(sims.df$accuracy),'med.accuracy')
remember(median(sims.df$accuracy.y0),'med.accuracy.y0')
remember(median(sims.df$accuracy.y1),'med.accuracy.y1')
remember(median(sims.df$fpr),'med.fpr')
remember(median(sims.df$fpr.y0),'med.fpr.y0')
remember(median(sims.df$fpr.y1),'med.fpr.y1')
remember(median(sims.df$fnr),'med.fnr')
remember(median(sims.df$fnr.y0),'med.fnr.y0')
remember(median(sims.df$fnr.y1),'med.fnr.y1')
remember(median(sims.df$cor.resid.w_pred),'cor.resid.w_pred')
#ggplot(df,aes(x=Bxy.est.mle)) + geom_histogram() + facet_grid(accuracy_imbalance_difference ~ Bzy)
## ## ## df[gmm.ER_pval<0.05]

73
simulations/plot_irr_dv_example.R Normal file
View File

@@ -0,0 +1,73 @@
source("RemembR/R/RemembeR.R")
library(arrow)
library(data.table)
library(ggplot2)
library(filelock)
library(argparser)
parser <- arg_parser("Simulate data and fit corrected models.")
parser <- add_argument(parser, "--infile", default="", help="name of the file to read.")
parser <- add_argument(parser, "--name", default="", help="The name to safe the data to in the remember file.")
args <- parse_args(parser)
source("summarize_estimator.R")
build_plot_dataset <- function(df){
x.true <- summarize.estimator(df, 'true','x')
z.true <- summarize.estimator(df, 'true','z')
x.naive <- summarize.estimator(df, 'naive','x')
z.naive <- summarize.estimator(df, 'naive','z')
x.loa0.feasible <- summarize.estimator(df, 'loa0.feasible','x')
z.loa0.feasible <- summarize.estimator(df,'loa0.feasible','z')
x.loa0.mle <- summarize.estimator(df, 'loa0.mle', 'x')
z.loa0.mle <- summarize.estimator(df, 'loa0.mle', 'z')
x.loco.feasible <- summarize.estimator(df, 'loco.feasible', 'x')
z.loco.feasible <- summarize.estimator(df, 'loco.feasible', 'z')
x.loco.mle <- summarize.estimator(df, 'loco.mle', 'x')
z.loco.mle <- summarize.estimator(df, 'loco.mle', 'z')
z.loco.amelia <- summarize.estimator(df, 'amelia.full', 'z')
x.loco.amelia <- summarize.estimator(df, 'amelia.full', 'x')
z.loco.zhang <- summarize.estimator(df, 'zhang', 'z')
x.loco.zhang <- summarize.estimator(df, 'zhang', 'x')
accuracy <- df[,mean(accuracy)]
plot.df <- rbindlist(list(x.true,z.true,x.loa0.feasible,z.loa0.feasible,x.naive,z.naive,x.loa0.mle,z.loa0.mle,x.loco.feasible, z.loco.feasible, z.loco.mle, x.loco.mle, x.loco.amelia, z.loco.amelia, z.loco.zhang, x.loco.zhang),use.names=T)
plot.df[,accuracy := accuracy]
plot.df <- plot.df[,":="(sd.est=sqrt(var.est)/N.sims)]
return(plot.df)
}
plot.df <- read_feather(args$infile)
print(unique(plot.df$N))
# df <- df[apply(df,1,function(x) !any(is.na(x)))]
if(!('Bzx' %in% names(plot.df)))
plot.df[,Bzx:=NA]
if(!('accuracy_imbalance_difference' %in% names(plot.df)))
plot.df[,accuracy_imbalance_difference:=NA]
unique(plot.df[,'accuracy_imbalance_difference'])
#plot.df <- build_plot_dataset(df[accuracy_imbalance_difference==0.1][N==700])
plot.df <- build_plot_dataset(plot.df)
change.remember.file("remember_irr.RDS",clear=TRUE)
remember(plot.df,args$name)

152
simulations/plot_irr_example.R Normal file
View File

@@ -0,0 +1,152 @@
source("RemembR/R/RemembeR.R")
library(arrow)
library(data.table)
library(ggplot2)
library(filelock)
library(argparser)
parser <- arg_parser("Simulate data and fit corrected models.")
parser <- add_argument(parser, "--infile", default="", help="name of the file to read.")
parser <- add_argument(parser, "--name", default="", help="The name to safe the data to in the remember file.")
args <- parse_args(parser)
source("summarize_estimator.R")
build_plot_dataset <- function(df){
x.true <- summarize.estimator(df, 'true','x')
z.true <- summarize.estimator(df, 'true','z')
x.naive <- summarize.estimator(df, 'naive','x')
z.naive <- summarize.estimator(df, 'naive','z')
x.loa0.feasible <- summarize.estimator(df, 'loa0.feasible','x')
z.loa0.feasible <- summarize.estimator(df,'loa0.feasible','z')
x.loa0.mle <- summarize.estimator(df, 'loa0.mle', 'x')
z.loa0.mle <- summarize.estimator(df, 'loa0.mle', 'z')
x.loco.feasible <- summarize.estimator(df, 'loco.feasible', 'x')
z.loco.feasible <- summarize.estimator(df, 'loco.feasible', 'z')
x.loco.mle <- summarize.estimator(df, 'loco.mle', 'x')
z.loco.mle <- summarize.estimator(df, 'loco.mle', 'z')
x.loco.mle <- summarize.estimator(df, 'loco.mle', 'x')
z.loco.amelia <- summarize.estimator(df, 'amelia.full', 'z')
x.loco.amelia <- summarize.estimator(df, 'amelia.full', 'x')
z.loco.zhang <- summarize.estimator(df, 'zhang', 'z')
x.loco.zhang <- summarize.estimator(df, 'zhang', 'x')
z.loco.gmm <- summarize.estimator(df, 'gmm', 'z')
x.loco.gmm <- summarize.estimator(df, 'gmm', 'x')
## x.mle <- summarize.estimator(df, 'mle', 'x')
## z.mle <- summarize.estimator(df, 'mle', 'z')
accuracy <- df[,mean(accuracy)]
plot.df <- rbindlist(list(x.true,z.true,x.loa0.feasible,z.loa0.feasible,x.loa0.mle,z.loa0.mle,x.loco.feasible, z.loco.feasible, x.loco.mle, z.loco.mle, x.loco.amelia, z.loco.amelia,x.loco.zhang, z.loco.zhang,x.loco.gmm, z.loco.gmm,x.naive,z.naive),use.names=T)
plot.df[,accuracy := accuracy]
plot.df <- plot.df[,":="(sd.est=sqrt(var.est)/N.sims)]
return(plot.df)
}
sims.df <- read_feather(args$infile)
print(unique(sims.df$N))
# df <- df[apply(df,1,function(x) !any(is.na(x)))]
if(!('Bzx' %in% names(sims.df)))
sims.df[,Bzx:=NA]
if(!('accuracy_imbalance_difference' %in% names(sims.df)))
sims.df[,accuracy_imbalance_difference:=NA]
unique(sims.df[,'accuracy_imbalance_difference'])
#plot.df <- build_plot_dataset(df[accuracy_imbalance_difference==0.1][N==700])
plot.df <- build_plot_dataset(sims.df)
change.remember.file("remember_irr.RDS",clear=TRUE)
remember(plot.df,args$name)
set.remember.prefix(gsub("plot.df.","",args$name))
remember(median(sims.df$loco.accuracy),'med.loco.acc')
#ggplot(df,aes(x=Bxy.est.mle)) + geom_histogram() + facet_grid(accuracy_imbalance_difference ~ Bzy)
## ## ## df[gmm.ER_pval<0.05]
## plot.df.test <- plot.df[,':='(method=factor(method,levels=c("Naive","Multiple imputation", "Multiple imputation (Classifier features unobserved)","Regression Calibration","2SLS+gmm","Bespoke MLE", "Feasible"),ordered=T),
## N=factor(N),
## m=factor(m))]
## plot.df.test <- plot.df.test[(variable=='x') & (method!="Multiple imputation (Classifier features unobserved)")]
## p <- ggplot(plot.df.test, aes(y=mean.est, ymax=mean.est + var.est/2, ymin=mean.est-var.est/2, x=method))
## p <- p + geom_hline(data=plot.df.test, mapping=aes(yintercept=0.1),linetype=2)
## p <- p + geom_pointrange() + facet_grid(N~m,as.table=F,scales='free') + scale_x_discrete(labels=label_wrap_gen(4))
## print(p)
## plot.df.test <- plot.df[,':='(method=factor(method,levels=c("Naive","Multiple imputation", "Multiple imputation (Classifier features unobserved)","Regression Calibration","2SLS+gmm","Bespoke MLE", "Feasible"),ordered=T),
## N=factor(N),
## m=factor(m))]
## plot.df.test <- plot.df.test[(variable=='z') & (method!="Multiple imputation (Classifier features unobserved)")]
## p <- ggplot(plot.df.test, aes(y=mean.est, ymax=mean.est + var.est/2, ymin=mean.est-var.est/2, x=method))
## p <- p + geom_hline(data=plot.df.test, mapping=aes(yintercept=-0.1),linetype=2)
## p <- p + geom_pointrange() + facet_grid(m~N,as.table=F,scales='free') + scale_x_discrete(labels=label_wrap_gen(4))
## print(p)
## x.mle <- df[,.(N,m,Bxy.est.mle,Bxy.ci.lower.mle, Bxy.ci.upper.mle, y_explained_variance, Bzx, Bzy, accuracy_imbalance_difference)]
## x.mle.plot <- x.mle[,.(mean.est = mean(Bxy.est.mle),
## var.est = var(Bxy.est.mle),
## N.sims = .N,
## variable='z',
## method='Bespoke MLE'
## ),
## by=c("N","m",'y_explained_variance', 'Bzx', 'Bzy','accuracy_imbalance_difference')]
## z.mle <- df[,.(N,m,Bzy.est.mle,Bzy.ci.lower.mle, Bzy.ci.upper.mle, y_explained_variance, Bzx, Bzy, accuracy_imbalance_difference)]
## z.mle.plot <- z.mle[,.(mean.est = mean(Bzy.est.mle),
## var.est = var(Bzy.est.mle),
## N.sims = .N,
## variable='z',
## method='Bespoke MLE'
## ),
## by=c("N","m",'y_explained_variance','Bzx')]
## plot.df <- z.mle.plot
## plot.df.test <- plot.df[,':='(method=factor(method,levels=c("Naive","Multiple imputation", "Multiple imputation (Classifier features unobserved)","Regression Calibration","2SLS+gmm","Bespoke MLE", "Feasible"),ordered=T),
## N=factor(N),
## m=factor(m))]
## plot.df.test <- plot.df.test[(variable=='z') & (m != 1000) & (m!=500) & (method!="Multiple imputation (Classifier features unobserved)")]
## p <- ggplot(plot.df.test, aes(y=mean.est, ymax=mean.est + var.est/2, ymin=mean.est-var.est/2, x=method))
## p <- p + geom_hline(aes(yintercept=0.2),linetype=2)
## p <- p + geom_pointrange() + facet_grid(m~Bzx, Bzy,as.table=F) + scale_x_discrete(labels=label_wrap_gen(4))
## print(p)
## ## ggplot(plot.df[variable=='x'], aes(y=mean.est, ymax=mean.est + var.est/2, ymin=mean.est-var.est/2, x=method)) + geom_pointrange() + facet_grid(-m~N) + scale_x_discrete(labels=label_wrap_gen(10))
## ## ggplot(plot.df,aes(y=N,x=m,color=p.sign.correct)) + geom_point() + facet_grid(variable ~ method) + scale_color_viridis_c(option='D') + theme_minimal() + xlab("Number of gold standard labels") + ylab("Total sample size")
## ## ggplot(plot.df,aes(y=N,x=m,color=abs(mean.bias))) + geom_point() + facet_grid(variable ~ method) + scale_color_viridis_c(option='D') + theme_minimal() + xlab("Number of gold standard labels") + ylab("Total sample size")

12
simulations/run_simulation.sbatch
View File

@@ -1,21 +1,23 @@
#!/bin/bash
#SBATCH --job-name="simulate measurement error models"
## Allocation Definition
#SBATCH --account=comdata
#SBATCH --partition=compute-bigmem
#SBATCH --account=comdata-ckpt
#SBATCH --partition=ckpt
## Resources
#SBATCH --nodes=1
## Walltime (12 hours)
#SBATCH --time=1:00:00
## Walltime (4 hours)
#SBATCH --time=4:00:00
## Memory per node
#SBATCH --mem=8G
#SBATCH --mem=4G
#SBATCH --cpus-per-task=1
#SBATCH --ntasks-per-node=1
#SBATCH --chdir /gscratch/comdata/users/nathante/ml_measurement_error_public/simulations
#SBATCH --output=simulation_jobs/%A_%a.out
#SBATCH --error=simulation_jobs/%A_%a.err
source ~/.bashrc
TASK_NUM=$(($SLURM_ARRAY_TASK_ID + $1))
TASK_CALL=$(sed -n ${TASK_NUM}p $2)
echo ${TASK_CALL}
${TASK_CALL}

233
simulations/simulation_base.R
View File

@@ -7,6 +7,7 @@ library(Zelig)
library(bbmle)
library(matrixStats) # for numerically stable logsumexps
source("pl_methods.R")
source("measerr_methods.R") ## for my more generic function.
## This uses the pseudolikelihood approach from Carroll page 349.
@@ -37,84 +38,6 @@ my.pseudo.mle <- function(df){
}
## model from Zhang's arxiv paper, with predictions for y
## Zhang got this model from Hausman 1998
### I think this is actually eqivalent to the pseudo.mle method
zhang.mle.iv <- function(df){
nll <- function(B0=0, Bxy=0, Bzy=0, sigma_y=0.1, ppv=0.9, npv=0.9){
df.obs <- df[!is.na(x.obs)]
df.unobs <- df[is.na(x.obs)]
## fpr = 1 - TNR
### Problem: accounting for uncertainty in ppv / npv
ll.w1x1.obs <- with(df.obs[(w_pred==1)], dbinom(x.obs,size=1,prob=ppv,log=T))
ll.w0x0.obs <- with(df.obs[(w_pred==0)], dbinom(1-x.obs,size=1,prob=npv,log=T))
## fnr = 1 - TPR
ll.y.obs <- with(df.obs, dnorm(y, B0 + Bxy * x + Bzy * z, sd=sigma_y,log=T))
ll <- sum(ll.y.obs)
ll <- ll + sum(ll.w1x1.obs) + sum(ll.w0x0.obs)
# unobserved case; integrate out x
ll.x.1 <- with(df.unobs, dnorm(y, B0 + Bxy + Bzy * z, sd = sigma_y, log=T))
ll.x.0 <- with(df.unobs, dnorm(y, B0 + Bzy * z, sd = sigma_y,log=T))
## case x == 1
lls.x.1 <- colLogSumExps(rbind(log(ppv) + ll.x.1, log(1-ppv) + ll.x.0))
## case x == 0
lls.x.0 <- colLogSumExps(rbind(log(1-npv) + ll.x.1, log(npv) + ll.x.0))
lls <- colLogSumExps(rbind(lls.x.1, lls.x.0))
ll <- ll + sum(lls)
return(-ll)
}
mlefit <- mle2(minuslogl = nll, control=list(maxit=1e6), lower=list(sigma_y=0.0001, B0=-Inf, Bxy=-Inf, Bzy=-Inf,ppv=0.00001, npv=0.00001),
upper=list(sigma_y=Inf, B0=Inf, Bxy=Inf, Bzy=Inf, ppv=0.99999,npv=0.99999),method='L-BFGS-B')
return(mlefit)
}
## this is equivalent to the pseudo-liklihood model from Carolla
zhang.mle.dv <- function(df){
nll <- function(B0=0, Bxy=0, Bzy=0, ppv=0.9, npv=0.9){
df.obs <- df[!is.na(y.obs)]
## fpr = 1 - TNR
ll.w0y0 <- with(df.obs[y.obs==0],dbinom(1-w_pred,1,npv,log=TRUE))
ll.w1y1 <- with(df.obs[y.obs==1],dbinom(w_pred,1,ppv,log=TRUE))
# observed case
ll.y.obs <- vector(mode='numeric', length=nrow(df.obs))
ll.y.obs[df.obs$y.obs==1] <- with(df.obs[y.obs==1], plogis(B0 + Bxy * x + Bzy * z,log=T))
ll.y.obs[df.obs$y.obs==0] <- with(df.obs[y.obs==0], plogis(B0 + Bxy * x + Bzy * z,log=T,lower.tail=FALSE))
ll <- sum(ll.y.obs) + sum(ll.w0y0) + sum(ll.w1y1)
# unobserved case; integrate out y
## case y = 1
ll.y.1 <- vector(mode='numeric', length=nrow(df))
pi.y.1 <- with(df,plogis(B0 + Bxy * x + Bzy*z, log=T))
## P(w=1| y=1)P(y=1) + P(w=0|y=1)P(y=1) = P(w=1,y=1) + P(w=0,y=1)
lls.y.1 <- colLogSumExps(rbind(log(ppv) + pi.y.1, log(1-ppv) + pi.y.1))
## case y = 0
ll.y.0 <- vector(mode='numeric', length=nrow(df))
pi.y.0 <- with(df,plogis(B0 + Bxy * x + Bzy*z, log=T,lower.tail=FALSE))
## P(w=1 | y=0)P(y=0) + P(w=0|y=0)P(y=0) = P(w=1,y=0) + P(w=0,y=0)
lls.y.0 <- colLogSumExps(rbind(log(npv) + pi.y.0, log(1-npv) + pi.y.0))
lls <- colLogSumExps(rbind(lls.y.1, lls.y.0))
ll <- ll + sum(lls)
return(-ll)
}
mlefit <- mle2(minuslogl = nll, control=list(maxit=1e6),method='L-BFGS-B',lower=list(B0=-Inf, Bxy=-Inf, Bzy=-Inf, ppv=0.001,npv=0.001),
upper=list(B0=Inf, Bxy=Inf, Bzy=Inf,ppv=0.999,npv=0.999))
return(mlefit)
}
## This uses the likelihood approach from Carroll page 353.
## assumes that we have a good measurement error model
my.mle <- function(df){
@@ -168,13 +91,25 @@ my.mle <- function(df){
run_simulation_depvar <- function(df, result, outcome_formula=y~x+z, proxy_formula=w_pred~y){
accuracy <- df[,mean(w_pred==y)]
(accuracy <- df[,mean(w_pred==y)])
result <- append(result, list(accuracy=accuracy))
(error.cor.z <- cor(df$z, df$y - df$w_pred))
(error.cor.x <- cor(df$x, df$y - df$w_pred))
(error.cor.y <- cor(df$y, df$y - df$w_pred))
result <- append(result, list(error.cor.x = error.cor.x,
error.cor.z = error.cor.z,
error.cor.y = error.cor.y))
model.null <- glm(y~1, data=df,family=binomial(link='logit'))
(model.true <- glm(y ~ x + z, data=df,family=binomial(link='logit')))
(lik.ratio <- exp(logLik(model.true) - logLik(model.null)))
true.ci.Bxy <- confint(model.true)['x',]
true.ci.Bzy <- confint(model.true)['z',]
result <- append(result, list(cor.xz=cor(df$x,df$z)))
result <- append(result, list(lik.ratio=lik.ratio))
result <- append(result, list(Bxy.est.true=coef(model.true)['x'],
Bzy.est.true=coef(model.true)['z'],
Bxy.ci.upper.true = true.ci.Bxy[2],
@@ -211,7 +146,7 @@ run_simulation_depvar <- function(df, result, outcome_formula=y~x+z, proxy_formu
naivecont.ci.Bxy <- confint(model.naive.cont)['x',]
naivecont.ci.Bzy <- confint(model.naive.cont)['z',]
## my implementatoin of liklihood based correction
## my implementation of liklihood based correction
temp.df <- copy(df)
temp.df[,y:=y.obs]
@@ -243,34 +178,28 @@ run_simulation_depvar <- function(df, result, outcome_formula=y~x+z, proxy_formu
Bzy.ci.lower.zhang = ci['Bzy','2.5 %']))
# amelia says use normal distribution for binary variables.
tryCatch({
amelia.out.k <- amelia(df, m=200, p2s=0, idvars=c('y','ystar','w'))
mod.amelia.k <- zelig(y.obs~x+z, model='ls', data=amelia.out.k$imputations, cite=FALSE)
(coefse <- combine_coef_se(mod.amelia.k, messages=FALSE))
est.x.mi <- coefse['x','Estimate']
est.x.se <- coefse['x','Std.Error']
result <- append(result,
list(Bxy.est.amelia.full = est.x.mi,
Bxy.ci.upper.amelia.full = est.x.mi + 1.96 * est.x.se,
Bxy.ci.lower.amelia.full = est.x.mi - 1.96 * est.x.se
))
est.z.mi <- coefse['z','Estimate']
est.z.se <- coefse['z','Std.Error']
amelia.out.k <- amelia(df, m=200, p2s=0, idvars=c('y','ystar','w'))
mod.amelia.k <- zelig(y.obs~x+z, model='ls', data=amelia.out.k$imputations, cite=FALSE)
(coefse <- combine_coef_se(mod.amelia.k, messages=FALSE))
est.x.mi <- coefse['x','Estimate']
est.x.se <- coefse['x','Std.Error']
result <- append(result,
list(Bxy.est.amelia.full = est.x.mi,
Bxy.ci.upper.amelia.full = est.x.mi + 1.96 * est.x.se,
Bxy.ci.lower.amelia.full = est.x.mi - 1.96 * est.x.se
))
result <- append(result,
list(Bzy.est.amelia.full = est.z.mi,
Bzy.ci.upper.amelia.full = est.z.mi + 1.96 * est.z.se,
Bzy.ci.lower.amelia.full = est.z.mi - 1.96 * est.z.se
))
},
error = function(e){
message("An error occurred:\n",e)
result$error <- paste0(result$error,'\n', e)
})
est.z.mi <- coefse['z','Estimate']
est.z.se <- coefse['z','Std.Error']
result <- append(result,
list(Bzy.est.amelia.full = est.z.mi,
Bzy.ci.upper.amelia.full = est.z.mi + 1.96 * est.z.se,
Bzy.ci.lower.amelia.full = est.z.mi - 1.96 * est.z.se
))
return(result)
@@ -278,11 +207,36 @@ run_simulation_depvar <- function(df, result, outcome_formula=y~x+z, proxy_formu
## outcome_formula, proxy_formula, and truth_formula are passed to measerr_mle
run_simulation <- function(df, result, outcome_formula=y~x+z, proxy_formula=w_pred~x, truth_formula=x~z){
run_simulation <- function(df, result, outcome_formula=y~x+z, proxy_formula=NULL, truth_formula=NULL){
accuracy <- df[,mean(w_pred==x)]
result <- append(result, list(accuracy=accuracy))
accuracy.y0 <- df[y<=0,mean(w_pred==x)]
accuracy.y1 <- df[y>=0,mean(w_pred==x)]
cor.y.xi <- cor(df$x - df$w_pred, df$y)
fnr <- df[w_pred==0,mean(w_pred!=x)]
fnr.y0 <- df[(w_pred==0) & (y<=0),mean(w_pred!=x)]
fnr.y1 <- df[(w_pred==0) & (y>=0),mean(w_pred!=x)]
fpr <- df[w_pred==1,mean(w_pred!=x)]
fpr.y0 <- df[(w_pred==1) & (y<=0),mean(w_pred!=x)]
fpr.y1 <- df[(w_pred==1) & (y>=0),mean(w_pred!=x)]
cor.resid.w_pred <- cor(resid(lm(y~x+z,df)),df$w_pred)
result <- append(result, list(accuracy=accuracy,
accuracy.y0=accuracy.y0,
accuracy.y1=accuracy.y1,
cor.y.xi=cor.y.xi,
fnr=fnr,
fnr.y0=fnr.y0,
fnr.y1=fnr.y1,
fpr=fpr,
fpr.y0=fpr.y0,
fpr.y1=fpr.y1,
cor.resid.w_pred=cor.resid.w_pred
))
result <- append(result, list(cor.xz=cor(df$x,df$z)))
(model.true <- lm(y ~ x + z, data=df))
true.ci.Bxy <- confint(model.true)['x',]
true.ci.Bzy <- confint(model.true)['z',]
@@ -319,8 +273,8 @@ run_simulation <- function(df, result, outcome_formula=y~x+z, proxy_formula=w_p
Bzy.ci.lower.naive = naive.ci.Bzy[1]))
tryCatch({
amelia.out.k <- amelia(df, m=200, p2s=0, idvars=c('x','w_pred'))
amelia.out.k <- amelia(df, m=200, p2s=0, idvars=c('x','w'))
mod.amelia.k <- zelig(y~x.obs+z, model='ls', data=amelia.out.k$imputations, cite=FALSE)
(coefse <- combine_coef_se(mod.amelia.k, messages=FALSE))
@@ -341,14 +295,7 @@ run_simulation <- function(df, result, outcome_formula=y~x+z, proxy_formula=w_p
Bzy.ci.lower.amelia.full = est.z.mi - 1.96 * est.z.se
))
},
error = function(e){
message("An error occurred:\n",e)
result$error <-paste0(result$error,'\n', e)
}
)
tryCatch({
temp.df <- copy(df)
temp.df <- temp.df[,x:=x.obs]
mod.caroll.lik <- measerr_mle(temp.df, outcome_formula=outcome_formula, proxy_formula=proxy_formula, truth_formula=truth_formula)
@@ -365,14 +312,6 @@ run_simulation <- function(df, result, outcome_formula=y~x+z, proxy_formula=w_p
Bzy.est.mle = coef['z'],
Bzy.ci.upper.mle = ci.upper['z'],
Bzy.ci.lower.mle = ci.lower['z']))
},
error = function(e){
message("An error occurred:\n",e)
result$error <- paste0(result$error,'\n', e)
})
tryCatch({
mod.zhang.lik <- zhang.mle.iv(df)
coef <- coef(mod.zhang.lik)
@@ -384,12 +323,6 @@ run_simulation <- function(df, result, outcome_formula=y~x+z, proxy_formula=w_p
Bzy.est.zhang = coef['Bzy'],
Bzy.ci.upper.zhang = ci['Bzy','97.5 %'],
Bzy.ci.lower.zhang = ci['Bzy','2.5 %']))
},
error = function(e){
message("An error occurred:\n",e)
result$error <- paste0(result$error,'\n', e)
})
## What if we can't observe k -- most realistic scenario. We can't include all the ML features in a model.
## amelia.out.nok <- amelia(df, m=200, p2s=0, idvars=c("x","w_pred"), noms=noms)
@@ -440,29 +373,29 @@ run_simulation <- function(df, result, outcome_formula=y~x+z, proxy_formula=w_p
Bzy.ci.lower.gmm = gmm.res$confint[2,1]))
tryCatch({
mod.calibrated.mle <- mecor(y ~ MeasError(w_pred, reference = x.obs) + z, df, B=400, method='efficient')
(mod.calibrated.mle)
(mecor.ci <- summary(mod.calibrated.mle)$c$ci['x.obs',])
result <- append(result, list(
Bxy.est.mecor = mecor.ci['Estimate'],
Bxy.ci.upper.mecor = mecor.ci['UCI'],
Bxy.ci.lower.mecor = mecor.ci['LCI'])
)
## tryCatch({
## mod.calibrated.mle <- mecor(y ~ MeasError(w_pred, reference = x.obs) + z, df, B=400, method='efficient')
## (mod.calibrated.mle)
## (mecor.ci <- summary(mod.calibrated.mle)$c$ci['x.obs',])
## result <- append(result, list(
## Bxy.est.mecor = mecor.ci['Estimate'],
## Bxy.ci.upper.mecor = mecor.ci['UCI'],
## Bxy.ci.lower.mecor = mecor.ci['LCI'])
## )
(mecor.ci <- summary(mod.calibrated.mle)$c$ci['z',])
## (mecor.ci <- summary(mod.calibrated.mle)$c$ci['z',])
result <- append(result, list(
Bzy.est.mecor = mecor.ci['Estimate'],
Bzy.ci.upper.mecor = mecor.ci['UCI'],
Bzy.ci.lower.mecor = mecor.ci['LCI'])
)
},
error = function(e){
message("An error occurred:\n",e)
result$error <- paste0(result$error, '\n', e)
}
)
## result <- append(result, list(
## Bzy.est.mecor = mecor.ci['Estimate'],
## Bzy.ci.upper.mecor = mecor.ci['UCI'],
## Bzy.ci.lower.mecor = mecor.ci['LCI'])
## )
## },
## error = function(e){
## message("An error occurred:\n",e)
## result$error <- paste0(result$error, '\n', e)
## }
## )
## clean up memory
## rm(list=c("df","y","x","g","w","v","train","p","amelia.out.k","amelia.out.nok", "mod.calibrated.mle","gmm.res","mod.amelia.k","mod.amelia.nok", "model.true","model.naive","model.feasible"))

50
simulations/summarize_estimator.R Normal file
View File

@@ -0,0 +1,50 @@
summarize.estimator <- function(df, suffix='naive', coefname='x'){
part <- df[,c('N',
'm',
'Bxy',
paste0('B',coefname,'y.est.',suffix),
paste0('B',coefname,'y.ci.lower.',suffix),
paste0('B',coefname,'y.ci.upper.',suffix),
'y_explained_variance',
'Bzx',
'Bzy',
'accuracy_imbalance_difference'
),
with=FALSE]
true.in.ci <- as.integer((part$Bxy >= part[[paste0('B',coefname,'y.ci.lower.',suffix)]]) & (part$Bxy <= part[[paste0('B',coefname,'y.ci.upper.',suffix)]]))
zero.in.ci <- as.integer(0 >= part[[paste0('B',coefname,'y.ci.lower.',suffix)]]) & (0 <= part[[paste0('B',coefname,'y.ci.upper.',suffix)]])
bias <- part[[paste0('B',coefname,'y')]] - part[[paste0('B',coefname,'y.est.',suffix)]]
sign.correct <- as.integer(sign(part$Bxy) == sign(part[[paste0('B',coefname,'y.est.',suffix)]]))
part <- part[,':='(true.in.ci = true.in.ci,
zero.in.ci = zero.in.ci,
bias=bias,
sign.correct =sign.correct)]
part.plot <- part[, .(p.true.in.ci = mean(true.in.ci),
mean.bias = mean(bias),
mean.est = mean(.SD[[paste0('B',coefname,'y.est.',suffix)]]),
var.est = var(.SD[[paste0('B',coefname,'y.est.',suffix)]]),
est.upper.95 = quantile(.SD[[paste0('B',coefname,'y.est.',suffix)]],0.975,na.rm=T),
est.lower.95 = quantile(.SD[[paste0('B',coefname,'y.est.',suffix)]],0.025,na.rm=T),
mean.ci.upper = mean(.SD[[paste0('B',coefname,'y.ci.upper.',suffix)]],na.rm=T),
mean.ci.lower = mean(.SD[[paste0('B',coefname,'y.ci.lower.',suffix)]],na.rm=T),
ci.upper.975 = quantile(.SD[[paste0('B',coefname,'y.ci.upper.',suffix)]],0.975,na.rm=T),
ci.upper.025 = quantile(.SD[[paste0('B',coefname,'y.ci.upper.',suffix)]],0.025,na.rm=T),
ci.lower.975 = quantile(.SD[[paste0('B',coefname,'y.ci.lower.',suffix)]],0.975,na.rm=T),
ci.lower.025 = quantile(.SD[[paste0('B',coefname,'y.ci.lower.',suffix)]],0.025,na.rm=T),
N.ci.is.NA = sum(is.na(.SD[[paste0('B',coefname,'y.ci.lower.',suffix)]])),
N.sims = .N,
p.sign.correct = mean(as.integer(sign.correct & (! zero.in.ci))),
variable=coefname,
method=suffix
),
by=c("N","m",'y_explained_variance','Bzx', 'Bzy', 'accuracy_imbalance_difference')
]
return(part.plot)
}