# com9b of February 2016, # exercise with a 2 x 2 table, data from Vermont # 44 women and 101 men have been asked about whether # they are positive to same-sex marriages or not, # for the US state Vermont. For the women, 35 said "yes", 9 "no". # For the men, 60 said "yes", 41 "no". # think in terms of X = (men,woman), Y = (yes,no) # yes no ## 60 41 men p11 p12 ## 35 9 women p21 p22 # here I compare pmen and pwomen, Pr(yes) # rho = pwomen/pmen NN <- rbind(c(60,41),c(35,9)) nn <- sum(NN) phat <- NN/nn nmen <- sum(NN[1, ]) # 101 nwomen <- sum(NN[2, ]) # 44 ymen <- NN[1,1] # 60 ywomen <- NN[2,1] # 35 pval <- seq(0,1,by=0.001) Cmenval <- 1 - pbinom(ymen,nmen,pval) + 0.5*dbinom(ymen,nmen,pval) Cwomenval <- 1 - pbinom(ywomen,nwomen,pval) + 0.5*dbinom(ywomen,nwomen,pval) ccmenval <- abs(1-2*Cmenval) ccwomenval <- abs(1-2*Cwomenval) matplot(pval,cbind(Cmenval,Cwomenval),type="l") matplot(pval,cbind(ccmenval,ccwomenval),type="l") logL <- function(pp) { pm <- pp[1] pw <- pp[2] # log(dbinom(ymen,nmen,pm)) + log(dbinom(ywomen,nwomen,pw)) } minuslogL <- function(pp) {-logL(pp)} nlm(minuslogL,c(0.5,0.5),hessian=T) # 0.594059 0.795454 for men, women ## now rho = pwomen / pmen rhoval <- seq(0.75,1.75,by=0.01) logLprofval <- 0*rhoval for (j in 1:length(rhoval)) { rhonow <- rhoval[j] logLnow <- function(pm) { log(dbinom(ymen,nmen,pm)) + log(dbinom(ywomen,nwomen,pm*rhonow)) } # minuslogLnow <- function(pm) {-logLnow(pm)} # nilsnow <- nlm(minuslogLnow,0.5,hessian=T) # logLprofval[j] <- logLnow(nilsnow$estimate) print(c(j,rhonow,logLprofval[j])) } devval <- 2*(max(logLprofval)-logLprofval) ccval <- pchisq(devval,1) matplot(rhoval,ccval,type="l") matlines(rhoval,0.95+0*rhoval,type="l",lty=2)