# # weil_eda_571.q # # -------------------------------------------------- data <- matrix(scan("weil.data"), ncol = 3, byrow = T) weil <- list( y = data[, 1], n = data[, 2], x = cbind(1, data[, 3]), z = matrix(1, nrow(data), 1), id = c(1:nrow(data)) ) # ##### define (y0,n0) and (y1,n1) # y0 <- weil\$y[ weil\$x[,2]==0 ] n0 <- weil\$n[ weil\$x[,2]==0 ] y1 <- weil\$y[ weil\$x[,2]==1 ] n1 <- weil\$n[ weil\$x[,2]==1 ] # ooo<- order(y0/n0) out0 <- t( cbind( y0, n0 )[ooo,] ) dimnames( out0 )[[2]]<- paste("R",c(1:16),sep="") out0 ooo<- order(y1/n1) out1 <- t( cbind( y1, n1 )[ooo,] ) dimnames( out1 )[[2]]<- paste("R",c(1:16),sep="") out1 # ### ### means and contrast with w_i = n_i ### # ctl.mean <- sum( y0 )/sum( n0 ) tx.mean <- sum( y1 )/sum( n1 ) # cat("[1] weights = n_i\n") cat(paste("control mean =", round( ctl.mean, 3 ), "\n")) cat(paste(" tx mean =", round( tx.mean, 3 ), "\n")) cat(paste(" contrast =", round( logit(tx.mean) - logit(ctl.mean), 3 ),"\n")) # ### ### means and contrast with w_i = 1 ### # ctl.mean <- mean( y0/n0 ) tx.mean <- mean( y1/n1 ) # cat("\n\n[2] weights = 1\n") cat(paste("control mean =", round( ctl.mean, 3 ), "\n")) cat(paste(" tx mean =", round( tx.mean, 3 ), "\n")) cat(paste(" contrast =", round( logit(tx.mean) - logit(ctl.mean), 3 ),"\n")) # ### ### weighted means as a function of alpha ### # alpha <- c(0:5)/5 m <- length(alpha) p0 <- sum( y0 )/sum( n0 ) p1 <- sum( y1 )/sum( n1 ) # cat("\nWeighted means as a function of alpha\n\n") for( j in 1:m ){ aaa <- alpha[j] w0 <- n0 / ( p0*(1-p0)*( 1 + aaa*(n0-1) ) ) w0 <- w0/sum(w0) mu0 <- sum( w0 * (y0/n0) ) w1 <- n1 / ( p1*(1-p1)*( 1 + aaa*(n1-1) ) ) w1 <- w1/sum(w1) mu1 <- sum( w1 * (y1/n1) ) cat( paste("alpha =", aaa, ": mu0 =", round(mu0,3), ", mu1 =", round(mu1,3),"\n") ) } # # end-of-file...