####### RQ1 solutions R version 3.2.2 (2015-08-14) -- "Fire Safety" # mercury from fish > library(sensitivitymw) ## mv not needed here, but... > install.packages("sensitivitymv") > library(sensitivitymv) ################################# part a > ## mercury in fish > ## vignette sec 3 matched sets with multiple controls # two controls for every treated (see question text) > data(mercury) > head(mercury) Treated Zero One 1 4.60 0.23 0.42 2 0.85 0.76 0.34 3 0.59 0.23 0.23 4 1.39 0.23 0.85 5 17.09 0.23 0.75 6 2.21 0.83 1.34 > fivenum(mercury$Treated) [1] 0.23 1.30 2.36 4.28 38.00 > dim(mercury) [1] 397 3 > length(unlist(mercury[,2:3])) #combined controls both cols [1] 794 > # so median 2.36 in treated vs .52 controls of mercury in blood, you care > fivenum(unlist(mercury[,2:3])) #combined controls both cols Zero1 Zero102 Zero284 One150 One13 0.23 0.33 0.52 0.88 14.60 > # matches page 8 vign > # try senmw two ways: 't' and now 'h' Huber, see snemw options in sensitivitymw manual > senmw(mercury, gamma = 1, method = "t")$pval [1] 0 > senmw(mercury, gamma = 1, method = "h")$pval [1] 0 ## big significant effect from data analysis (equiv to gamma = 1) ## how big does gamma need to be to negate this highly significant effect? # remember gamma = 3 is a pretty big distortion of propensity matching # here it takes gamma in the teens to have any consequence, wow # note below are one-sided p-vals, double to be careful # gamma <10 has no real consequence, either t or h methods > senmw(mercury, gamma = 4, method = "t")$pval [1] 4.628853e-12 > senmw(mercury, gamma = 4, method = "h")$pval [1] 0 > senmw(mercury, gamma = 6, method = "t")$pval [1] 1.855396e-07 > senmw(mercury, gamma = 6, method = "h")$pval [1] 1.092164e-10 > senmw(mercury, gamma = 8, method = "t")$pval [1] 3.413761e-05 > senmw(mercury, gamma = 8, method = "h")$pval [1] 1.092262e-06 > senmw(mercury, gamma = 10, method = "t")$pval [1] 0.0007221081 > senmw(mercury, gamma = 10, method = "h")$pval [1] 0.0001911545 > senmw(mercury, gamma = 11, method = "t")$pval [1] 0.002131966 > senmw(mercury, gamma = 11, method = "h")$pval [1] 0.001105712 > senmw(mercury, gamma = 12, method = "t")$pval [1] 0.005167374 > senmw(mercury, gamma = 12, method = "h")$pval [1] 0.004432061 > senmw(mercury, gamma = 13, method = "t")$pval [1] 0.01076347 > senmw(mercury, gamma = 13, method = "h")$pval [1] 0.01342352 > senmw(mercury, gamma = 14, method = "t")$pval [1] 0.01990541 > senmw(mercury, gamma = 14, method = "h")$pval [1] 0.03266895 > senmw(mercury, gamma = 16, method = "t")$pval [1] 0.05208875 > senmw(mercury, gamma = 16, method = "h")$pval [1] 0.1187912 ###################################################### # also take a look at point estimates and CI # another way to calibrate > # effect of gamma on point est, two-sided CI > senmwCI(mercury, gamma = 5, method = "h", one.sided = FALSE) $PointEstimate minimum maximum 0.93 4.09 $Confidence.Interval minimum maximum 0.74 4.66 #smallest point est in range still well above 0 for a big gamma > senmwCI(mercury, gamma = 8, method = "h", one.sided = FALSE) $PointEstimate minimum maximum 0.62 5.12 $Confidence.Interval minimum maximum 0.40 5.97 # still outside 0 for really big gamma > senmwCI(mercury, gamma = 10, method = "h", one.sided = FALSE) $PointEstimate minimum maximum 0.47 5.71 $Confidence.Interval minimum maximum 0.24 6.67 > # not even gamma = 10 will move point est or CI to zero; consistent with p-values > # sec 3.1 3,2 of vignette ################### part b > # section 6.3 dispersion weighting # method w weights rows (subclasses) with larger discrepancies > senmw(mercury, method = "w", gamma = 16)$pval [1] 0.009092439 > senmw(mercury, method = "w", gamma = 19)$pval [1] 0.03644606 > senmw(mercury, method = "w", gamma = 17)$pval [1] 0.01540211 > senmw(mercury, method = "w", gamma = 18)$pval [1] 0.02437851 > ## get greater 'power' in the sense of you can tolerate even a bigger gamma before # p-values get too large (here 18). # or for same gamma, resulting p-values may be smaller with weighting than without ### end RQ1