#### Week 6, Sensitivity calculations, DOS Ch 3, Rosenbaum vignette


> install.packages("sensitivitymw")
> data(erpcp) # Welders DNA damage
> dim(erpcp) # # data are outcomes in wide form; each row is a subclass
[1] 39  2
> head(erpcp)
  welder control
1  0.899   0.751
2  1.233   0.875
3  1.733   0.161
4  3.156   0.630
5  1.749   1.462
6  0.431   0.702
> attach(erpcp)
> boxplot(welder,control)# matches DOS Fig 3.1

> t.test(erpcp$welder, erpcp$control)
        Welch Two Sample t-test
data:  erpcp$welder and erpcp$control
t = 5.1442, df = 54.368, p-value = 3.785e-06
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 0.3502495 0.7974940
sample estimates:
mean of x mean of y 
1.3957436 0.8218718 

> wilcox.test(erpcp$welder, erpcp$control)
        Wilcoxon rank sum test with continuity correction
data:  erpcp$welder and erpcp$control
W = 1251.5, p-value = 9.497e-07
alternative hypothesis: true location shift is not equal to 0
> wilcox.test(erpcp$welder - erpcp$control)
        Wilcoxon signed rank test
data:  erpcp$welder - erpcp$control
V = 715, p-value = 6.247e-07
alternative hypothesis: true location is not equal to 0

#### p.4 Gamma and p-values
>  senmw(erpcp, gamma = 1, method = "t")$pval
[1] 2.048115e-05
>  senmw(erpcp, gamma = 2, method = "t")$pval
[1] 0.003737467
>  senmw(erpcp, gamma = 3, method = "t")$pval
[1] 0.02275942
>  senmw(erpcp, gamma = 4, method = "t")$pval
[1] 0.0579339
> # I think doubling 'p' is right



> #### now to senmwCI page 5
>  senmwCI(erpcp, gamma = 1, method = "t", one.sided = TRUE)
$PointEstimate
minimum maximum 
 0.5739  0.5739 
$Confidence.Interval
minimum maximum 
  0.394     Inf 

>  senmwCI(erpcp, gamma = 1, method = "t", one.sided = FALSE) # get two-sided CI
$PointEstimate
minimum maximum 
 0.5739  0.5739 
$Confidence.Interval
minimum maximum 
 0.3561  0.7916 

>  senmwCI(erpcp, gamma = 2, method = "t", one.sided = TRUE)
$PointEstimate
minimum maximum 
 0.4167  0.7487 
$Confidence.Interval
minimum maximum 
 0.2081     Inf 

>  senmwCI(erpcp, gamma = 2, method = "t", one.sided = FALSE)
$PointEstimate
minimum maximum 
 0.4167  0.7487 
$Confidence.Interval
minimum maximum 
 0.1552  1.0414 

### now try the "amplify stuff" (2009 JASA paper)
> install.packages("sensitivitymv")
##### page 6  sec 2.4
### Gamma = f(Lamda [bias in assignment], Delta [bias in outcome])
### amplify takes arguments Gamma, Lamda(values >Gamma)  produces Delta
> amplify(3,c(4:7))
        4         5         6         7 
11.000000  7.000000  5.666667  5.000000 
## discussion p.7




> ## mercury in RQ