R/sim.item.R

In psych: Procedures for Psychological, Psychometric, and Personality Research

"sim.item" <-
function (nvar = 72 ,nsub = 500, 
    circum = FALSE, xloading =.6, yloading = .6, gloading=0, xbias=0,  ybias = 0,categorical=FALSE, low=-3,high=3,truncate=FALSE,threshold=NULL) 
	{ 
	avloading <- (xloading+yloading)/2
	
	errorweight <- sqrt(1-(avloading^2  + gloading^2))  #squared errors and true score weights add to 1
    g <- rnorm(nsub) 
	truex <- rnorm(nsub)* xloading  +xbias #generate normal true scores for x + xbias
	truey <- rnorm(nsub) * yloading + ybias #generate normal true scores for y + ybias
	 if(!is.null(threshold)) {if (length(threshold) < nvar) threshold <- sample(threshold, nvar, replace=TRUE)}
   

	if (circum)  #make a vector of radians (the whole way around the circle) if circumplex
	{radia <- seq(0,2*pi,len=nvar+1)  
      rad <- radia[which(radia<2*pi)]        #get rid of the last one
     } else   rad <- c(rep(0,nvar/4),rep(pi/2,nvar/4),rep(pi,nvar/4),rep(3*pi/2,nvar/4))  #simple structure 
        
	error<- matrix(rnorm(nsub*(nvar)),nsub)    #create normal error scores

	#true score matrix for each item reflects structure in radians
	trueitem <- outer(truex, cos(rad)) + outer(truey,sin(rad)) 

	item<- gloading * g +  trueitem  + errorweight*error   #observed item = true score + error score 
    if (categorical) {
        if(is.null(threshold)){
    	item = round(item)       #round all items to nearest integer value
		item[(item<= low)] <- low     
		item[(item>high) ] <- high  } else {
		  i <- 1:nvar
		  item <- t(t(item [,i]) > threshold[i]) + 0  }  
		}
	#if (truncate) {item[item < cutpoint] <- 0  }  #replaced with threshold
	colnames(item) <- paste("V",1:nvar,sep="")
	return (item) 
	}  

"sim.spherical" <-
function (simple=FALSE, nx=7,ny=12 ,nsub = 500,  xloading =.55, yloading = .55, zloading=.55, gloading=0, xbias=0,  ybias = 0, zbias=0,categorical=FALSE, low=-3,high=3,truncate=FALSE,threshold=NULL) 
	{ 

	nvar <- nx * (ny-1)
	 if(!is.null(threshold)) {if (length(threshold) < nvar) threshold <- sample(threshold, nvar, replace=TRUE)}
	errorweight <- sqrt(1-(xloading^2  + yloading^2 + zloading^2+  gloading^2))  #squared errors and true score weights add to 1
    g <- rnorm(nsub) 
	truex <- rnorm(nsub)  +xbias #generate normal true scores for x + xbias
	truey <- rnorm(nsub)  + ybias #generate normal true scores for y + ybias
	truez <- rnorm(nsub)  + zbias #generate normal true scores for y + ybias
    true <- matrix(c(g,truex,truey,truez),nsub)
    
	 #make a vector of radians (the whole way around the circle) if circumplex
	if(!simple) {f1  <- rep(cos(seq(0,pi,length.out = nx)),each=ny-1)*xloading
	 f2 <- rep(cos(seq(0,2*pi*(ny-1)/ny,length.out=ny-1)),nx)*yloading
    f3 <- rep(sin(seq(0,2*pi*(ny-1)/ny,length.out=ny-1)),nx)*zloading
    f2 <- f2 * (1-f1^2)  #added to make the cylinder a sphere
    f3 <- f3 * (1-f1^2)
    f <- matrix(c(rep(gloading,(ny-1)*nx),f1,f2,f3),(ny-1)*nx)
      
     } else {    #simple structure -- 
     nvar <- 4*ny
     f1  <- rep(c(1,0,-1,0),each=ny)
	 f2 <- rep(c(1,0,-1,0),ny)
     f3 <- rep(c(0,1,0,-1),ny)
    f <- matrix(c(rep(gloading,ny*4),f1,f2,f3),4*ny)
     }
        
	error<- matrix(rnorm(nsub*(nvar)),nsub)    #create normal error scores

	#true score matrix for each item reflects structure in radians
	

	item <- true %*% t(f) + errorweight*error   #observed item = true score + error score 
    if (categorical) {
     if(is.null(threshold)){
        
    	item = round(item)       #round all items to nearest integer value
		item[(item<= low)] <- low     
		item[(item>high) ] <- high   
		} else {
		  i <- 1:nvar
		  item <- t(t(item [,i]) > threshold[i]) + 0  }  
		  }
	#if (truncate) {item[item < cutpoint] <- 0  }
	colnames(item) <- paste("V",1:nvar,sep="")
	return (item) 
	}  



"sim.rasch" <-
function (nvar = 5 , n = 500, low=-3,high=3,d=NULL, a=1,mu=0,sd=1) 
	{ 
	if(is.null(d)) {d <- seq(low,high,(high-low)/(nvar-1))}
	theta <- rnorm(n,mu,sd)
	item <- matrix(t(1/(1+exp(a*t(-theta %+% t( d))))),n,nvar)
	#now convert these probabilities to outcomes
    item[] <- rbinom(n*nvar, 1, item)  
   colnames(item) <- paste("V",1:nvar,sep="")
   result <- list(items=item,tau=d,theta=theta)
	return (result) 
	}  
	
	
"sim.irt" <- 
function (nvar = 5 ,n = 500,low=-3,high=3,a=NULL,c=0,z=1,d=NULL, mu=0,sd=1,mod="logistic",theta=NULL) 
	{ 
	if(mod=="logistic") {result <- sim.npl(nvar,n,low,high,a,c,z,d,mu,sd,theta)} else {result <- sim.npn(nvar,n,low,high,a,c,z,d,mu,sd,theta)}
	return (result) 
	}  

"sim.npn" <- 
function (nvar = 5 ,n = 500, low=-3,high=3,a=NULL,c=0,z=1,d=NULL,mu=0,sd=1,theta=NULL) 
	{ 
	if(is.null(d)) {d <- seq(low,high,(high-low)/(nvar-1))} else {if(length(d)==1) d <- rep(d,nvar)}
	if(is.null(a)) {a <- rep(1,nvar)}
	if(is.null(theta)) {theta <- rnorm(n,mu,sd)} # the latent variable
	
	item <- matrix(t(c+(z-c)*pnorm(a*t(theta %+% t(- d)))),n,nvar)  #need to transpose and retranpose to get it right
	#now convert these probabilities to outcomes
	
    item[] <- rbinom(n*nvar, 1, item)
	colnames(item) <- paste("V",1:nvar,sep="")        
   result <- list(items=item,discrimination=a,difficulty=d,gamma=c,zeta=z,theta=theta)
	return (result) 
	}  
	
	
"sim.npl" <- 
function (nvar = 5 ,n = 500, low=-3,high=3,a=NULL,c=0,z=1,d=NULL,mu=0,sd=1,theta=NULL) 
	{ 
	if(is.null(d)) {d <- seq(low,high,(high-low)/(nvar-1))} else {if(length(d)==1) d <- rep(d,nvar)}
	if(is.null(a)) {a <- rep(1,nvar)}
	if(is.null(theta)) {theta <- rnorm(n,mu,sd)} 
	item <- matrix(t(c+(z-c)/(1+exp(a*t((-theta %+% t( d)))))),n,nvar)
    item[] <- rbinom(n*nvar, 1, item) #now convert these probabilities to outcomes	  
    colnames(item) <- paste("V",1:nvar,sep="")
    result <- list(items=item,discrimination=a,difficulty=d,gamma=c,zeta=z,theta=theta)
	return (result) 
	} 


"sim.poly" <- 
function (nvar = 5 ,n = 500,low=-2,high=2,a=NULL,c=0,z=1,d=NULL, mu=0,sd=1,cat=5,mod="logistic",theta=NULL) 
	{ 
	if(mod=="normal") {result <- sim.poly.npn(nvar,n,low,high,a,c,z,d,mu,sd,cat,theta)} else {result <- sim.poly.npl(nvar,n,low,high,a,c,z,d,mu,sd,cat,theta)}
	return (result) 
	}  	
	
"sim.poly.npn" <- 
function (nvar = 5 ,n = 500, low=-2,high=2,a=NULL,c=0,z=1,d=NULL,mu=0,sd=1,cat=5,theta=NULL) 
	{ cat <- cat - 1
	if(is.null(d)) {d <- seq(low,high,(high-low)/(nvar-1))} else {if(length(d)==1) d <- rep(d,nvar)}
	if(is.null(a)) {a <- rep(1,nvar)}
	if(is.null(theta)) {theta <- rnorm(n,mu,sd)} # the latent variable
	
	item <- matrix(t(c+(z-c)*pnorm(a*t(theta %+% t(- d)))),n,nvar)  #need to transpose and retranpose to get it right
	#now convert these probabilities to outcomes
	item[] <- rbinom(n*nvar, cat, item)
   
	colnames(item) <- paste("V",1:nvar,sep="")        
   result <- list(items=item,discrimination=a,difficulty=d,gamma=c,zeta=z,theta=theta)
	return (result) 
	}  	
	
"sim.poly.npl" <- 
function (nvar = 5 ,n = 500, low=-2,high=2,a=NULL,c=0,z=1,d=NULL,mu=0,sd=1,cat=5,theta=NULL) 
	{ cat <- cat - 1
	if(is.null(d)) {d <- seq(low,high,(high-low)/(nvar-1))} else {if(length(d)==1) d <- rep(d,nvar)}
	if(is.null(a)) {a <- rep(1,nvar)}
	if(is.null(theta)) {theta <- rnorm(n,mu,sd)}
	item <- matrix(t(c+(z-c)/(1+exp(a*t((-theta %+% t( d)))))),n,nvar)
	item[] <- rbinom(n*nvar, cat, item)

      
    colnames(item) <- paste("V",1:nvar,sep="")
     
    result <- list(items=item,discrimination=a,difficulty=d,gamma=c,zeta=z,theta=theta)
	return (result) 
	} 
	
"sim.poly.ideal" <- 
function (nvar = 5 ,n = 500,low=-2,high=2,a=NULL,c=0,z=1,d=NULL, mu=0,sd=1,cat=5,mod="logistic") 
	{ 
	if(mod=="normal") {result <- sim.poly.ideal.npn(nvar,n,low,high,a,c,z,d,mu,sd,cat)} else {result <- sim.poly.ideal.npl(nvar,n,low,high,a,c,z,d,mu,sd,cat)}
	return (result) 
	}  	

	
"sim.poly.ideal.npl.absolute" <- 
function (nvar = 5 ,n = 500, low=-2,high=2,a=NULL,c=0,z=1,d=NULL,mu=0,sd=1,cat=5) 
	{ cat <- cat -1  #binomial is based upon one fewer than categories
	if(is.null(d)) {d <- seq(low,high,(high-low)/(nvar-1))} else {if(length(d)==1) d <- rep(d,nvar)}
	if(is.null(a)) {a <- rep(1,nvar)}
	theta <- rnorm(n,mu,sd)
	item <- matrix(t(c+(z-c)/(1+exp(a*t(abs(-theta %+% t( d)))))),n,nvar)
	item[] <- rbinom(n*nvar, cat, item)

      
    colnames(item) <- paste("V",1:nvar,sep="")
     
    result <- list(items=item,discrimination=a,difficulty=d,gamma=c,zeta=z,theta=theta)
	return (result) 
	} 
	
"sim.poly.ideal.npl" <- 
function (nvar = 5 ,n = 500, low=-2,high=2,a=NULL,c=0,z=1,d=NULL,mu=0,sd=1,cat=5,theta=NULL) 
	{ cat <- cat -1  #binomial is based upon one fewer than categories
	if(is.null(d)) {d <- seq(low,high,(high-low)/(nvar-1))} else {if(length(d)==1) d <- rep(d,nvar)}
	if(is.null(a)) {a <- rep(1,nvar)}
	if(is.null(theta)) {theta <- rnorm(n,mu,sd)}
	item <- 2*matrix((c+(z-c)*exp(a*t(-theta %+% t( d)))/(1+2*exp(a*t(-theta %+% t( d))) + exp(a*t(2*(-theta %+% t( d)))))),n,nvar,byrow=TRUE)
    p <- item
	item[] <- rbinom(n*nvar, cat, item)

      
    colnames(item) <- paste("V",1:nvar,sep="")
    result <- list(p=p,items=item,discrimination=a,difficulty=d,gamma=c,zeta=z,theta=theta)
	return (result) 
	} 
	


"sim.poly.ideal.npn" <- 
function (nvar = 5 ,n = 500, low=-2,high=2,a=NULL,c=0,z=1,d=NULL,mu=0,sd=1,cat=5) {
warning("Not ready for prime time")
	 cat <- cat -1  #binomial is based upon one fewer than categories
	if(is.null(d)) {d <- seq(low,high,(high-low)/(nvar-1))} else {if(length(d)==1) d <- rep(d,nvar)}
	if(is.null(a)) {a <- rep(1,nvar)}
	theta <- rnorm(n,mu,sd) # the latent variable
	
	item <- matrix(t(c+(z-c)*pnorm(abs(a*t(theta %+% t(- d))))),n,nvar)  #need to transpose and retranpose to get it right
	#now convert these probabilities to outcomes
	item[] <- rbinom(n*nvar, cat, item)
   
	colnames(item) <- paste("V",1:nvar,sep="")        
   result <- list(items=item,discrimination=a,difficulty=d,gamma=c,zeta=z,theta=theta)
	return (result) 
	}  	

"sim.poly.npl.m" <- 
function (nvar = 5 ,n = 500, low=-2,high=2,a=NULL,c=0,z=1,d=NULL,mu=0,sd=1,cat=5,theta=NULL) 
	{ cat <- cat - 1
	if(is.null(d)) {d <- seq(low,high,(high-low)/(nvar-1))} else {if(length(d)==1) d <- rep(d,nvar)}
	
	if(is.null(a)) {a <- rep(1,nvar)}
	if(is.null(theta)) {theta <- rnorm(n,mu,sd)}
	item <- matrix(t(c+(z-c)/(1+exp(a*t((-theta %+% t( d)))))),n,nvar)
	item[] <- rbinom(n*nvar, cat, item)

      
    colnames(item) <- paste("V",1:nvar,sep="")
     
    result <- list(items=item,discrimination=a,difficulty=d,gamma=c,zeta=z,theta=theta)
	return (result) 
	} 

#simulate a particular polychoric structure R with specified marginals (m)
#Modified December 17,2013 to actually work
"sim.poly.mat" <- function(R,m,n) {
 e <- eigen(R)
 v <- pmax(e$values,0)
 nvar <- ncol(R)
 ncat <- ncol(m) 
X <- matrix(rnorm(nvar*n),n)
X <- t(e$vectors %*% sqrt(diag(v)) %*% t(X))
marg <- m
Y <- matrix(0,ncol=n,nrow=nvar)
for(i in 1:(ncat)) {
Y[t(X) > marg[,i]] <- i }
return(t(Y))
}