1"sim.item" <- 2function (nvar = 72 ,nsub = 500, 3 circum = FALSE, xloading =.6, yloading = .6, gloading=0, xbias=0, ybias = 0,categorical=FALSE, low=-3,high=3,truncate=FALSE,cutpoint=0) 4 { 5 avloading <- (xloading+yloading)/2 6 7 errorweight <- sqrt(1-(avloading^2 + gloading^2)) #squared errors and true score weights add to 1 8 g <- rnorm(nsub) 9 truex <- rnorm(nsub)* xloading +xbias #generate normal true scores for x + xbias 10 truey <- rnorm(nsub) * yloading + ybias #generate normal true scores for y + ybias 11 12 if (circum) #make a vector of radians (the whole way around the circle) if circumplex 13 {radia <- seq(0,2*pi,len=nvar+1) 14 rad <- radia[which(radia<2*pi)] #get rid of the last one 15 } else rad <- c(rep(0,nvar/4),rep(pi/2,nvar/4),rep(pi,nvar/4),rep(3*pi/2,nvar/4)) #simple structure 16 17 error<- matrix(rnorm(nsub*(nvar)),nsub) #create normal error scores 18 19 #true score matrix for each item reflects structure in radians 20 trueitem <- outer(truex, cos(rad)) + outer(truey,sin(rad)) 21 22 item<- gloading * g + trueitem + errorweight*error #observed item = true score + error score 23 if (categorical) { 24 25 item = round(item) #round all items to nearest integer value 26 item[(item<= low)] <- low 27 item[(item>high) ] <- high 28 } 29 if (truncate) {item[item < cutpoint] <- 0 } 30 colnames(item) <- paste("V",1:nvar,sep="") 31 return (item) 32 } 33 34 35"sim.rasch" <- 36function (nvar = 5 , n = 500, low=-3,high=3,d=NULL, a=1,mu=0,sd=1) 37 { 38 if(is.null(d)) {d <- seq(low,high,(high-low)/(nvar-1))} 39 theta <- rnorm(n,mu,sd) 40 item <- matrix(t(1/(1+exp(a*t(-theta %+% t( d))))),n,nvar) 41 #now convert these probabilities to outcomes 42 item[] <- rbinom(n*nvar, 1, item) 43 colnames(item) <- paste("V",1:nvar,sep="") 44 result <- list(items=item,tau=d,theta=theta) 45 return (result) 46 } 47 48 49"sim.irt" <- 50function (nvar = 5 ,n = 500,low=-3,high=3,a=NULL,c=0,z=1,d=NULL, mu=0,sd=1,mod="logistic") 51 { 52 if(mod=="logistic") {result <- sim.npl(nvar,n,low,high,a,c,z,d,mu,sd)} else {result <- sim.npn(nvar,n,low,high,a,c,z,d,mu,sd)} 53 return (result) 54 } 55 56"sim.npn" <- 57function (nvar = 5 ,n = 500, low=-3,high=3,a=NULL,c=0,z=1,d=NULL,mu=0,sd=1) 58 { 59 if(is.null(d)) {d <- seq(low,high,(high-low)/(nvar-1))} else {if(length(d)==1) d <- rep(d,nvar)} 60 if(is.null(a)) {a <- rep(1,nvar)} 61 theta <- rnorm(n,mu,sd) # the latent variable 62 63 item <- matrix(t(c+(z-c)*pnorm(a*t(theta %+% t(- d)))),n,nvar) #need to transpose and retranpose to get it right 64 #now convert these probabilities to outcomes 65 66 item[] <- rbinom(n*nvar, 1, item) 67 colnames(item) <- paste("V",1:nvar,sep="") 68 result <- list(items=item,discrimination=a,difficulty=d,gamma=c,zeta=z,theta=theta) 69 return (result) 70 } 71 72 73"sim.npl" <- 74function (nvar = 5 ,n = 500, low=-3,high=3,a=NULL,c=0,z=1,d=NULL,mu=0,sd=1) 75 { 76 if(is.null(d)) {d <- seq(low,high,(high-low)/(nvar-1))} else {if(length(d)==1) d <- rep(d,nvar)} 77 if(is.null(a)) {a <- rep(1,nvar)} 78 theta <- rnorm(n,mu,sd) 79 item <- matrix(t(c+(z-c)/(1+exp(a*t((-theta %+% t( d)))))),n,nvar) 80 item[] <- rbinom(n*nvar, 1, item) #now convert these probabilities to outcomes 81 colnames(item) <- paste("V",1:nvar,sep="") 82 result <- list(items=item,discrimination=a,difficulty=d,gamma=c,zeta=z,theta=theta) 83 return (result) 84 } 85 86 87"sim.poly" <- 88function (nvar = 5 ,n = 500,low=-2,high=2,a=NULL,c=0,z=1,d=NULL, mu=0,sd=1,cat=5,mod="logistic") 89 { 90 if(mod=="normal") {result <- sim.poly.npn(nvar,n,low,high,a,c,z,d,mu,sd,cat)} else {result <- sim.poly.npl(nvar,n,low,high,a,c,z,d,mu,sd,cat)} 91 return (result) 92 } 93 94"sim.poly.npn" <- 95function (nvar = 5 ,n = 500, low=-2,high=2,a=NULL,c=0,z=1,d=NULL,mu=0,sd=1,cat=5) 96 { cat <- cat - 1 97 if(is.null(d)) {d <- seq(low,high,(high-low)/(nvar-1))} else {if(length(d)==1) d <- rep(d,nvar)} 98 if(is.null(a)) {a <- rep(1,nvar)} 99 theta <- rnorm(n,mu,sd) # the latent variable 100 101 item <- matrix(t(c+(z-c)*pnorm(a*t(theta %+% t(- d)))),n,nvar) #need to transpose and retranpose to get it right 102 #now convert these probabilities to outcomes 103 item[] <- rbinom(n*nvar, cat, item) 104 105 colnames(item) <- paste("V",1:nvar,sep="") 106 result <- list(items=item,discrimination=a,difficulty=d,gamma=c,zeta=z,theta=theta) 107 return (result) 108 } 109 110"sim.poly.npl" <- 111function (nvar = 5 ,n = 500, low=-2,high=2,a=NULL,c=0,z=1,d=NULL,mu=0,sd=1,cat=5) 112 { cat <- cat - 1 113 if(is.null(d)) {d <- seq(low,high,(high-low)/(nvar-1))} else {if(length(d)==1) d <- rep(d,nvar)} 114 if(is.null(a)) {a <- rep(1,nvar)} 115 theta <- rnorm(n,mu,sd) 116 item <- matrix(t(c+(z-c)/(1+exp(a*t((-theta %+% t( d)))))),n,nvar) 117 item[] <- rbinom(n*nvar, cat, item) 118 119 120 colnames(item) <- paste("V",1:nvar,sep="") 121 122 result <- list(items=item,discrimination=a,difficulty=d,gamma=c,zeta=z,theta=theta) 123 return (result) 124 } 125 126"sim.poly.ideal" <- 127function (nvar = 5 ,n = 500,low=-2,high=2,a=NULL,c=0,z=1,d=NULL, mu=0,sd=1,cat=5,mod="logistic") 128 { 129 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)} 130 return (result) 131 } 132 133 134"sim.poly.ideal.npl.absolute" <- 135function (nvar = 5 ,n = 500, low=-2,high=2,a=NULL,c=0,z=1,d=NULL,mu=0,sd=1,cat=5) 136 { cat <- cat -1 #binomial is based upon one fewer than categories 137 if(is.null(d)) {d <- seq(low,high,(high-low)/(nvar-1))} else {if(length(d)==1) d <- rep(d,nvar)} 138 if(is.null(a)) {a <- rep(1,nvar)} 139 theta <- rnorm(n,mu,sd) 140 item <- matrix(t(c+(z-c)/(1+exp(a*t(abs(-theta %+% t( d)))))),n,nvar) 141 item[] <- rbinom(n*nvar, cat, item) 142 143 144 colnames(item) <- paste("V",1:nvar,sep="") 145 146 result <- list(items=item,discrimination=a,difficulty=d,gamma=c,zeta=z,theta=theta) 147 return (result) 148 } 149 150"sim.poly.ideal.npl" <- 151function (nvar = 5 ,n = 500, low=-2,high=2,a=NULL,c=0,z=1,d=NULL,mu=0,sd=1,cat=5,theta=NULL) 152 { cat <- cat -1 #binomial is based upon one fewer than categories 153 if(is.null(d)) {d <- seq(low,high,(high-low)/(nvar-1))} else {if(length(d)==1) d <- rep(d,nvar)} 154 if(is.null(a)) {a <- rep(1,nvar)} 155 if(is.null(theta)) {theta <- rnorm(n,mu,sd)} 156 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) 157 p <- item 158 item[] <- rbinom(n*nvar, cat, item) 159 160 161 colnames(item) <- paste("V",1:nvar,sep="") 162 browser() 163 result <- list(p=p,items=item,discrimination=a,difficulty=d,gamma=c,zeta=z,theta=theta) 164 return (result) 165 } 166 167 168 169"sim.poly.ideal.npn" <- 170function (nvar = 5 ,n = 500, low=-2,high=2,a=NULL,c=0,z=1,d=NULL,mu=0,sd=1,cat=5) { 171warning("Not ready for prime time") 172 cat <- cat -1 #binomial is based upon one fewer than categories 173 if(is.null(d)) {d <- seq(low,high,(high-low)/(nvar-1))} else {if(length(d)==1) d <- rep(d,nvar)} 174 if(is.null(a)) {a <- rep(1,nvar)} 175 theta <- rnorm(n,mu,sd) # the latent variable 176 177 item <- matrix(t(c+(z-c)*pnorm(abs(a*t(theta %+% t(- d))))),n,nvar) #need to transpose and retranpose to get it right 178 #now convert these probabilities to outcomes 179 item[] <- rbinom(n*nvar, cat, item) 180 181 colnames(item) <- paste("V",1:nvar,sep="") 182 result <- list(items=item,discrimination=a,difficulty=d,gamma=c,zeta=z,theta=theta) 183 return (result) 184 } 185