danielcfurr / edstan Goto Github PK

The stan program should use the Cholesky decomposition, perhaps with a weak prior on the covariance matrix. Also, a better example is needed for the help documentation--the four-item spelling data is not a good example.

Dear readers,
I just worked through your case study "Two-Parameter Logistic Item Response Model" (https://mc-stan.org/users/documentation/case-studies/tutorial_twopl.html) and ran into an issue adapting your odds ratio function. I had some compares where some of the pairs counted to zero and table() gave back a frequency table smaller that 2x2. It ttok my some hours of work to come up with this (probably not optimized) code snipped (see below).

Maybe you know a better way to handle this issue but i think it might be handy for future readers of your case study to bypass this error faster.

Sincerelly Simon

  n <- nrow(X)
  p <- ncol(X)
  ind <- t(combn(p, 2))
  nind <- nrow(ind)
  OR <- numeric(nind)
  for (i in 1:nind) {
    # print(i)
    xtab <- my.freq.table(X[, ind[i, 1]], X[, ind[i, 2]]) 
    n00 <- xtab[1, 1]
    n01 <- xtab[1, 2]
    n10 <- xtab[2, 1]
    n11 <- xtab[2, 2]
    OR[i] <- (n00 * n11)/(n01 * n10)
  }
  pw.OR <- data.frame(ind, OR)
  names(pw.OR) <- c("item_i", "item_j", "OR")
  
  return(pw.OR)
}

my.freq.table <- function(x,y) {
  xtab <- table(x,y)
  
# handles frequency tables with lower rank
  if((nrow(xtab) != ncol(xtab)) | (nrow(xtab) == 1 & nrow(xtab) == 1)) {
    xtab <- table(x,y) %>% as.data.frame() %>% defactor()
    dat <- data.frame(c(0,0,0),
                      c(1,0,0),
                      c(0,1,0),
                      c(1,1,0)) %>% t() %>% data.frame() %>% setNames(c('x','y','Freq'))
    rownames(dat) <- NULL
    
    dat <- dat %>% mutate(Freq = ifelse(x == xtab$x & y == xtab$y, xtab$Freq, Freq))
    
    xtab <- xtabs(Freq ~ x+y, data = dat)
  } 
  
  return(xtab)
}

# I don't can handle factors well yet so I just pull the numbers out
defactor <- function(d) {
  dat <- data.frame(matrix(NA, nrow = nrow(d), ncol = ncol(d))) %>% setNames(colnames(d))
  
  for (i in 1:nrow(d)) {
    for (j in 1:ncol(d)) {
      if (is.factor(d[i,j])) {dat[i,j] <- as.numeric(levels(d[i,j]))[d[i,j]]}
      else {dat[i,j] <- d[i,j]}
    }
  }
  
  return(dat)
}

The output for the 2PL results() method is using odd choices for number of decimal places. Also, n_eff should be rounded to whole numbers and rhat should be set to 1 if < 1.

In the documentation file for 2PL regression, the formula used for the logistic function is:

However, the formula in the stan file is:
y ~ bernoulli_logit(alpha[ii].*theta[jj] - beta[ii]);

Which, if I am understanding the stan syntax correctly and considering matrix multiplication between alpha and theta, is:

Could you verify if this is correct?

Either (1) integrate the latent reg code into edstan by generalizing it to multiple person covariates and then writing a wrapper for it, or (2) store the code outside of edstan.

Add "probs" option to all $show() methods. This is just like the probs option in print() for stanfit objects.

First of all, thank you very much for great package ! Incredibly useful. I was wondering if there are any functions to generate trace plots or information plots for item sets (such as the ones provided in other packages like mirt). I understand with edstan we would have entire posterior distributions for parameters instead of point estimates. Can you provide some hints or divert me to any existing resources that would make such plots from edstan fits easier? Thanks a lot.

Methods $best_logLik() and $DIC() both rely on getting likelihood for given parameter values. Verify that the method used to do this is correct. Some built in rstan functions may provide this functionality, but are not well-documented. May have to check with Columbia team.

For the 2PL, write a method to create a data frame that provides all the information in print or $results() but as a data frame. Should include every value of theta I think -- this is equivalent to EAP estimates.