How to write a function that rejects a Z-score in R?

Question

How can I write a function that returns "reject" when the the zscore>=qnorm(1-alpha/2) for 10 simulations for alpha=0.05 and a sample size 10. I wrote the following code but I'm not getting the corret output. "zscore" is the test statistic and t is normal with mean and standard deviation 6/n. sims corresponds to number of simulations to perform.This function should imitate a Monte Carlo assessment.

testsk=function(n,alph,sims){
    t=numeric(sims)
for (i in 1:sims) {
    x=rnorm(n)
t[i]=skewness(x)
zscore=t/(6/n)
return(zscore)
}
if(zscore>=qnorm(1-alph/2)){
print("REJECT")
}
}


testsk(10,0.05,10)

Thanks!

Answer 1

Following your edit I believe what you want to do is to see how many times over sims trials the skewness calculated from a size n sample taken from a normal distribution would be rejected as too skewed by a significance test at the alph level.

You have several coding problems

1. You want to do the z score test for every trial, but the test is outside the loop.
2. The z score is calculated using the vector t, but you want to calculate it using the scalar t[i].

3. There is a return statement inside the loop which will cause the function to terminate in the first iteration of the loop, returning the z score. For reason no 2., the z score is a vector, but its second to last values are all zero because you've only run one iteration, so a typical output from the function is as follows

   0.003623371 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000

Fixing these immediate problems results in the following code

library(e1071)
testsk=function(n,alph,sims) {
    t=numeric(sims)
    for (i in 1:sims) {
        x=rnorm(n)
        t[i]=skewness(x)
        zscore=t[i]/(6/n)
        if(zscore>=qnorm(1-alph/2)){
            print("REJECT")
        }
    }
}

However this stil suffers from some problems:

1. From a programming point of view
   • printing out "REJECT" gives immediate feedback but isn't very scalable. If you had sims=1000 you would be better off returning the number of rejects, nr. And if you still wanted to print out "REJECT" nr times you could do so :)
   • Also the code could be simpler, and written in more of an R style, vectorizing rather than using loops. This would also have the advantage of being much faster. Because R is an interpreted language, vectorizing makes a huge difference because the number crunching can happen under the hood without R having to step through your for loop over and over.
2. Perhaps more seriously, there are some statistical problems:
   • 6/n is an estimate of the variance of the skewness (wikipedia), but you want the standard deviation so you need to take the square root of 6/n.
   • The code rejects if the z score is greater than the 1-alph/2th quantile but it should also reject if the z score is less than the alph/2th quantile. As it stands your rejection region is alph/2 not alph.
   • There may be other issues too but these are the main ones in my opinion. (I assume you're aware that 6/n is only a good estimate of the variance for large samples.)

A program that is along the right lines is as follows

library(e1071)
testsk=function(n,alph,sims) {
    # Generate random numbers in a matrix, each trial is a row
    X=matrix(rnorm(sims*n), ncol=n)
    # Get skewnesses, 1 means apply to rows
    skews=apply(X,1,skewness)
    # Calculate z score vector and rejection vector
    zscore=skews/sqrt(6/n)
    reject=!(qnorm(alph/2) < zscore & zscore < qnorm(1-alph/2))
    # Return the number of rejections
    sum(reject)
}

You ought to be able to modify it to suit your purposes, but I can clarify if necessary.

Answer 2

Not sure what are you trying to achieve but here is one way you can do this

testsk <- function(n, alph, sims){
  for (i in 1:sims){
  x <- rnorm(n)
  zscore <- skewness(x)/(6/n)
  cat(paste0("Simulation #", i,":"), ifelse(zscore >= qnorm(1 - alph/2), "REJECT", "Don't REJECT"), "\n")
  }
}

n <- 10
alph <- .05
sims <- 10
testsk(n, alph, sims)

#Simulation #1: Don't REJECT 
#Simulation #2: REJECT 
#Simulation #3: Don't REJECT 
#Simulation #4: Don't REJECT 
#Simulation #5: Don't REJECT 
#Simulation #6: Don't REJECT 
#Simulation #7: Don't REJECT 
#Simulation #8: Don't REJECT 
#Simulation #9: Don't REJECT 
#Simulation #10: Don't REJECT 

Answer 3

You are incorrectly looping over sims. Can you explain what you are trying to do?

testsk <- function(n,alph,sims) {
  t <- numeric(sims)
  for (i in seq_along(sims)) {
    x <- rnorm(n)
    t[i] <- skewness(x)
  }
  zscore <- t/(6/n)
  if (any(zscore>=qnorm(1-alph/2))) {
    print("REJECT")
  }
  return(zscore)
}

testsk(10,0.05,10)