# Simulate Godfrey Thomson's "sampling model" of mental abilities, and perform # factor analysis on the resulting test scores. # MASS package used to make random multivariate normal vectors require(MASS) # stats package used to do factor analysis, but if it's missing you've got big # problems! require(stats) # Set the mean and standard deviation of each elementary ability to that of an # IQ test - not actually important! ability.mean = 100 ability.sd = 15 thomson.model <- function(n,d,a,q,output="variance") { # Using incomprehensible parameter names is bad # number of testees = n # number of tests = d # number of shared abilities = a # max. number of specific abilities per test = q # assign abilities to tests general.per.test = floor(runif(d,min=0,max=a)) + 1 specifics.per.test = floor(runif(d,min=0,max=q)) + 1 # Define the matrix assigning abilities to tests general.to.tests = matrix(0,a,d) # The use of a for loop here is maybe a little slower than some sort # of vectorization, but it's sanity-preserving; so is using the temporary # variable "abilties" to hold the sample. for (i in 1:d) { abilities = sample(1:a,size=general.per.test[i],replace=FALSE) general.to.tests[abilities,i] = 1 } # Covariance matrix of the general abilities sigma = matrix(0,a,a) diag(sigma) = (ability.sd)^2 mu=rep(ability.mean,a) x = mvrnorm(n,mu,sigma) # person-by-abilities matrix of abilities # The "general" part of the tests general.tests = x %*% general.to.tests specific.tests = matrix(0,n,d) noisy.tests = matrix(0,n,d) # Each test gets its own specific abilities, which are independent for each # person # Again, I could probably vectorize, but d is small and this is saner for (i in 1:d) { # Each test has noises.per.test disturbances, each of which has s.d. 15, # since these are all independent their variances add j = specifics.per.test[i] specifics = rnorm(n,mean=100*j,sd=ability.sd*sqrt(j)) specific.tests[,i] = general.tests[,i] + specifics # Finally, for extra realism, some mean-zero trial-to-trial noise, so # that if we re-use this combination of general and specific ability # scores, we won't get the exact same test scores twice noises = rnorm(n,mean=0,sd=ability.sd) noisy.tests[,i] = specific.tests[,i] + noises } # do the factor analysis my.fa = factanal(noisy.tests,1) # extract the loadings v = as.vector(my.fa\$loadings) # compute how much of the variance (really, correlation) the factor # describes varexp = sum(v^2)/d # output: should do a return() instead to be completely pukka switch(output, variance = varexp, everything = list(general.ability.pattern = general.to.tests, numbers.of.specifics = specifics.per.test, ability.matrix = x, specific.tests = specific.tests, data = noisy.tests, loadings = v, varexp = varexp), "I can't tell what you want me to output") }