# # Lord's Unidimensional Normal-Ogive Latent Trait Model # (also estimated with 1 and 2 parameter IRT models and # compared to the FA on the tetrachoric correlations; # conversion formulaes for the parameters are given). # # Set seed so we can repeat same simulation each time # simulation is run set.seed(123456) # Load libraries library(polycor) library(ltm) # Number of subjects n<-1000 # Create ability scores for n individuals # Assuming N(0,1) f1<-rnorm(n) # Create n response thresholds (t) on each of 6 items (t1-t6) # for n individuals assuming N(0,1) for t1-t6. # # Each item has a mean threshold across individuals of mean=0, sd=1 # Thresholds for items generated for each individual on each item t1<-rnorm(n,0,1) t2<-rnorm(n,0,1) t3<-rnorm(n,0,1) t4<-rnorm(n,0,1) t5<-rnorm(n,0,1) t6<-rnorm(n,0,1) # Loadings or Discriminations (slopes) b1<-.6 b2<-.6 b3<-.6 b4<-.6 b5<-.6 b6<-.6 # Item errors (based on identity: uniqueness = 1- lambda^2) e1<-rnorm(n,0,sqrt(1-b1^2)) e2<-rnorm(n,0,sqrt(1-b2^2)) e3<-rnorm(n,0,sqrt(1-b3^2)) e4<-rnorm(n,0,sqrt(1-b4^2)) e5<-rnorm(n,0,sqrt(1-b5^2)) e6<-rnorm(n,0,sqrt(1-b6^2)) # Linear function to generate the continuous "unobserved" trait # based on the item function (conditioned on, or given f); f are # the latent trait scores; b is discriminiation or loading (slope) # for an item; and e is item error; Note that the difficulites are # set to 0. ability<-function(b,f1,e){ Z<-b*f1+e return(Z) } # We need to threshold the unobserved trait scores based on item function (Z), # so that scores become an observed score value of 1 or 0 (y1, y2...etc). # Compare Z (trait) and t (threshold) for each item); # If f<t generate observed score of 0; # If f>t generate observed score of 1. y1<-ifelse(ability(b=b1, e=e1, f1) < t1 , 0 ,1) y2<-ifelse(ability(b=b2, e=e2, f1) < t2 , 0 ,1) y3<-ifelse(ability(b=b3, e=e3, f1) < t3 , 0 ,1) y4<-ifelse(ability(b=b4, e=e4, f1) < t4 , 0 ,1) y5<-ifelse(ability(b=b5, e=e5, f1) < t5 , 0 ,1) y6<-ifelse(ability(b=b6, e=e6, f1) < t6 , 0 ,1) # Gather data into a single data frame f1.data<-data.frame(cbind(y1, y2, y3, y4, y5, y6)) # Display first six and last six records head(f1.data) tail(f1.data) # Rasch model - discrimination assumed to be equal f1.data.rasch.fit<-rasch(f1.data) # Two parameter IRT - discrimination allowed to vary f1.data.ltm.fit<-ltm(f1.data~z1) # Summarize model and plot item characteristic curves # # Note: Coefficients in the IRT parameterization, if divided # by D=1.701, give coefficient values that can be used in # the formula below, to calculate the "factor loading" # values that would be obtained by performing a linear 1-factor # common trait model on the "tetrachoric" correlation matrix. # # For example: estimated IRT b = .7288 (slope or discrim. from rasch fit) # new_b = .7288/1.701 = .4285 # # A comparable estimate of lambda = .4284/sqrt(1+.4285^2) # = .393 # # A factor analysis on the tetrachoric correlation matrix gives: # # lambda = .391 # # .393 is fairly close to .391 # # summary(f1.data.rasch.fit) plot(f1.data.rasch.fit) GoF.rasch(f1.data.rasch.fit) summary(f1.data.ltm.fit) plot(f1.data.ltm.fit) # Estimate and print out ability scores for n respondents f1.data.ltm.scores<-factor.scores(f1.data.ltm.fit) f1.data.rasch.scores<-factor.scores(f1.data.rasch.fit) head(f1.data.ltm.scores) head(f1.data.rasch.scores) # Calculate polychoric correlation matrix # Create ordinal data from numeric data # Make 1 and 0's categorical f1.data$y1<-as.factor(f1.data$y1) f1.data$y2<-as.factor(f1.data$y2) f1.data$y3<-as.factor(f1.data$y3) f1.data$y4<-as.factor(f1.data$y4) f1.data$y5<-as.factor(f1.data$y5) f1.data$y6<-as.factor(f1.data$y6) library(polycor) f1.poly.cor<-hetcor(f1.data,bins=10) f1.poly.cor f1.poly.cor<-f1.poly.cor$correlations f1.poly.cor # FA Loadings based on Polychoric correlation # # Note that these loadings can be used to calculate the # discrimination parameters from an IRT parameterization # whose coefficients have been divided by D=1.701 (see earlier) # discussion. # # For example: # # estimated lamda loading = .391 # # comparable IRT loading = .391/sqrt(1-.391^2) # = .4248 # # comparable estimate of IRT loading = .4248 * 1.701 = .7226 # # A Two parameter IRT analysis gives .6901 # # .7266 is close to .6901 # # Still very close.... f1.poly.fact<-factanal(covmat=f1.poly.cor, factors=1) f1.poly.fact ########################################################### # # Omega Reliability calculated from FA loadings based on # polychoric correlations. # # Calculate Omega Coefficient PartA<-(sum(f1.poly.fact$loadings))^2 PartB<-sum(1-f1.poly.fact$loadings^2) omega.poly<-PartA/(PartA+PartB) omega.poly