Fun with Bordered Cubes

I am interested in generating 3D reasoning items in R. To this end I have adapted some of the awesome functions built in the rgl library to my ends. My new function is 'cube' and it takes position and automatically sizes itself as a 1x1x1 cube though this can be adjusted through the scale argument.

Find the code on github 

# Easy Bordered Cubes in R
 
library('rgl'); library('magrittr')
 
cube <- function(x=0,y=0,z=0, bordered=TRUE, 
                 filled = TRUE, lwd=2, scale=1,
                 fillcol = gray(.95),
                 bordercol ='black', ...) {
 
  mycube <- cube3d()
 
  # Reduce size to unit
  mycube$vb[4,] <- mycube$vb[4,]/scale*2
 
  for (i in 1:length(x)) {
    # Add cube border
    if (bordered) {
      bcube <- mycube
      bcube$material$lwd <- lwd
      bcube$material$front <- 'line'
      bcube$material$back <- 'line'
      bcube %>% translate3d(x[i], y[i], z[i]) %>% shade3d
    }
    # Add cube fill
    if (filled) {
      fcube <- mycube
      fcube$vb[4,] <- fcube$vb[4,]*1.01
      fcube$material$col <- fillcol
      fcube %>% translate3d(x[i], y[i], z[i]) %>% shade3d
    }
  }
}
 
clear3d()
cube(0,0,0)
cube(1,0,0, filled=F)
cube(-1,0,0, bordered=F)
movie3d(spin3d(axis=c(0,0,1), rpm=20), duration=2.95) 

  
# I mapped R using an excel spreadsheet which 
# translated Xs into 2D locations points
clear3d()
y <- c(1,1,1,1,1,1,2,2,2,2,3,3,3,3,4,4,4,4,4,5,5,
       5,5,6,6,6,6,7,7,7,7,7,8,8,8,8,9,9,9,9,10,
       10,10,10,11,11,11,11,11,12,12,12,12,12)
 
x <- c(8,7,6,3,2,1,7,6,3,2,6,5,3,2,6,5,4,3,2,5,4,
       3,2,5,4,3,2,6,5,4,3,2,7,6,3,2,7,6,3,2,7,6,
       3,2,6,5,4,3,2,5,4,3,2,1)
 
z <- cummax(y)*.5
 
length(x)==length(y)
cube(x,y,z)
movie3d(spin3d(axis=c(0,0,1), rpm=20), duration=2.95) 

 # Let's see how sin and cos can work together
z <- seq(0,6,.1)
x <- sin(pi*z)*z
y <- cos(pi*z)*z
 
clear3d()
cube(x,y,z*2, scale=.75)
movie3d(spin3d(axis=c(0,0,1), rpm=20), duration=2.95) 

# Now let's look at some of the prototyping for
# my 3D reasoning items.
clear3d()
vlist <- c(0,0,0)
for (i in 1:15) {
  cube(vlist[1],vlist[2],vlist[3])
  step <- sample(1:3, 1)
  vlist[step] <- vlist[step]+(-1)^rbinom(1,1,.25)
}
rgl.material(shininess=1)
bg3d("white")
clear3d(type = "lights")
 
rgl.viewpoint(theta = 90, phi = 0, fov = 0)
rgl.snapshot("2014-09-09angle1.png") 

rgl.viewpoint(theta = 0, phi = 0, fov = 0)
rgl.snapshot("2014-09-09angle2.png") 

rgl.viewpoint(theta = 0, phi = 90, fov = 0)
rgl.snapshot("2014-09-09angle3.png")
 


rgl.light()
rgl.viewpoint(theta = 180, phi = 20, fov = 60)
rgl.snapshot("2014-09-09cubes3d1.png") 

rgl.viewpoint(theta = -20, phi = -20, fov = 60)
rgl.snapshot("2014-09-09cubes3d2.png")

Created by Pretty R at inside-R.org

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

Modified Bin and Union Method for Item Pool Design

# Reckase (2003) proposes a method for designing an item pool for a computer
# adaptive test that has been known as the bin and union method.  This method
# involves drawing a subject from a distribution of abilities. Then selecting
# the item that maximizes that subject's information from the possible set of 
# all items given a standard CAT proceedure. This is repeated until the test
# reaches the predifined stopping point.
 
# Then then next subject is drawn and a new set of items is drawn.  Items are
# divided into bins such that there is a kind of rounding.  Items which are
# sufficiently close to other items it terms of parameter fit are considered
# the same item and the two sets are unionized together into a larger pool.
 
# As more subjects are added more items are collected though at a decreasing
# rate as fewer new items become neccessary.
 
# In the original paper he uses a fixed length test though in a forthcoming 
# paper he and his student Wei He is also using a variable length test.
 
# I have modified his proceedure slightly in this simulation.  Rather than
# selecting optimal items for each subject based from the continuous pool
# of possible items I have the test look within the already constructed pool
# to see if any items are within bin length of the subject's estimated ability.
 
# If there is no item then I add an item that perfectly matches the subject's
# estimated ability. The reason I prefer this method is that I think it better
# represents the process that a CAT test typically must go through with items
# close to but rarely exactly at the level of the subjects. Thus the information
# for each subject will be slightly less as a result of this modified method
# relative to the original.
 
# As with the new paper this simulation uses a variable length test. My stopping
# rule is simple.  Once the test achieves a sufficiently high level of 
# information, then it stops.
 
# I have constructed this simulation as one with three nested loops.
# Over subjects within the item pool construction.
 
# It simulates the item pool construction a number of times to get the
# average number of items after each subject as well as a histogram
# of average number of items required at each difficulty level.
 
# I have also included a control for item exposure.  This control 
# dicatates that as the acceptable exposure rate is reduced, more items
# will be required since some are too frequently exposed.
 
# Overall this method is seems pretty great to me.  It allows for
# item selection criteria, stopping rules, and exposure controls
# to be easily modified to accomidate most any CAT design.
 
require("catR")
 
# Variable Length Test
 
# The number of times to repeat the simulation
nsim <- 10
 
# The number of subjects to simulate
npop <- 1000
 
# The maximum number of items
max.items <- 5000
 
# Maximum exposure rate of individual item
max.exposure <- .2
 
# Stop the test when information reaches this level
min.information <- 10
 
# How far away will the program reach for a new item (b-b_ideal)
bin.width <- .25
 
expect.a <- 1
 
p <- function(theta, b) exp(theta-b)/(1+exp(theta-b))
info <- function(theta, b, a=expect.a) p(theta,b)*(1-p(theta,b))*a^2
 
info(0,0)
 
# The choose.item funciton takes an input thetahat and searches
# available items to see if any already exist that can be used
# otherwise it finds a new item.
choose.item <- function(thetahat, item.b, items.unavailable, bin.width) {
  # Construct a vector of indexes of available items
  avail.n <- (1:length(item.b))
 
  # Remove any already make unusuable
  if (length(items.unavailable)>0) 
    avail.n <- (1:length(item.b))[-items.unavailable]
 
  # If there are no items available then generate the next item
  # equal to thetaest.
  if (length(avail.n)==0) 
    return(c(next.b=thetahat, next.n=length(item.b)+1))
 
  # Figure out how far each item is from thetahat
  avail.dist <- abs(item.b[avail.n]-thetahat)
 
  # Reorder the n's and dist in terms of proximity
  avail.n <- avail.n[order(avail.dist)]
  avail.dist <- sort(avail.dist)
 
  # If the closest item is within the bin width return it
  if (avail.dist[1]<bin.width) 
    return(c(next.b=item.b[avail.n[1]], next.n=avail.n[1]))
  # Otherwise generate a new item
  if (avail.dist[1]>=bin.width) 
    return(c(next.b=thetahat, next.n=length(item.b)+1))
}
 
# Define the simulation level vectors which will become matrices
Tnitems <- Ttest.length <- Titems.taken.N <- Titem.b<- NULL
 
# Loop through the number of simulations
for (j in 1:nsim) {
 
 
  # Seems to be working well
  choose.item(3, c(0,4,2,2,3.3), NULL, .5)
 
  # This is the initial item pool
  item.b <- 0
 
  # This is the initial number of items taken
  items.taken.N <- rep(0,max.items)
 
  # A vector to record the individual test lengths 
  test.length <- NULL
 
  # Number of total items after each individual
  nitems <- NULL
 
  # Draw theta from a population distribution
  theta.pop <- rnorm(npop)
 
  # Start the individual test
  for (i in 1:npop) {
    # The this person has a theta of:
    theta0 <- theta.pop[i]
 
    # Our initial guess at theta = 0
    thetahat <- 0
 
    print(paste("Subject:", i,"- Item Pool:", length(item.b)))
    response <- items.taken <- NULL
 
    # Remove any items that would have been overexposed
    items.unavailable <- (1:length(item.b))[!(items.taken.N < max.exposure*npop)]
 
    # The initial imformation on each subject is zero
    infosum <- 0
 
    # Loop through each subject
    while(infosum < min.information) {
 
      chooser <- choose.item(thetahat, item.b, items.unavailable, bin.width)
 
      nextitem <- chooser[2]
      nextb <- chooser[1]
        names(nextitem) <- names(nextb) <- NULL
 
      items.unavailable <- c(items.unavailable,nextitem)
      item.b[nextitem] <- nextb
 
      response <- c(response, runif(1)<p(theta0, nextb))
 
      items.taken <- c(items.taken, nextitem)
 
      it <- cbind(1, item.b[items.taken], 0,1)
 
      thetahat <- thetaEst(it, response)
 
      infosum <- infosum+info(theta0, nextb)
    }
 
    # Save individual values
    nitems <- c(nitems, length(item.b))
 
    test.length <- c(test.length, length(response))
 
    items.taken.N[items.taken] <- items.taken.N[items.taken]+1
  }
 
  # Save into matrices the results of each simulation
  Titem.b <- c(Titem.b, sort(item.b))
  Tnitems <- cbind(Tnitems, nitems)
  Ttest.length <- cbind(Ttest.length, test.length)
  Titems.taken.N <- cbind(Titems.taken.N, items.taken.N)
 
}
 
plot(apply(Tnitems, 1, max), type="n",
     xlab = "N subjects", ylab = "N items",
     main = paste(nsim, "Different Simulations"))
for (i in 1:nsim) lines(Tnitems[,i], col=grey(.3+.6*i/nsim))
 

# We can see that the number of items is a function of the number of

# subjects taking the exam.  This relationship becomes relaxed

# when the number of subjects becomes large and the exposure controls

# are removed.

hist(Titem.b, breaks=30)

 
hist(Ttest.length, breaks=20)

Created by Pretty R at inside-R.org