## This is a single time point model. The basic idea here is that we ## can reuse the general framework to fit general IRT and CTT models. ## This model assumes that there is a common mean and standard ## deviation for all of the data points. ## This is loosely based on is.positive.definite from the package ## matrixcalc. However, the version I was using has a bug, and ## returns an error is situtations for which I want it to quietly ## return false. is.pdm <- function (m, tol = sqrt(.Machine$double.eps)) { if (!is.matrix(m) || !is.numeric(m) || any(is.na(m))) return(FALSE) if (ncol(m) != nrow(m)) return(FALSE) ## Test for symmetry if (any(abs(m-t(m)) > tol)) return (FALSE) lambda <- eigen(m,symmetric=TRUE,only.values=TRUE)$values return (all(lambda >= tol)) } TimelessNormalParam <- setClass("TimelessNormalParam", slots = c(mu="vector",Sigma="matrix", is.pd="logical", iSig="array"), ## prototype(mu=numeric(),Sigma=matrix(),is.pd=as.logical(NA), ## iSigma=matirx(NA)), contains="ParameterObject") "TimelessNormalParam" <- function (mu=numeric(),Sigma=matrix(),is.pd=is.pdm(Sigma), iSig=solve(Sigma)) { if (missing(Sigma) && !missing(iSig)) { is.pd <- is.dpm(iSig) if (is.pd) { Sigma <- solve(iSig) } else { Sigma <- matrix(NA) } } if (missing(iSig) && !is.pd) { iSig <- matrix(NA) } new("TimelessNormalParam",mu=mu, Sigma=Sigma, is.pd=is.pd, iSig=iSig) } setGeneric("getmu",function(param) standardGeneric("getmu")) setGeneric("getmu<-",function(param,value) standardGeneric("getmu<-")) setGeneric("getSigma",function(param) standardGeneric("getSigma")) setGeneric("getSigma<-",function(param,value) standardGeneric("getSigma<-")) setGeneric("is.pd",function(param) standardGeneric("is.pd")) setGeneric("getiSigma",function(param) standardGeneric("getiSigma")) setMethod("getmu","TimelessNormalParam", function(param) param@mu) setMethod("getmu<-","TimelessNormalParam", function(param,value) { param@mu <- value param }) setMethod("getSigma","TimelessNormalParam", function(param) param@Sigma) setMethod("getSigma<-","TimelessNormalParam", function(param,value) { param@Sigma <- value param@is.pd <- is.pdm(value) if (param@is.pd) { param@iSig <- solve(value) } else { param@iSig <- matrix(NA) } param }) setMethod("getiSigma","TimelessNormalParam", function(param) { if (any(is.na(param@iSig))) { solve(param@Sigma) } else { param@iSig } }) setMethod("is.pd","TimelessNormalParam", function (param) { if (is.na(param@is.pd)) { is.pdm(getSigma(param)) } else { param@is.pd } }) setMethod("pvec","TimelessNormalParam", function(param) c(param@mu,diag(param@Sigma),param@Sigma[lower.tri(param@Sigma)])) setMethod("pvec<-","TimelessNormalParam", function(param,value) { p <- length(param@mu) if (length(value) != p+p+p*(p-1)/2) stop("Parameter vector is incorrect length.") param@mu <- value[1:p] S <- diag(value[p+(1:p)]) S[lower.tri(S)] <- value[(2*p+1):length(value)] S <- t(S) S[lower.tri(S)] <- value[(2*p+1):length(value)] getSigma(param) <- S param }) ## Priors ## mu ~ N(muMean,Sigma/varWeight) ## Sigma ~ Wishart(precScale,Sdf) TimelessNormalPM <- setClass("TimelessNormalPM", slots =c(muMean="vector",varWeight="numeric", precScale="matrix",Sdf="numeric", preciScale="matrix"), ## prototype(param = "TimelessNormalParam", paramtype="TimelessNormalParam", ## muMean=numeric(), varWeight=NA_real_, ## precScale=matrix(), Sdf=NA_real_,preciScale=matrix()), contains="ProficiencyModel") TimelessNormalPM <- function (muMean=numeric(),varWeight=length(muMean), precScale=diag(length(muMean)), Sdf=length(muMean)+1, param=NullParameter()) { ## Check for consistent dimensions d <- length(muMean) if (length(varWeight) != 1L) { stop("VarWeight must be of length 1") } if (nrow(precScale) !=d || ncol(precScale) !=d) { stop("Dimensions of precision don't match dimensions of mean vector.") } if (!is.pdm(precScale)) { stop("Precision scale is not positive definite.") } if (!is(param,"TimelessNormalParam") && !is(param,"NullParameter")) { stop("Expected a TimelessNormal Parameter.") } if (Sdf <= d) { stop("Improper prior distribution for Sigma, df less than matrix rank.") } new("TimelessNormalPM",muMean=muMean,varWeight=varWeight, precScale=precScale, Sdf=Sdf, preciScale=solve(precScale), paramtype="TimelessNormalParam",param=param) } setMethod("nlatent","TimelessNormalPM", function(pm) { length(pm@muMean) }) setMethod("lpriorPMParam", "TimelessNormalPM", function(pm,param=parameters(pm),background=NA) { if (is.pd(param)) { dmvnorm(getmu(param),pm@muMean,getSigma(param)/pm@varWeight, log=TRUE) + log(dwish(getiSigma(param),pm@Sdf,pm@precScale)) } else { -Inf } }) setMethod("drawPMParam","TimelessNormalPM", function(pm,background=NA) { iSig <- riwish(pm@Sdf,pm@precScale) Sigma <- solve(iSig) mu <- as.vector(rmvnorm(1,pm@muMean,Sigma/pm@varWeight)) TimelessNormalParam(mu=mu, Sigma=Sigma, is.pd=TRUE, iSig=iSig) }) setMethod("drawInitialLatent", "TimelessNormalPM", function(pm,ncases=1,background=NA,param=parameters(pm)) { if (!is.pd(param)) { stop("Singular covariance matrix.") } rmvnorm(ncases,getmu(param),getSigma(param)) }) setMethod("lpriorLatent","TimelessNormalPM", function(pm,latent,background=NA,param=parameters(pm)) { if (is.pd(param)) { dmvnorm(latent,getmu(param),getSigma(param),log=TRUE) } else { -Inf } }) setMethod("optimalPMParams", "TimelessNormalPM", function (pm,latent,weights=1, background=NA,param=parameters(pm), control=list(),maxit=100,nowarn=FALSE, Bayes=TRUE) { result <- NULL if (length(weights)==1L) { weights <- rep(weights,nrow(latent)) } n <- sum(weights) xbar <- apply(sweep(latent,1,weights,"*"),2,sum)/n res <- sweep(sweep(latent,2,xbar),1,weights,"*") S <- t(res)%*%(res/weights) if (Bayes) { dfnew <- pm@Sdf + n wnew <- pm@varWeight + n munew <- (pm@varWeight*pm@muMean + n*xbar)/wnew lambda <- pm@preciScale+S+ (pm@varWeight*n)/wnew*outer(xbar-pm@muMean,xbar-pm@muMean,"*") result <- TimelessNormalParam(mu=munew, Sigma=lambda/dfnew) } else { result <- TimelessNormalParam(mu=xbar, Sigma=S/(n-1)) } dev <- -2*sum(weights*lpriorLatent(pm,latent,param=result)) if (Bayes) { dev <- dev -2*lpriorPMParam(pm,result) } list(param=result,deviance=dev,converged=TRUE) })