Raising Statistical Standards Effect on Sample Size

The failure of mainstream research to consistently reproduce results have led many to look for the faults in current methodologies.

One of these potential faults identified is that the significance levels of current standards is too high.  A standard rejection rate of either .05 or .01 is not high enough.

Statistician Valen Johnson recently released an article in The National Proceedings Academy of Sciences which reccomends more appropriate standard of rejection being .005 or .001.

In this post I will attempt to examine how such a change could affect the required sample size for studies.

Before initiating an experimental study often a power analysis is done.  When possible it is using a relatively simple closed form numerical statistic. This is the case when relatively simple methods are intended to be used.  More complex methods often require the use of simulations to do a power analysis.

Usually the logic of a power analysis (as far as I know) goes something like.  Let's say the possible effect size is tau and we know the conditional distribution of outcomes (from previous work) has a standard deviation of SD. How many people would be need to reject the null at our intended level.

Referencing wikipedia:

With Pi being power and tau being effect size and alpha being rejection power.  Assuming lots of normality:

Pi(tau) = 1 - PDFNORMAL(Z(alpha)-tau*N^.5/SD)

We want to solve for N:

# PDFNORMAL^(-1)(1 - Pi) =  Z(alpha)-tau*N^.5/SD

# tau*N^.5/SD = PDFNORMAL^(-1)(1 - Pi) - Z(alpha)

# N = ((SD/tau)*(PDFNORMAL^(-1)(1 - Pi) - Z(alpha)))^2

samp.power <- alpha="" function="" p="" pi="" tau="">

  ceiling(((SD/tau)*(pnorm(1 - pi) - qnorm(alpha)))^2)

# added the ceiling function to round up.

Let's see it in action!

Let's say we have an outcome which we know has a SD=2 and we hope our effect will have at least a size of tau=1. Following standard practices we require a detection rate of 80% for our power analysis.  Let's see what happens when we vary our alpha rate!

samp.power(SD=2, tau=1, pi=.8, alpha=.05)

# 20

samp.power(SD=2, tau=1, pi=.8, alpha=.01)

# 34

samp.power(SD=2, tau=1, pi=.8, alpha=.005)

# 40

samp.power(SD=2, tau=1, pi=.8, alpha=.001)

# 54

We can see that increasing our rejection standards from .05 to .001 we are basically increasing our required sample pool by 2.7. Which is really not that bad.  Looking at a smaller effect size does not change things except by a squared factor s. tau'=tau/s

# Looking at our equation for N

# N = ((SD/tau)*(pnorm(1 - pi) - qnorm(alpha)))^2

# N = (1/tau * H)^2

# where H = (SD)*(pnorm(1 - pi) - qnorm(alpha))

# Thus substituting in tau'

# N = (1/tau/s * H)^2 = s^2 * (H/tau)^2

# So let's say s=10 then N(tau')=N(tau)*100

samp.power(SD=2, tau=1/10, pi=.8, alpha=.001)

# 5387

Not to make light of this new proposed rejection rate and its larger sample size.  By increasing the sample size by roughly a factor of 2.7 the cost of a study might easily double.

However, a researcher saying that the chance of an outcome occurring randomly going from 1 out of 20 to 1 out of 1000 might easily be worth the additional cost.

The nice thing about this new approach would be that it would still allow for less strong rejections even when the effect size is smaller than expected or when there is more noise in the sample than expected.

Well, that is what I have to say in support of the idea.  I also have some reservations.  If the cost of the larger rejection rates is really doubling the cost of the study then why not do two studies?  Assuming the outcomes of each study are random and iid the likelihood of rejecting the null 1 out of 20 (5%) of the time given two studies is 1 out of 20^2 or 1 out of 400.

.25% (1/400) is not as good as as .1% but it is still pretty strong.

It makes sense in a world where we do not know what really has an effect to have more smaller studies which we follow up with larger studies when we do find an effect.  In addition, there might be factors unique to individual studies which are for whatever reason unobservable and nonreproducable, driving the results.

Say, the researchers introduced bias without intending to.  Scaling up the project might have no effect on removing or controlling that bias.  However, having two different studies run by different research groups is less likely to reintroduce the same bias.

Overall, I think it may be useful to introduce higher standards into social science research especially in non-experimental data in which numerous potential researchers are looking at the data with different hypothesizes.  It is improbable that if there is enough researchers looking at the data from enough angles that there will not be at least a few that reject at a 5% level.  Imagine that you have 20 research teams each picking 5 different angles that is 1000 different draws. 

Assuming they are independent the likelihood of rejecting at a 5% level would lead researchers to falsely reject 50 null hypothesizes on average.  That is a lot of false rejections.  Choosing a level of .1% however would lead only one research team to reject one null falsely on average.  This is a pretty appealing change for a conservative statistician.   However, once again there will be many times in which we fail to reject the null when our hypothesizes are in fact true but our data or effect set are insufficiently large.  Which would we rather we have?

Simulation of Blackjack: the odds are not with you

I often want to simulate outcomes varying across a set of 

parameters. In order to accomplish this in an efficient 

manner I have coded up a little function that takes parameter

vectors and produces results. First I will show how to set it

up with some dummy examples and next I will show how it can 

be used to select the optimal blackjack strategy.

 
SimpleSim <- function(..., fun, pairwise=F) {
  # SimpleSim allows for the calling of a function varying 

  # multiple parameters entered as vectors. In pairwise form 

  # it acts much like apply. In non-paiwise form it makes a 

  # combination of each possible parameter mix

  # in a manner identical to block of nested loops.
 
  returner <- NULL
  L        <- list(...)
  # Construct a vector that holds the lengths of each object
  vlength <- unlist(lapply(L, length)) 
  npar    <- length(vlength)
  CL      <- lapply(L, "[", 1) # Current list is equal to the first element
 
 # Pairwise looping
 if (pairwise) {
  # If pairwise is selected than all elements greater than 1 must be equal.
  # Checks if all of the elements of a vector are equal
  if (!(function(x) all(x[1]==x))(vlength[vlength>1])) {
   print(unlist(lapply(L, length)))
   stop("Pairwise: all input vectors must be of equal length", call. =F)
  }
  for (i in 1:max(vlength)) { # Loop through calling the function
   CL[vlength>1]  <- lapply(L, "[", i)[vlength>1] # Current list
   returner <- rbind(returner,c(do.call(fun, CL),pars="", CL))
  }
 } # End Pairwise
 
 # Non-pairwise looping
 if (!pairwise) {
  ncomb <- prod(vlength) # Calculate the number of combinations
  print(paste(ncomb, "combinations to loop through"))
  comb <- matrix(NA, nrow=prod(vlength), ncol=npar+1)
  comb[,1] <- 1:prod(vlength) # Create an index value
  comb <- as.data.frame(comb) # Converto to data.frame
  colnames(comb) <- c("ID", names(CL))
  for (i in (npar:1)) { # Construct a matrix of parameter combinations
   comb[,i+1] <- L[[i]] # Replace one column with values
   comb<-comb[order(comb[,(i+1)]),] # Reorder rows
  }
  comb<-comb[order(comb[,1]),]
  for (i in 1:ncomb) {
   for (ii in 1:npar) CL[ii] <- comb[i,ii+1]
   returner <- rbind(returner,c(do.call(fun, CL),pars="", CL))
  }
 } # End Non-Pairwise
 
 return(returner)
 
} # END FUNCTION DEFINITION
 
# Let's first define a simple function for demonstration
minmax <- function(...) c(min=min(...),max=max(...))
 
# Pairwise acts similar to that of a multidimensional apply across columns 
SimpleSim(a=1:20,b=-1:-20,c=21:40, pairwise=T, fun="minmax")
# The first set of columns are those of returns from the function "fun" called.
# The second set divided by "par" are the parameters fed into the function.
SimpleSim(a=1:20,b=-1:-20,c=10, pairwise=T, fun="minmax")
 
# Non-pairwise creates combinations of parameter sets.
# This form is much more resource demanding.
SimpleSim(a=1:5,b=-1:-5,c=1:2, pairwise=F, fun="minmax")
 
# Let's try something a little more interesting.
 
# Let's simulate a game of black jack strategies assuming no card counting is possible.
blackjack <- function(points=18, points.h=NULL, points.ace=NULL, 
                      cards=10, cards.h=NULL, cards.ace=NULL,
                      sims=100, cutoff=10) {
  # This function simulates a blackjack table in which the player
  # has a strategy of standing (not asking for any more cards)
  # once he has either recieved a specific number of points or 
  # a specific number of cards.  This function repeates itself sims # of times.
  # This function allows for up to three different strategies to be played.
  # 1. If the dealer's hole card is less than the cuttoff
  # 2. If the dealer's hole card is greater than or equal to the cuttoff
  # 3. If the dealer's hole card is an ace
  # In order to use 3 level strategies input parameters as .h and .ace
 
  # It returns # of wins, # of losses, # of pushes (both player and dealer gets 21)
  # and the number of blackjacks.
 
  # This simulation assumes the number of decks used is large thus
  # the game is like drawing with replacement.
 
  if (is.null(points.h))   points.h   <- points
  if (is.null(points.ace)) points.ace <- points.h
  if (is.null(cards.h))    cards.h    <- cards
  if (is.null(cards.ace))  cards.ace  <- cards.h
 
  bdeck <- c(11,2:9,10,10,10,10) # 11 is the ace
 
  bdresult <- c(ppoints=NULL, pcards=NULL, dpoints=NULL, dcards=NULL)
 
  for (s in 1:sims) {
   dhand <- sample(bdeck,1) # First draw the deal's revealed card
   phand <- sample(bdeck,2, replace=T)
 
   # Specify target's based on dealer's card
   if (dhand<cutoff) {
     pcuttoff <- points
     ccuttoff <- cards
   }
   if (dhand>=cutoff) {
     pcuttoff <- points.h
     ccuttoff <- cards.h
   }
   if (dhand==11) {
     pcuttoff <- points.ace
     ccuttoff <- cards.ace
   }
 
   # player draws until getting above points or card count
   while ((sum(phand)<pcuttoff)&(length(phand)<ccuttoff)){
     phand <- c(phand, sample(bdeck,1))
       # If player goes over then player may change aces to 1s 
       if (sum(phand)>21) phand[phand==11] <- 1
   }
 
   # Dealer must always hit 17 so hand is predetermined
   while (sum(dhand)<17) {
     dhand <- c(dhand, sample(bdeck,1))
     # If dealer goes over then dearler may change aces to 1s
     if (sum(dhand)>21) dhand[dhand==11] <- 1
   }
   bdresult <- rbind(bdresult, 
        c(ppoints=sum(phand), pcards=length(phand), 
          dpoints=sum(dhand), dcards=length(dhand)))
  }
 
  # Calculate the times that the player wins, pushes (ties), and loses
  pbj <- (bdresult[,1]==21) & (bdresult[,2]==2)
  dbj <- (bdresult[,3]==21) & (bdresult[,4]==2)
  pwins <- ((bdresult[,1] >  bdresult[,3]) & (bdresult[,1] <  22)) | (pbj & !dbj)
  push  <- (bdresult[,1] == bdresult[,3]) | (pbj & dbj)
  dwins <- !(pwins | push)
 
  # Specify the return.
  c(odds=sum(pwins)/sum(dwins), 
    pwins=sum(pwins), 
    dwins=sum(dwins), 
    push=sum(push), 
    pcards=mean(bdresult[,2]), 
    dcards=mean(bdresult[,4]),
    pblackjack=sum(pbj),
    dblackjack=sum(dbj))
}
 
blackjack(points=18, sims=4000)
# We can see unsurprisingly, that the player is not doing well.
 
blackjack(points=18, points.h=19, sims=4000)
# We can see that by adopting a more aggressive strategy for when
# the dealer has a 10 point card or higher, we can do slightly better.
# But overall, the dealer is still winning about 3x more than us.
 
# We could search through different parameter combinations manually to
# find the best option.  Or we could use our new command SimpleSim!
 
MCresults <- SimpleSim(fun=blackjack, points=15:21, points.h=18:21, 
                       points.ace=18:21, cutoff=9:10, cards=10, sims=100)
 
# Let's now order our results from the most promising.
MCresults[order(-unlist(MCresults[,1])),]
 
# By the simulation it looks like we have as high as a 50% ratio of loses to wins.
# Which means for every win there are 2 loses.
# However, I don't trust it since we only drew 100 simulations.
# In addition, this is the best random draw from all 224 combinations which each
# have different probabilities.
 
# Let's do the same simulation but with 2000 draws per.
# This might take a little while.
MCresults <- SimpleSim(fun=blackjack, points=15:21, points.h=18:21, 
                       points.ace=18:21, cutoff=9:10, cards=10, sims=5000)
 
# Let's now order our results from the most promising.
MCresults[order(-unlist(MCresults[,1])),]

hist(unlist(MCresults[,1]), main="Across all combinations\nN(Win)/N(Loss)", 

     xlab = "Ratio", ylab = "Frequency")

# The best case scenario 38% win to loss ratio appears around were we started, 
# playing to hit 18 always and doing almost as well when the dealer is high
# (having a 10 or ace) then playing for 19.
 
# Overall, the odds are not in our favor.  For every win we expect 1/.38 (2.63) loses.

Highlighted by Pretty R at inside-R.org

MSimulate Command

CLT is a powerful property

* Stata has a wonderfully effective simulate function that allows users to easily simulate data and analysis in a very rapid fashion.

* The only drawback is that when you run it, it will replace the data in memory with the simulated data automatically.

* Which is not a big problem if you stick a preserve in front of your simulate command.

* However, you may want to run sequential simulates and keep the data form all of the simulations together rather than only temporarily accessed.

* Fortunately we can accomplished this task by writing a small program.

cap program drop msim
program define msim

  * Gettoken will split the arguments fed into msim into those before colon and those after.
  gettoken before after : 0 , parse(":")

  * I really like this feature of Stata!

  * First let's strip the colon.  The 1 is important since we want to make sure to only remove the first colon.
  local simulation = subinstr("`after'", ":", "", 1)

  * Now what I propose is that the argument in `before' is used as an extension for names of variables created by the simulate command.

  * First let's save the current data set.


  * Generate an id that we will later use to merge in more data
  cap gen id = _n

  * Save the current data to a temporary location
  tempfile tempsave
  save `tempsave'

  * Now we run the simulation which wipes out the current data.
  `simulation'

  * First we will rename all of the variables to have an extension equal to the first argument
  foreach v of varlist * {
    cap rename `v' `v'`before'
  }

  * Now we need to generate the ID to merge into

  cap gen id = _n

  merge 1:1 id using `tempsave'
  * Get rid of the _merge variable generated from the above command.
  drop _merge
end


* Let's write a nice little program that we would like to simulate.
cap program drop simCLT
program define simCLT

  clear
  set obs `1'
  * 1 is defined as the first argument of the program sim

  * Let's say we would like to see how many observations we need for the central limit theorem (CLT) to make the means of a bernoulli distribution look normal.  Remember, so long as the mean and variance is defined the generally central limit theorem will eventually force any random distribution of means to approximate a normal distribution as the number of observations gets large.
  gen x = rbinomial(1,.25)

  sum x
end

* So let's see first how the simulate command works initially

simulate, rep(200): simCLT 100

* The simulate command will automatically save the returns from the sum command as variables (at least in version 12)

hist mean, kden
* The mean is looking good but not normal

* Now normally what we need to do would be to run simulate again with a different argument.

* But instead let's try our new command with 200!


* But instead let's try our new command!
clear
* Clear out the old results
msim 100: simulate, rep(200): simCLT 100

msim 200: simulate, rep(200): simCLT 200

* Looks good!

msim 400: simulate, rep(200): simCLT 400

msim 1000: simulate, rep(200): simCLT 1000

msim 10000: simulate, rep(200): simCLT 10000

msim 100000: simulate, rep(200): simCLT 100000

msim 1000000: simulate, rep(200): simCLT 100000

* The next two commands can take a little while.
msim 10000000: simulate, rep(200): simCLT 1000000

msim 100000000: simulate, rep(200): simCLT 10000000

sum mean*
* We can see that the standard deviations are getting smaller with a larger sample size.

* How is the histograms looking?
foreach v in 100 200 400 100 1000 10000 100000 1000000 10000000 {
  hist mean`v', name(s`v', replace)  nodraw  title(`v') kden

}

graph combine s100 s200 s400 s1000 s1000 s10000 s100000 s1000000 s10000000 ///
  , title("CLT between 100 and 10,000,000 observations")

* We can see that the distribution of means approximates the normal distribution as the number of draws in each sample gets large.

* This is one of the fundamental findings of statistics and pretty awesome if you think about it.

Checking for differences in estimates by group

Stata do file

* Imagine that you have two groups in your sample and you are wondering if both groups respond the explanatory variables in the same manner.

clear
set obs 1000

gen group=rbinomial(1,.5)

gen x = rnormal()

gen y = 2 + 3*x + rnormal()*2 if group==0
replace y = 1 + 0*x + rnormal()*2 if group==1

* There is several ways to test if the two sets of estimators are the same.

* One method is using interaction dummies.

gen x_group = x*group

* The technique is to include an interaction term in the model for each coefficient that you would like to check if it is equal.

* Group counts as an interaction term between the constant and the group indicator.
* x_group is the interaction between the x variable and the group variable.

reg y group x_group x

* In order to interpret the coefficients in the above regression.

* First set group equal to zero.

* In that case all that is in the regression is: y = 1*b0 + x*b1 which is the estimator for the group 0.

* We can see the constant estimate is close to 2 and the coefficient on x is close to 3 which are our values.

* Now the tricky thing is to interpret the remaining coefficients.

* Because the group coefficient is effectively interacted with the constant then in order to estimate the constant in the group=1 set just add _b[group- to _b[cons].

di "Group 2 constant estimate is " _b[group]+ _b[_cons] " which we can see is close to 1"

* We do a similar procedure with the slope.

di "Group 2 slope estimate is " _b[x]+ _b[x_group] " which we can see is close to 0"

* The nice thing about this formulation is that the coefficients are already included in the estimation.

* Thus we can easily test the hypothesis that both are the same manner.

test x_group=group=0

* If both groups were the same then they would have the same mean (same constants) and respond the same the the treatment x.

* An alternative way of testing each coefficient individually is:

reg y x if group == 0

local b_cong0 = _b[_cons]
local b_xg0 = _b[x]
local se_cong0 = _se[_cons]
local se_xg0 = _se[x]

reg y x if group == 1

* The t-test that they are the same is b[group=0]-b[group=1]=0
* var(b[group=0]-b[group=1])= var(b[group=0]) + var(b[group=1])

local joint_var_cons = `se_cong0'^2 + _se[_cons]^2

* And the t test is (b[group=0]-b[group=1])/(var(b[group=0]) + var(b[group=1]))^.5

di "t_constants = " (`b_cong0'- _b[_cons] )/(`joint_var_cons')^.5

local joint_var_x= `se_xg0'^2 + _se[x]^2

di "t_constants = " (`b_xg0'-_b[x])/`joint_var_x'^.5

* This formation assumes that the estimates vary independently.

* If done properly the t-stats are almost indentical between formulations.

* By the way, please fill out my survey posted on the last post.

Path Analysis

Stata do file

* Path analysis is an interesting statistical method that can be used to indentify complex relationships beween variables and an outcome variable.

* As with all statistical methods the modelling framework is essential to derive reasonable results.

* Conviently, I am only interested in simulating data so as usual my data will perfectly conform to the model's specifications.

* Imagine the following model.

* All of the boxes are observable variables.  The arrows indicate the causal direction of the effects.

* There are two exogenous variables: A and D.   These variables are not influenced by any other variables in the model.

* All other variables are endogenous.

* Each of the variables represents a direct effect of a one unit change in one variable on that of the other variable.

* This framework is convient because it allows us to indentify a "total effect" which is a combined result of both the direct and indirect effects of variables on the outcome variable.

* The variable of primary interest in explaining is H.

* The variable G has only a direct effect on H (pHG).

* While the variable C only has an indirect effect on H (pFC*pHF).

* The reason the indirect effect is a product is because C has a pFC effect on F, and F has a pHF effect on H, thus a change in H as a result of a change in C is how much F changes as a result of C and how much that change effects H.

* Variables can have both and indirect and direct effect.

* B for instance has the direct effect: pHB
* Indirect effects: pCB*pFC*pHF + pEB*pHE
* Total effect: pHB + pCB*pFC*pHF + pEB*pHE

* The key feature about this particular example is that all of the arrows are one directional.

* Making a great deal of inference possible that otherwise would not be possible.

* Usually we cannot say that when trying to explain H with explanatory variables A through G that A causes B and B causes H.

* However, if we do the work to indentify reasonable pathways then this type of analysis could be quite interesting.

* Let's generate out data.

clear

* Let's imagine 6000 youth in our sample.

set obs 6000

* Let's first specify our effects

* pEA = .3
* pEB = .13

* pHA = .2
* pHB = .2
* pHE = .3
* pHG = 1.1
* pHF = .2

* pBA = .5

* pCB = .2
* pCD = .1

* pGD = .2

* pFC = .76
* pFB = .4

* For B we can calculate our true effects:

* B Direct: pHB = .2

* Indirect effects: pCB*pFC*pHF + pEB*pHE
* Indirect effects:.2*.76*.2 + .13*.3 = .0694

* Total effect: .0694+.2 = .2694

gen A =                rnormal()
gen B = A*.5 +         rnormal()
gen D =                rnormal()
gen C = B*.2 + D*.1 +  rnormal()
gen E = A*.3 + B*.13 + rnormal()
gen F = B*.4 + C*.76 + rnormal()
gen G = D*.2 +         rnormal()
gen H = E*.3 + A*.2 + B*.2 + F*.2 + G*1.1 + rnormal()

* Simualtion Done

* In order to generate our different effects we simply run OLS for each endogenous variable.

reg A B
  local pBA = _b[B]

reg C B D
  local pCB = _b[B]
  local pCD = _b[D]

reg C B D
  local pCB = _b[B]
  local pCD = _b[D]

reg G D
  local pGD = _b[D]

reg F C B
  local pFB = _b[B]
  local pFC = _b[C]

reg E A B
  local pEA = _b[A]
  local pEB = _b[B]

reg H A B E F G
  local pHA = _b[A]
  local pHB = _b[B]
  local pHE = _b[E]
  local pHF = _b[F]
  local pHG = _b[G]

* In order to estimate the indirect effect say of B on H.
* We just plug our estimates into the equation.

* B direct effect: pHB
* Indirect effects: pCB*pFC*pHF + pEB*pHE
* Total effect: pHB + pCB*pFC*pHF + pEB*pHE

di "B's estimated indirect effect = `pCB'*`pFC'*`pHF' + `pEB'*`pHE'"
di "B's estimated indirect effect = " `pCB'*`pFC'*`pHF' + `pEB'*`pHE'

* Which turns out to be close to our true value.

di "B's total estimated effect on H is " `pHB' + `pCB'*`pFC'*`pHF' + `pEB'*`pHE'

* It is possible to use the user written command pathreg to make things easier.

* Install it by typing the following command. findit pathreg
pathreg (H E B F G) (G D) (C B D) (B A) (E A B) (F B C)

* This command does not currently calculate out all of the indirect and direct effects.

* I am not sure the best way to calculate the standard errors of the different effect estimates.

* My guess is that since this is just a series of fast OLS regressions the easiest thing to do would be to boostrap the entire process.

* This would require slightly more code but definitely easy to do from this point.

R-squared

Original Code

* R-Squared

* R-squared and pseudo r-squared is a useful statistics produced by most regression type estimation routines.

* R-squared (R2) is a measure of how much of the variance in y is explained by the model.

* Thus a model with only an intercept has an R2 of 0.

set seed 101
clear
set obs 10000

gen y1=rnormal()

reg y1

* While in the opposite extreme a model which does not have any unexplained variance has an r2 of 1.

reg y1 y1, noconstant

* Technically this regression should not work but Stata does the math and produces the results.

* Let's see how well R2 approximates explainable variation.

gen x=rnormal()
gen u=rnormal()

gen y2 = (1)*x + u
* The variance of the model is equal to 1 (from 1^2 * var(x))
* The variance of the unexplained error is equal to 1 (var(u))
* Thus our true explained variance should be equal to var(x)/(1^2 * var(x)+var(u)) = 1/2

reg y2 x
* Thus we can see our R2 estimate of the explained variance is very close to the true which is .5

* I made enphasis on noting the coefficient on the x.

* That coefficient significantly scales explainable variation.

* Thus:
gen y3 = 2*x + u

* Should have a much larger R2 because model variance = 2^2*varx = 4
* Var(u) = 1
* R2 = 4/(4+1)=.8

reg y3 x

* If we were to add multiple xs the calculation is similar though if there is correlation between the xs then that will factor into the model.

gen x1 = rnormal()
gen x2 = rnormal()

gen y4 = x1 + x2 + u

* R2 = var(x1) + var(x2) / (var(x1) + var(x2) + var(u) = 2/3
reg y4 x1 x2

* If there is correlation between the xs then that can substantially throw off the calculations.

* In the extreme cases corr(x1, x2)=1 then we are back to the same scenario as y3

* y = x1 + x2 + u (if x1~N(0,1) and x2~N(0,1)) then x1=x2

* y = 2*x1 + u

* Thus R2 = .8

* In the other extreme corr(x1,x2)=-1

* Then, given that they are both N~(0,1), x2=-x1

* y = x1 + x2 + u = x1 - x1 + u = u

* Which is the same as y1

* R2 = 0

* R-squared can also be thought of as the square of the correlation between the predicted values and the observed.

reg y4 x1 x2
predict y4hat

corr y4hat y4
* Thus we can see that there is an 81% correlation between yhat and y observed.

* A high correlation would indicate that our model have done well at predicting observable characteristics.

di r(rho)^2

* A brief note on adjusted R2.

* R2 is known to always be larger the more variables are in your model.


gen z1 = rnormal()
gen z2 = rnormal()
gen z3 = rnormal()

reg y4 x? z?

* Thus: the R2 moved from   R-squared     =  0.6610 to
*                           R-squared     =  0.6611

* This factor being known researchers have developed the Adj-R2 which slightly penalizes the R2 for including more variables.

* Thus Adj R-squared =  0.6609

* This might be appropriate given known facts, however it is trivial and almost always worth ignoring.

* I generally don't pay attention to the AR2 and I don't know anybody else who does either.

* A .0001 difference in R2 is so unimportant as to be completely ignorable without significant loss of content.

A better way of updating election outcomes

Written 12:30 Tuesday evening.

I don’t understand why pollsters never do this calculation.  At this moment Nevada is officially undecided.  However 80% of votes are cast with Obama leading 52 to Romney 46.  From these numbers we can easily calculate what percentage of voters Romney needs to take Nevada.  Let’s first calculate the number he needs to match Obama.

We can thing of this as the weighted average between the votes already counted and those yet to be counted.  This can be written out as the cast votes .8 and the uncast votes .2.  a is the proportion who vote Romney, (1-a) is the proportion that vote Obama.  Everybody votes in this approximation.

.8 * .46 + .2 * a = .8 *.52 + .2*(1-a)

Now we need only solve for a.  

P*R + (1-P)*a = P*O + (1-P)*(1-a)

P*R - P*O = (1-P)*(1-a) - (1-P)*a = (1-P)*(1-2a)

(P*R - P*O)/(1-P) = 1-2a

a = ½  -  (R - O)/2   * P/(1-P)

R is Romney voters, O Obama and P is the proportional that has already voted.

Using the previous example:

a = .5 - (.46-.52) * .8/.2                                                  

The above formula can be adjusted by making the totals into proportions.  I get Romney would need to win 74 percent of the remaining Nevada votes to win.  Probably unlikely.

A general formula can be defined as:

P*R + (1-P)*a = P*O + (1-P)*(1-a)

Now let’s look at Florida.

Wednesday morning. This state is very close, just using a "1 percent" margin is not sufficiently accurate for me.

PTotal = 4129360+4083321

R = 4083321 / PTotal

O = 4129360/ PTotal

Now we can easily calculate the P value from this.  .  Also, we have 94% of districts reporting. Thus :

P=.97

So Romney only has a P value of .496175.  Finally let’s plug in the values above. 

a = .5  -  (R - O)/2   * P/(1-P)

a = 1/2  - (R-O)/2 * P/(1-P)

Thus we get that Romney would need win 59% of the remaining votes to take Florida.  This as opposed to Nevada is possible though at this point unlikely.

Power Analysis

# Power analysis is a frequently used tool to calculate the minimum sample size neccessary in order to show that a hypothesized result exists given some specified rate of type I and type II error.  A typical level of error for type I is alpha=5%.  Type II error for the purpose of power analysis is often specified at a beta=20% level.

# Let's first look at how selecting a sample corresponds to rejection rates and type I errors. We can do this based on the normality assumption of the underlying distribution of errors.  If we assume that we know that the underlying errors are normally distributed then we can say with exact precision how frequently a difference in means will reject the null for a particular sample size.

# A critical feature of power analysis is the specification of the expected effect size.  The larger the specified effect size the smaller the statistical sample we need in order to reject the null.  Symmetrically, the larger the effect size you assume the more likely you will fail to reject the null in practice as your analysis fails to have sufficient power.  Thus, choosing expected effect size is a non-trivial task.  Sometimes, discipline specific requirements will inform your decision as to the a-priori effect size.  For instance in educational interventions, a new technology that is expensive may not be considered economically feasible unless it increases student performance by .1 on overall GPA.
d = .1

# It is also important to either have information about the population's sample variance or to make some assumption about the populations variance.

# I assume that the sample's standard deviation of GPAs is 1.
s = 1

# I would like to first look at the number of students that may be part of the study.  In this case anywhere between 2 and 500.
n = 2:500

# The standard errror is:
se = s/(n^.5)

# The resulting t score is:
t = d/se

# Probability of rejecting the null is drawn from the Student t-distribution because we are assuming a two sided rejection rate.
p <- pt(t,n)

plot(n, p, ylim=c(.5,1), ylab="Probability", col="blue",
          xlab="Sample Size", type="l", yaxs = "i", xaxs = "i" ,
          lwd=2, main="Probability of Rejecting the Null")
# I choose .975 because 5% significance level on a two sided test leaves only 2.5% in each tail.
lines(y=c(0.975,0.975), x=c(2,500), col="red", lwd=2)

# We can see that for an effect size that is only 10% of the standard deviation of the error we would need nearly 400 people before we can be confident of rejecting the null 95% of the time.  I believe there is a symetric argument for type I errors as well.  I will attempt to address this more later.

# At this point one may wonder if we really need to hassle with the t-statistic.  For the most part you would be right.  The difference between the t and the z is extremely small.

z = d/se
pz <- pnorm(z)
plot(n, p, ylim=c(.5,1), ylab="Probability", col="blue",
          xlab="Sample Size", type="l", yaxs = "i", xaxs = "i" ,
          lwd=8, main=c("The difference between t and z distributions is", "not driving results: t=blue, z=green"))
lines(n, pz, col="green", lwd=3)
lines(y=c(0.975,0.975), x=c(2,500), col="red", lwd=2)

# Thus, for simplicity, we will use the z distribution.

# Let's look briefly at what happens when we vary effect size.

# Start with a blank graph with a range of n from 2 to 100
n<- 2:100

# we have to redefine the standard error vector
se = s/(n^.5)

plot(n, n*0, ylim=c(.5,1), ylab="Probability", col="blue",
          xlab="Sample Size", type="n", yaxs = "i", xaxs = "i" ,
          lwd=2, main="Probability of Rejecting the Null")

# Now let's loop through the difference in power as a function of effect size:
for (d in seq(.1,1,.1)) {
  z = d/se
  pz <- pnorm(z)
  lines(n, pz, col=gray(d), lwd=2)
}  
lines(y=c(0.975,0.975), x=c(2,200), col="orange", lwd=2)

# I will now strengthen the axes since they seem to have been somewhat overwritten by the other lines on the graph.
lines(y=c(1,1), x=c(2,200), lwd=2)      
lines(y=c(0,1), x=c(2,2), lwd=2)
# Thus we can see that if the effect is close to one standard deviation, then a very small sample size is all that is needed to identify the effect.

# Overall the construction of power analysis is found by rearranging the common form of the t-test.

# t = (M0 - M1)/(SE)
# SE = s/n^.5
# t=n^.5 * (M0-M1)/s
# (t*s/(M0-M1))^2 = n
# d=M0-M1
# n = t^2 (var/d^2)

# In order to find t you find what t value is required in order to reject at the desired level.

# T~1.96 at 5% rejection rate.

# Let's imagine our previous example
d=.1
s=1

# Let's see how our required sample size increases as our desired rejection rate also increases.
alpha = seq(.5,.975,.025)

target.z <- qnorm(alpha)

nreq = (target.z)^2*s^2/d^2

plot(nreq, alpha, type="l", main="As the required rejection rate increases, patient count increases")

# Thus on a target difference in means of .1 standard deviation, there need be at least:
qnorm(.975)^2*s^2/d^2

# Around 384 observations.  (This is of course because the effect size is very small.)

# Statisticians often do not stop there.  They often also assume that the probabiliy of failing to reject the null when there actually is an effect (Type II error) is independent of Type I errors and thus need to be taken into account.  This is done by saying that we are willing to accept a 20% possibility that there is an effect but we are not detecting it.  We can adjust the requirement simply by adding an additional z score:
(qnorm(.975)+qnorm(.8))^2*s^2/d^2

# This nearly doubles the required sample.

# At this point I would like to make the caveat that most of my posts I simply pull out of my head.  Thus, if there is something wrong I welcome corrections and apoligize.  I say this now because I am generally unfamiliar with power analysis and am just piecing together this simulation as a learning device.