More Explorations with catR

# For the purposes of simulating computerized adaptive tests
# the R package catR is unparallelled.
 
# catR is an excellent tool for students who are curious about
# how a computerized adaptive test might work. It is also useful
# for testing companies that are interested in seeing how
# their choices of number of items, or model, stopping rule,
# or quite a few of the other options which are available
# when designing a specific computerized adaptive test.
 
# In this post I will explore some of the features of the 
# function randomCAT, an extremely powerful function 
# that simulates an entire response pattern for an individual.
 
# In a previous post I explore  some of the other function 

# in catR in order to step by step demonstrate how to use 

# the package to simulate a test.
 
library("catR")
 
# First let's generate an item bank. 
# Items specifies how many items to generate
 
# Model specifies which model to use in generating the items
# a,b,c Priors are specifying distributions to draw
# the parameters from for each item.
 
# The final set of arguments is for specifying
# what range of theta values the bank will initially
# draw item parameters for.  Theta values are the typical
# latent traits for which item response theory is concerned
# with estimating.
Bank <- createItemBank(items = 500, model = "3PL", 
                       aPrior=c("norm",1,0.2), 
                       bPrior=c("norm",0,1), 
                       cPrior=c("unif",0,0.25),
                       thMin = -4, thMax = 4,
                       step = 0.05)
 
# We may want to examine the object we have created called "Bank"
attributes(Bank)
 
# Within the Bank object of class "itBank" there is three named
# attributes.
 
# itemPar lists the item parameters for those items which have been
# generated.  We could see a histogram of difficulty parameters (b) by
# targeting within the Bank object:
 
hist(Bank$itemPar[,2], breaks=30, 
     main="Distribution of Item Difficulties",
     xlab="b parameter")


# We can also see how much information a particular item would add
# accross a range of ability levels.  This information is already
# available within the Bank object under the names infoTab and
# theta.
 
# Plot the first item's information
plot(rep(Bank$theta,1),Bank$infoTab[,1], 
     type="l", main="Item 1's information", 
     xlab="Ability (theta)", ylab="Information")
 
# Plot the first 3 items
# By specifying type = "n" this plot is left empty
nitems = 3
plot(rep(Bank$theta,nitems),Bank$infoTab[,1:nitems], type="n", 
     main=paste0("First ",nitems," items' information"), 
     xlab="Ability (theta)", ylab="Information")
# Now we plot the 
for (i in 1:nitems) lines(Bank$theta,Bank$infoTab[,i], 
                          col=grey(.8*i/nitems))


# We can see how different items can have information that 
# spans different ability estimates as well as some items
# which just have more information than other items.
 
 
# Plotting all 500 items (same code as previously but now
# we specify the number of items as 500)
nitems = 500
plot(rep(Bank$theta,nitems),Bank$infoTab[,1:nitems], type="n", 
    main=paste0("First ",nitems," items' information"), 
    xlab="Ability (theta)", ylab="Information")
for (i in 1:nitems) lines(Bank$theta,Bank$infoTab[,i], 
                          col=grey(.8*i/nitems)) 
# This plot may look nonsensical at first.  Be it actually
# provides some useful information.  From it you can see the
# maximum amount of information available for any one
# item at different levels of ability.  In the places where
# there is only one very tall item standing out we may be 
# concerned about item exposure since subjects which seem to
# be in the area of that item are disproportionately more likely
# to get the same high info item than other other subjects
# in which the next highest item is very close in information
# to the max item.
 
# To see the max information for each ability we can add a line.
lines(Bank$theta,apply(Bank$infoTab, 1, max), col="blue", lwd=2)
 
# We might also be interested in seeing how much information
# on average a random item chosen from the bank would provide
# or in other words what is the expected information from a
# random item drawn from the bank at different ability levels.
lines(Bank$theta,apply(Bank$infoTab, 1, mean), col="red", lwd=2)
 
# Or perhaps we might want to see what the maximum average information
# for a 20 item test might be. So we calculate the average information
# for the top 20 items at different ability levels.
maxmean <- function(x, length=20) mean(sort(x, decreasing=T)[1:length])
maxmean(1:100) # Returns 90.5, seems to be working properly
 
lines(Bank$theta,apply(Bank$infoTab, 1, maxmean), col="orange", lwd=3)
 


# Now this last line is very interesting because it reflects
# per item the maximum amount of information this bank can provide
# given a fixed length of 20. Multiply this curve by 20 and it will give
# us the maximum information this bank can provide given a 20 item test
# and a subject's ability.
 
# This can really be thought of as a theoretical maximum for which
# any particular CAT test might attempt to meet but on average will
# always fall short.
 
# We can add a lengend
legend(-4.2, .55, c("max item info", "mean(info)", 
                    "mean(top items)"),
       lty = 1, col = c("blue","red","orange"),  adj = c(0, 0.6))
 
 
 
library("reshape")
library("ggplot2")
 
# Let's seperate info tab
infoTab <- Bank$infoTab
 
# Let's add three columns to info tab for max, mean, and mean(top 20)
infoTab <- cbind(infoTab, 
                 apply(Bank$infoTab, 1, max),
                 apply(Bank$infoTab, 1, mean),
                 apply(Bank$infoTab, 1, maxmean))
 
 
# Melt will turn the item information array into a long object
items.long <- melt(infoTab)
 
# Let's assign values to the first column which are thetas
items.long[,1] <- Bank$theta
 
# Now we are ready to name the different columns created by melt
names(items.long) <- c("theta", "item", "info")
 
itemtype <- factor("Item", c("Item","Max", "Mean", "Mean(Max)"))
items.long <- cbind(items.long, type=itemtype)
items.long[items.long$item==501,4] <- "Max"
items.long[items.long$item==502,4] <- "Mean"
items.long[items.long$item==503,4] <- "Mean(Max)"
 
# Now we are ready to start plotting
# Assign the data to a ggplot object
a <- ggplot(items.long, aes(x=theta, y=info, group=item))
 
# Plot a particular instance of the object
a + geom_line(colour = gray(.2)) +
    geom_line(aes(colour = type), size=2 , 
            subset = .(type %in% c("Max", "Mean", "Mean(Max)")))
 


# Now let's look at how the randomCAT function works.
# There are a number of arguments that the randomCAt function
# can take.  They can be defined as lists which are fed
# into the function.
 
# I will specify only that the stoping rule is 20 items.
# By specifying true Theta that is telling random CAT what the
# true ability level we are estimating.
res <- randomCAT(trueTheta = 3, itemBank = Bank, 
                 test=list(method = "ML"),
                 stop = list(rule = "length", thr = 20))
# I specify test (theta estimator) as using ML because the
# default which is Bayesian model is strongly centrally
# biased in this case.
 
# Let's examine what elements are contained with the object "res"
attributes(res)
 
# We can see our example response pattern.
thetaEst <- c(0, res$thetaProv)
 
plot(1:21, thetaEst, type="n",
     xlab="Item Number",
     ylab="Ability Estimate", 
     main="Sample Random Response Pattern")
# Add true ability line
  abline(h=3, col="red", lwd=2, lty=2)
# Add a line connecting responses
  lines(1:21, thetaEst, type="l", col=grey(.8))
# Add the response pattern to 
  text(1:21, thetaEst, c(res$pattern, "X"))
# Add the legend
  legend(15,1,"True Ability", col="red", lty=2, lwd=2)
 


# Plot the sample item information from the set of items selected.
plot(rep(Bank$theta,20),Bank$infoTab[,res$testItems], type="n", 
     main="High information items are often selected", 
     xlab="Ability (theta)", ylab="Information")
for (i in 1:500) lines(Bank$theta,Bank$infoTab[,i], col=grey(.75)) 
# Now we plot the 
for (i in res$testItems) lines(Bank$theta,Bank$infoTab[,i], 
                               lwd=2, col=grey(.2))
 


# Now let's see how randomCat performs with a random draw
# of 150 people with different ability estimates.
 
npers <- 150   # Specify number of people to simulate
 
theta <- rnorm(npers)
# Draw a theta ability level vector
 
thetaest <- numeric(npers)
# Creates an empty vector of zeros to hold future estimates
# of theta
 
# Create an empty item object
items.used <- NULL
 
# Create an empty object to hold b values for items used
b.values <- NULL
 
 
for (i in 1:npers) {
  # Input the particular theta[i] ability for a particular run.
  res <- randomCAT(trueTheta = theta[i], 
                   itemBank = Bank, 
                   test=list(method = "ML"),
                   stop = list(rule = "length", thr = 20))
  # Save theta final estimates
  thetaest[i] <- res$thFinal
  # Save a list of items selected in each row of items.used
  items.used <- rbind(items.used, res$testItems)
  # Save a list of b values of items selected in each row
  b.values <- rbind(b.values, res$itemPar[,2])
 
}
 
# Let's see how our estimated theta's compare with our true
plot(theta, thetaest, 
     main="Ability plotted against ability estimates",
     ylab="theta estimate") 

# To get a sense of how much exposure our items get 
itemTab <- table(items.used)
 
length(itemTab)
# We can see we only have 92 items used for all 150 subjects
# taking the cat exam.
 
mean(itemTab)
# On average each item used is exposed 32 times which means
mean(itemTab)/150
# over a 20% exposure rate on average in addition to some items
# have much higher exposure rates.

Created by Pretty R at inside-R.org

Results of an Informal Survey of R users

# This post does some basic correlation analysis between responses
# to the survey I recently released through R Shiny at:
http://www.econometricsbysimulation.com/2013/11/Shiny-Survey-Tool.html
 
# I have saved the data from the survey after completing the survey myself.
# This data is incomplete because the survey has been running since I
# saved the survey and because shinyApp.io server automatically resets
# data every so often. Thus survey results are lost sometimes.
# (Something to be remedied at a future time.)
 
# Let's get our data
Rsurvey <- read.csv(paste0(
  "https://raw.github.com/EconometricsBySimulation/",
  "Shiny-Demos/master/Survey/sample-results.csv"))
 
summary(Rsurvey)
# Looking at the data we have 321 responses in total though
# on average it looks like we have closer to 250 responses
# to work with.
 
# The majority of respondents consider their knowledge of R
# to be either Advanced (119) or Moderate (92).
 
# The majority of users are aged 26-35 (121) or 36-45 (79).
 
# The frequency that respondents read R bloggers is
# most frequently daily (164) or weekly (75).
 
# The frequency of respondents read my blog never (149)
# or monthly (48).
 
# The self-reported technical knowledge of most users in a
# theoretic stastics/econometrics/pschometrics field was
# most frequently reported as either moderate (103) or
# advanced (97).
 
# The favorite colors of people was most frequently blue (114)
# and green (66).
 
# The vast majority of respondents were male (243) compared
# with only 22 females. Looks like R bloggers is not going
# to become a dating website in the near future.
 
# Of those who selected an area of research the majority
# chose data analysis (150) followed by (statistics).
 
# As for the user's knowledge of Shiny few had much at all with
# Basic (118) and Non (76) being the most frequent responses.
 
# Finally, as to the question of "What is the air speed velocity
# of an unladen swallow" the majority of respondents chose
# the correct response to the Monte Python reference (129)
# while the next largest group indicated that "they did not 
# know" (98).
 
# Well let's see if there is any correlation between particular
# outcomes of interest.
 
# Let's see if there is a correlation between knowledge of R
# and frequency of reading R bloggers.
knowledgeR <- rep(NA, nrow(Rsurvey))
 
knowledgeR[Rsurvey[,2]=="None"] <- 0
knowledgeR[Rsurvey[,2]=="Basic"] <- 1
knowledgeR[Rsurvey[,2]=="Moderate"] <- 2
knowledgeR[Rsurvey[,2]=="Advanced"] <- 3
knowledgeR[Rsurvey[,2]=="Expert"] <- 4
 
frequency <- rep(NA, nrow(Rsurvey))
 
# Convert these rates to number of days per year
# reading R bloggers.
frequency[Rsurvey[,5]=="None"] <- 0
frequency[Rsurvey[,5]=="daily"] <- 360
frequency[Rsurvey[,5]=="weekly"] <- 50
frequency[Rsurvey[,5]=="monthly"] <- 12
 
cor(frequency, knowledgeR, use="pairwise.complete.obs")
# There seems to be modest correlation beteen # of days spent
# reading R bloggers and self assessment of R expertise.
 
# Let's see if coldness or warmth allong the color
# spectrum is a useful variable.
warmth <- rep(NA, nrow(Rsurvey))
 
warmth[Rsurvey[,8]=="blue"] <- 0
warmth[Rsurvey[,8]=="green"] <- 1
warmth[Rsurvey[,8]=="orange"] <- 2
warmth[Rsurvey[,8]=="red"] <- 3
 
cor(warmth, knowledgeR, use="pairwise.complete.obs")
# There is a slight negative correlation between
# the favorite color warmth of users and self-reported
# knoweldge of R.
 
# Finally, let's look at success on the Monte Python
# trivia question.
monte <- rep(NA, nrow(Rsurvey))
 
monte[Rsurvey[,12]=="I don't know!"] <- 0
monte[Rsurvey[,12]=="~50 MPH"] <- 0
monte[Rsurvey[,12]==
  "What do you mean? An African or European swallow?"] <- 1
 
cor(monte, knowledgeR, use="pairwise.complete.obs")
# We see a modest positive correlation between
# knowledge of R and being able to answer Monte Python
# trivia.
 
knowledgeStats <- rep(NA, nrow(Rsurvey))
 
knowledgeStats[Rsurvey[,6]=="None"] <- 0
knowledgeStats[Rsurvey[,6]=="Basic"] <- 1
knowledgeStats[Rsurvey[,6]=="Moderate"] <- 2
knowledgeStats[Rsurvey[,6]=="Advanced"] <- 3
knowledgeStats[Rsurvey[,6]=="Expert"] <- 4
 
cor(knowledgeStats, knowledgeR, use="pairwise.complete.obs")
# There seems to be a very stong correlation with
# self reported knowledge of Statstics and knowedge
# of R.
 
# In order to see the data points a little more clearly
# I will add some tiny noise to both knowledge sets
knowledgeStatsN <- knowledgeStats + .15*rnorm(nrow(Rsurvey))
knowledgeRN <- knowledgeR + .15*rnorm(nrow(Rsurvey))
 
plot(knowledgeRN, knowledgeStatsN, 
    main="Knowledge of Stats against that of R",
    xlab="R", ylab="Stats")
# We can see the diagnol elements 2,2 and 3,3 have the most
# frequency. 

  
summary(lm(knowledgeR~warmth+frequency+monte+knowledgeStats))
# Overall, we can see that knowledge of Statistics
# stongly positively predicts knowledge of R.
# Frequency of reading R bloggers also seems to
# have an effect size significant nearly at the
# %5 level.
 
# The coefficient on frequency is very small but that is because
# the scale is quite large (from 0 to 360).  However
# if someone where to start reading R bloggers daily
# (assuming R-bloggers -> R knowledge one directionally)
# Then we would expect a change of R knowledge of:
.0006994*360
# 0.251784
 
# Not as large a predictor as knowledge of stats but certainy
# there exists some relationship.
 
# Of course little causal relationships can be inferred from the
# data.  We cannot expect reading R bloggers to be independent
# of self-assessed knowledge of R any more than knowledge of 
# statistics to be uncorrelated with knowledge of R since
# many users, learn R simultaneously with statistics.

Created by Pretty R at inside-R.org

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?

A Survey Tool Designed Entirely in Shiny Surveying Users of R

I have written a very basic survey tool built entirely in the Shiny package of R.  I hope the tool is useful.  Modifying the survey for your own purposes is trivially easy (I hope).


I have not commented my code so it is pretty messy right now.  You can find the source on GitHub.

In order to run your own survey all you need do is edit the Qlist.csv file.  It contains one row for each question and one column for the question text and additional columns for potential answers (which will enter the radio buttons).  All of the questions will have the default option of "Prefer not to answer" selected.  It should be easy to modify this option as well since you need only find where this text appears in the server.R code.

It should be equally easy to change the introductory message and the accreditation to me, though it would be nice if you left a link somewhere back to me.

I will comment the GitHub files as soon as I have chance.  If people find this tool useful perhaps I will work on improving it.  Please tell me if you find it useful or have suggestions for future improvements!

Please take note that it is really not up in running in a sustainable way since I was recently told that sometimes the ShinyApps.io server must be restarted which means that all of the files are restored to those when the app was originally set up.


http://econometricsbysimulation.shinyapps.io/Survey/