stat_data_analysis_rcode

Statistical Data Analysis - Fall 2018
Log | Files | Refs | README

R-IC08.R (4781B)

      1 #*******************************************************#
      2 #*******************************************************#
      3 # 		         STAT320 - IC06 Solutions	      	  	    #
      4 #*******************************************************#
      5 #*******************************************************#
      6 #***************************************************#
      7 #           PVal/Power for Two Sample Proportions         #  
      8 #***************************************************#
      9 
     10 #***************************************************#
     11 #***************************************************#
     12 #***************************************************#
     13 #***************************************************#
     14 #***************************************************#
     15 
     16 #***************************************************#
     17 # Means P-value
     18 #***************************************************#
     19 #set.seed(300)
     20 #***************************************************#
     21 
     22 two.prop.permute.pval.f = function(xA,nA,xB,nB,LowTail = TRUE,NSim=100000)
     23 {
     24   xobs = xA/nA - xB/nB
     25   # Reconstruct sample A
     26   xsampA = rep(0,nA)
     27   xsampA[1:xA] = rep(1,xA)
     28   
     29   # Reconstruct sample B
     30   xsampB = rep(0,nB)
     31   xsampB[1:xB] = rep(1,xB)
     32   
     33   # Combine Samples to one super sample
     34   supersample = c(xsampA,xsampB)
     35   
     36   # Create object where simulated DIFFERENCES of averages will be stored
     37   ODiffProps = numeric(NSim)
     38   
     39   # Repeat simulations NSim times
     40   for(i in 1:NSim)
     41   {
     42     
     43     # Shuffle the elements inthe supersample
     44     shuffledsupersample = sample(supersample)
     45     
     46     # Generate sample 1 by resampling super Sample WITHOUT replacement
     47     # in practice we select the first nA elements of the shuffled supersample
     48     NewSampleA = shuffledsupersample[1:nA]
     49     
     50     # Assign remaining values in supersample to Sample B
     51     # that is the elements (nA+1) through (nA+nB)
     52     NewSampleB  = shuffledsupersample[nA+1:nB]
     53     
     54     # Save difference in means simulated samples
     55     ODiffProps[i] = mean(NewSampleA) - mean(NewSampleB )
     56   }
     57   # Compute p-value
     58   if(LowTail==TRUE)
     59   {
     60     # if alternative Ha is mu < mu0
     61     p = sum(ODiffProps<=xobs)/NSim
     62   }else
     63   {
     64     # if alternative Ha is mu > mu0
     65     p = sum(ODiffProps>=xobs)/NSim
     66   }
     67   return(p) 
     68 }
     69 
     70 two.prop.bootstrap.pval.f = function(xA,nA,xB,nB,LowTail = TRUE,NSim=100000)
     71 {
     72   xobs = xA/nA - xB/nB
     73   # Reconstruct sample A
     74   xsampA = rep(0,nA)
     75   xsampA[1:xA] = rep(1,xA)
     76   
     77   # Reconstruct sample B
     78   xsampB = rep(0,nB)
     79   xsampB[1:xB] = rep(1,xB)
     80   
     81   # Combine Samples to one super sample
     82   supersample = c(xsampA,xsampB)
     83   
     84   # Create object where simulated DIFFERENCES of averages will be stored
     85   ODiffProps = numeric(NSim)
     86   
     87   # Repeat simulations NSim times
     88   for(i in 1:NSim)
     89   {
     90     
     91     # Generate sample 1 by resampling super Sample WITH replacement
     92     NewSampleA = sample(supersample,size = nA,replace = TRUE) 
     93     
     94     # Generate sample 2 by resampling super Sample WITH replacement
     95     NewSampleB = sample(supersample,size = nB,replace = TRUE) 
     96     
     97     # Save difference in means simulated samples
     98     ODiffProps[i] = mean(NewSampleA) - mean(NewSampleB )
     99   }
    100   # Compute p-value
    101   if(LowTail==TRUE)
    102   {
    103     # if alternative Ha is mu < mu0
    104     p = sum(ODiffProps<=xobs)/NSim
    105   }else
    106   {
    107     # if alternative Ha is mu > mu0
    108     p = sum(ODiffProps>=xobs)/NSim
    109   }
    110   return(p) 
    111 }
    112 
    113 
    114 #***************************************************#
    115 # Means CI
    116 #***************************************************#
    117 
    118 two.prop.CI.f = function(xA,nA,xB,nB,alpha=.05,NSim=100000)
    119 {
    120   # Reconstruct sample A
    121   xsampA = rep(0,nA)
    122   xsampA[1:xA] = rep(1,xA)
    123   
    124   # Reconstruct sample B
    125   xsampB = rep(0,nB)
    126   xsampB[1:xB] = rep(1,xB)
    127   
    128   # Create object where simulated DIFFERENCES of averages will be stored
    129   ODiffProps = numeric(NSim)
    130   
    131   # Repeat simulations NSim times
    132   for(i in 1:NSim)
    133   {
    134     # Generate sample 1 by resampling original Sample A WITH replacement
    135     NewSampleA = sample(xsampA,size = nA,replace = TRUE) 
    136     
    137     # Generate sample 2 by resampling original Sample B WITH replacement
    138     NewSampleB = sample(xsampB,size = nB,replace = TRUE) 
    139     
    140     # Save difference in means simulated samples
    141     ODiffProps[i] = mean(NewSampleA) - mean(NewSampleB)
    142   }
    143   
    144   # Compute CI
    145   ci =quantile(ODiffProps,probs = c(alpha/2,1-alpha/2))
    146    return(ci) 
    147 }
    148 
    149 #***************************************************#
    150 # Questions
    151 #***************************************************#
    152 
    153 xA = 5
    154 nA = 119
    155 
    156 xB = 19
    157 nB = 115
    158 
    159 LowTail = TRUE
    160 NSim = 100000
    161 
    162 alpha=.3
    163 
    164 two.prop.bootstrap.pval.f(xA, nA, xB, nB, LowTail, NSim)
    165 two.prop.permute.pval.f(xA, nA, xB, nB, LowTail, NSim)
    166 two.prop.CI.f(xA, nA, xB, nB, alpha, NSim)