version 0.2

cran · Jul 9, 2012 · f47450e · f47450e
commit f47450e
Show file tree

Hide file tree

Showing 38 changed files with 8,043 additions and 0 deletions.
diff --git a/DESCRIPTION b/DESCRIPTION
@@ -0,0 +1,17 @@
+Package: pln
+Type: Package
+Title: Polytomous logit-normit (graded logistic) model estimation
+Version: 0.2
+Date: 2012-07-09
+Author: Carl F. Falk and Harry Joe
+Maintainer: Carl F. Falk <cffalk@gmail.com>
+Description: Performs bivariate composite likelihood and full
+        information maximum likelihood estimation for polytomous
+        logit-normit (graded logistic) item response theory (IRT)
+        models.
+Depends: R (>= 2.9.0)
+License: GPL-3
+LazyLoad: yes
+Packaged: 2012-07-10 21:34:10 UTC; cf
+Repository: CRAN
+Date/Publication: 2012-07-19 05:27:08
diff --git a/MD5 b/MD5
@@ -0,0 +1,37 @@
+dc9f45c6ff325941703c0f88f503cefd *DESCRIPTION
+96fa608e6e90c4841d918a82f251a679 *NAMESPACE
+2d975a3685008f1ac01c892441c15638 *R/m2pln.R
+37ef7cab2eab33af2d0463eb1b39aecc *R/m2rasch.R
+cdd695d35730325c7f28882f810f41c8 *R/nrbcpln.R
+33cb9672f058b204fe5b69ad99dee2f1 *R/nrmlepln.R
+ace929c3a7765da94de074c175fa76e4 *R/nrmlerasch.R
+00eb6eb752264562e9a791c06a8fb029 *R/startvals.R
+78d2795209981141b205250790f3e0af *R/utils.R
+3568c0662d33c4d872ba8c966eed100b *TODO
+64aa82663accb1e792cbaacbbcf63481 *data/item5fr.tab.gz
+2ddbafc19058e8dc4357259223a5601a *data/item9cat5.tab.gz
+4fc32462d0ab9e04508b457407d1edd5 *man/item5fr.Rd
+4dc73a12e4fc7bc6cde546dc02223f38 *man/item9cat5.Rd
+07482eead9613a69f5cee84d42ffecc7 *man/nrmlepln.Rd
+5e05a853686bc3da796a788f61f6016c *man/pln-package.Rd
+522ab09f7a4d7d5b7c7e81975a515dbe *man/startalphas.Rd
+8d55714705b3deaf64460c8e0ff83bbc *src/Makevars
+8d55714705b3deaf64460c8e0ff83bbc *src/Makevars.win
+4bb93d284e3331fe960746e32cb14113 *src/amatrix.c
+e54c39ac8ca747ecf6c85ec613c281b0 *src/bcplncov.c
+85dcc527be0e0f723254ca42182ef018 *src/cmata.c
+4e463971b445bd2a23437f9cfc79729a *src/d2v.c
+7c6dcc562853e7b820d5f0c78f1515e4 *src/geppldet.c
+fa22309f82c41a6c963e35ec321f5c23 *src/m2.c
+f9b36ed04f9d897008925092a711d593 *src/m2rasch.c
+d8b0735f16f82e5f5a6dba42c98fd535 *src/nrbcpln.c
+4df5fe880b205029ac724007988ab409 *src/nrmin.c
+5303dbd7798443133612c581886f2cd2 *src/nrminbcl.c
+31c20de51c661d80059524d3f687aa87 *src/nrmlepln.c
+a3b5b67ae55af26beac105579bbc34ed *src/nrmlerasch.c
+a868e69935aebdc03fb6574560fab511 *src/polyder2.c
+fe28a9b49338890a45f692d0d031a33f *src/simulpln.c
+f1965d423e20f916fc1399b998d284db *src/ssgauher.c
+e6bb3740d031f440c779e1caf4bbe488 *src/startpln.c
+f4d009060a4053bf7d034c1c6564cddf *src/x2statb.c
+772c1068f26ddd8395d16bcc5516b499 *tests/Examples/pln-Ex.Rout.save
diff --git a/NAMESPACE b/NAMESPACE
@@ -0,0 +1,6 @@
+useDynLib(pln)
+export(nrbcpln)
+export(nrmlepln)
+export(nrmlerasch)
+export(startalphas)
+export(startbetas)
diff --git a/R/m2pln.R b/R/m2pln.R
@@ -0,0 +1,29 @@
+m2pln.func <- function(nitem, ncat, nrec, myX, alphas, betas, nq, iprint){
+
+  ## prep return variables
+  g<-nitem*(ncat-1)
+  for (i in 1:(nitem-1))
+  {
+    for (j in 0:(i-1))
+    {
+      for (k in 1:(ncat-1))
+      {
+        for (l in 1:(ncat-1))
+        {
+          g<-g+1
+        }
+      }
+    }
+  }
+
+  samplemout=rep(0,g)
+  m2statout=0
+  dfout=0
+
+  ## TO DO: add package arument here?
+  out <- .C("Rm2",
+    as.integer(nitem), as.integer(ncat), as.integer(nrec), as.double(myX),
+    as.double(alphas), as.double(betas), samplemout=as.double(samplemout),
+    m2statout=as.double(m2statout), dfout=as.double(dfout), as.integer(nq))  
+  list(samplem=out$samplemout, m2stat=out$m2statout, df=out$dfout)  
+}
diff --git a/R/m2rasch.R b/R/m2rasch.R
@@ -0,0 +1,29 @@
+m2rasch.func <- function(nitem, ncat, nrec, myX, alphas, betas, nq, iprint){
+
+  ## prep return variables
+  g<-nitem*(ncat-1)
+  for (i in 1:(nitem-1))
+  {
+    for (j in 0:(i-1))
+    {
+      for (k in 1:(ncat-1))
+      {
+        for (l in 1:(ncat-1))
+        {
+          g<-g+1
+        }
+      }
+    }
+  }
+
+  samplemout=rep(0,g)
+  m2statout=0
+  dfout=0
+
+  ## TO DO: add package arument here?
+  out <- .C("Rm2rasch",
+    as.integer(nitem), as.integer(ncat), as.integer(nrec), as.double(myX),
+    as.double(alphas), as.double(betas), samplemout=as.double(samplemout),
+    m2statout=as.double(m2statout), dfout=as.double(dfout), as.integer(nq))  
+  list(samplem=out$samplemout, m2stat=out$m2statout, df=out$dfout)  
+}
diff --git a/R/nrbcpln.R b/R/nrbcpln.R
@@ -0,0 +1,67 @@
+nrbcpln<- function (x, ncat, nitem=NULL, alphas=NULL, betas=NULL, abound=c(-10,10),
+    bbound=c(-1,10), nq=48, mxiter=200, se=TRUE, iprint=FALSE) {
+
+  myInput<-check.input(x, ncat, nitem, nq, mxiter, iprint)
+
+  ## get starting values if not present already
+  if(!check.alphas(alphas, myInput$nitem, myInput$ncat)){
+    alphas<-startpln.func(myInput$nitem, myInput$ncat, myInput$nrec, myInput$myX)$alphas
+  }
+
+  ## prep betas
+  if(!check.betas(betas, myInput$nitem)){
+    betas<-startbetas.func(myInput$myX)
+  }
+
+  ## check bounds
+  abound<-check.bounds(alphas, abound)
+  bbound<-check.bounds(betas, bbound)
+
+  nrbcplnout <- nrbcpln.func(myInput$ncat, myInput$nitem, myInput$nrec, myInput$myX, alphas,
+        betas, abound, bbound, myInput$nq, myInput$mxiter, myInput$iprint)
+
+  alphas<-nrbcplnout$bcpln[1:((myInput$ncat-1)*myInput$nitem)]
+  betas<-nrbcplnout$bcpln[((myInput$ncat-1)*myInput$nitem+1):(myInput$ncat*myInput$nitem)]
+
+  out<-list(alphas=alphas,betas=betas,nllk=nrbcplnout$nrbcpln,conv=nrbcplnout$iconv)
+
+  if(se){
+    nrbcplncov<-nrbcplncov.func(myInput$nitem, myInput$ncat, myInput$nrec, alphas, betas,
+      myInput$N, myInput$nq, myInput$iprint)
+
+    V<-matrix(nrbcplncov$V, nrow=myInput$nitem*myInput$ncat,ncol=myInput$nitem*myInput$ncat)
+    seVec<-sqrt(diag(V))
+
+    sealphas<-seVec[1:((myInput$ncat-1)*myInput$nitem)]
+    sebetas<-seVec[((myInput$ncat-1)*myInput$nitem+1):(myInput$ncat*myInput$nitem)]
+
+    out<-append(out,list(sealphas=sealphas, sebetas=sebetas, vcov=V)) 
+  }
+
+  return(out)
+
+}
+
+nrbcpln.func<-function(ncat, nitem, nrec, myX, alphas, betas, abound, bbound, nq,
+    mxiter, iprint){
+  nrbcplnout<-0
+  iconv<-0
+  np<-ncat*nitem
+  bcplnout<-rep(0,np)
+  out <- .C("Rnrbcpln",
+    as.integer(nitem), as.integer(ncat), as.integer(nrec), as.double(myX),
+    as.double(alphas), as.double(betas), as.double(abound), as.double(bbound),
+    nbcplnout=as.double(nrbcplnout), bcplnout=as.double(bcplnout), as.integer(nq),
+    as.integer(mxiter), iconv=as.integer(iconv), as.integer(iprint))
+  list(nrbcpln=out$nbcplnout,bcpln=out$bcplnout, iconv=out$iconv)
+}
+
+## asymptotic covariance matrix
+nrbcplncov.func<-function(nitem, ncat, nrec, alphas, betas, N, nq, iprint) {
+  V<-rep(0,nitem*ncat*nitem*ncat)
+  params<-c(alphas,betas)
+  out<-.C("Rbclcov",
+          as.integer(nitem), as.integer(ncat), as.integer(N),
+          as.double(params), V=as.double(V), as.integer(nq), as.integer(iprint))
+  return(out)
+}
diff --git a/R/nrmlepln.R b/R/nrmlepln.R
@@ -0,0 +1,67 @@
+nrmlepln <- function (x, ncat, nitem=NULL, alphas=NULL, betas=NULL, abound=c(-10,10),
+    bbound=c(-1,10), nq=48, mxiter=200, m2=TRUE, iprint=FALSE) {
+
+  myInput<-check.input(x, ncat, nitem, nq, mxiter, iprint)
+
+  ## get starting values if not present already
+  if(!check.alphas(alphas, myInput$nitem, myInput$ncat)){
+    alphas<-startpln.func(myInput$nitem, myInput$ncat, myInput$nrec, myInput$myX)$alphas
+  }
+
+  ## prep betas
+  if(!check.betas(betas, myInput$nitem)){
+    betas<-startbetas.func(myInput$myX)
+  }
+
+  ## check bounds
+  abound<-check.bounds(alphas, abound)
+  bbound<-check.bounds(betas, bbound)
+
+  nrmleplnout <- nrmlepln.func(myInput$nitem, myInput$ncat, myInput$nrec, myInput$myX, alphas,
+        betas, abound, bbound, myInput$nq, myInput$mxiter, myInput$iprint)
+  alphas<-nrmleplnout$mlePlnOut[1:((myInput$ncat-1)*myInput$nitem)]
+  betas<-nrmleplnout$mlePlnOut[((myInput$ncat-1)*myInput$nitem+1):(myInput$ncat*myInput$nitem)]
+
+  V<-matrix(nrmleplnout$invHesOut, nrow=myInput$nitem*myInput$ncat,ncol=myInput$nitem*myInput$ncat)
+  seVec<-nrmleplnout$seVecOut
+
+  sealphas<-seVec[1:((myInput$ncat-1)*myInput$nitem)]
+  sebetas<-seVec[((myInput$ncat-1)*myInput$nitem+1):(myInput$ncat*myInput$nitem)]
+
+
+  out<-list(alphas=alphas,betas=betas,nllk=nrmleplnout$nllkOut,sealphas=sealphas,
+    sebetas=sebetas,invhes=V)
+
+  if(m2){
+    m2out<-m2pln.func(myInput$nitem, myInput$ncat, myInput$nrec, myInput$myX, alphas, betas,
+      myInput$nq, myInput$iprint)
+    pval<-pchisq(m2out$m2stat, m2out$df, lower.tail=FALSE)
+      ##out<-append(out, list(samplemout=m2out$samplem, m2stat=m2out$m2stat, df=m2out$df, pval=pval))      
+    out<-append(out, list(teststat=m2out$m2stat, df=m2out$df, pval=pval))
+  }
+
+  return(out)
+}
+
+nrmlepln.func <- function(nitem, ncat, nrec, myX, alphas, betas, abound, bbound, nq, mxiter,
+    iprint){
+
+  ## prep return variables
+  nllkOut<-0
+  iconv<-0
+  np<-ncat*nitem
+  mlePlnOut<-rep(0,np)
+  seVecOut<-rep(0,np)
+  invHesOut<-matrix(0,nrow=np,ncol=np)
+
+  ## TO DO: add PACKAGE argument here
+  out <- .C("Rnrmlepln",
+            as.integer(nitem), as.integer(ncat), as.integer(nrec), as.double(myX),
+            as.double(alphas), as.double(betas), as.double(abound), as.double(bbound),
+            nllkOut=as.double(nllkOut), mlePlnOut=as.double(mlePlnOut),
+            seVecOut=as.double(seVecOut), invHesOut=as.double(invHesOut),
+            as.integer(nq), as.integer(mxiter), iconv=as.integer(iconv),
+            as.integer(iprint))
+  return(out)
+
+}
diff --git a/R/nrmlerasch.R b/R/nrmlerasch.R
@@ -0,0 +1,72 @@
+nrmlerasch <- function (x, ncat, nitem=NULL, alphas=NULL, abound=c(-10,10),
+    bbound=c(-1,10), nq=48, mxiter=200,  m2=TRUE, iprint=FALSE) {
+
+  myInput<-check.input(x, ncat, nitem, nq, mxiter, iprint)
+
+  ## get starting values if not present already
+  if(!check.alphas(alphas, myInput$nitem, myInput$ncat)){
+    alphas<-startpln.func(myInput$nitem, myInput$ncat, myInput$nrec, myInput$myX)$alphas
+  }
+
+  ## prep betas
+#  if(is.null(betas) || betas%%1!=0 || length(betas)>1){
+    betas<-1
+#  }
+
+  ## check bounds
+  abound<-check.bounds(alphas, abound)
+  bbound<-check.bounds(betas, bbound)      
+
+  nrmleraschout <- nrmlerasch.func(myInput$nitem, myInput$ncat,
+        myInput$nrec, myInput$myX, alphas, betas, abound, bbound,
+        myInput$nq, myInput$mxiter, myInput$iprint)
+
+  alphas<-nrmleraschout$mlePlnOut[1:((myInput$ncat-1)*myInput$nitem)]
+  betas<-nrmleraschout$mlePlnOut[((myInput$ncat-1)*myInput$nitem)+1]
+
+  V<-matrix(nrmleraschout$invHesOut, nrow=((myInput$ncat-1)*myInput$nitem)+1,
+        ncol=((myInput$ncat-1)*myInput$nitem)+1)
+  seVec<-nrmleraschout$seVecOut
+
+  sealphas<-seVec[1:((myInput$ncat-1)*myInput$nitem)]
+  sebetas<-seVec[((myInput$ncat-1)*myInput$nitem)+1]
+
+
+  out<-list(alphas=alphas,betas=betas,nllk=nrmleraschout$nllkOut,sealphas=sealphas,
+    sebetas=sebetas,invhes=V)
+
+  if(m2){
+    m2out<-m2rasch.func(myInput$nitem, myInput$ncat, myInput$nrec, myInput$myX,
+      alphas, betas, myInput$nq, myInput$iprint)
+    pval<-pchisq(m2out$m2stat, m2out$df, lower.tail=FALSE)
+      ##out<-append(out, list(samplemout=m2out$samplem, m2stat=m2out$m2stat, df=m2out$df,pval=pval))
+    out<-append(out, list(teststat=m2out$m2stat, df=m2out$df,pval=pval))
+  }
+
+  return(out)
+
+}
+
+nrmlerasch.func <- function(nitem, ncat, nrec, myX, alphas, betas, abound, bbound,
+    nq, mxiter,iprint){
+
+  ## prep return variables
+  nllkOut<-0
+  iconv<-0
+  np<-ncat*nitem-nitem+1
+  mlePlnOut<-rep(0,np)
+  seVecOut<-rep(0,np)
+  invHesOut<-matrix(0,nrow=np,ncol=np)
+
+  ## TO DO: add PACKAGE argument here
+  ## TO DO: coerce invHes into a matrix?
+  out <- .C("Rnrmlerasch",
+            as.integer(nitem), as.integer(ncat), as.integer(nrec), as.double(myX),
+            as.double(alphas), as.double(betas), as.double(abound),
+            as.double(bbound),nllkOut=as.double(nllkOut),
+            mlePlnOut=as.double(mlePlnOut), seVecOut=as.double(seVecOut),
+            invHesOut=as.double(invHesOut), as.integer(nq), as.integer(mxiter),
+            iconv=as.integer(iconv), as.integer(iprint))
+  return(out)
+
+}
diff --git a/R/startvals.R b/R/startvals.R
@@ -0,0 +1,54 @@
+startalphas<- function (x, ncat, nitem=NULL) {
+
+  ## input checking & data prep
+  myInput<-check.input(x, ncat, nitem)
+
+  out <- startpln.func(myInput$nitem, myInput$ncat, myInput$nrec, myInput$myX)
+
+  return(out$alphas)
+
+}
+
+startpln.func<-function(nitem, ncat, nrec, myX){
+
+  ## prep return variables
+  alphas=rep(0,nitem*(ncat-1))
+
+  ## TO DO: add PACKAGE argument here
+  out <- .C("Rstartpln",
+            as.integer(nitem), as.integer(ncat), as.integer(nrec), as.double(myX),
+            alphas=as.double(alphas))
+  return(out)
+}
+
+startbetas<-function(x, ncat, nitem=NULL){
+  myInput<-check.input(x, ncat, nitem)
+  x<-myInput$myX
+  betas<-startbetas.func(x)
+  return(betas)
+}
+
+startbetas.func<-function(x){
+  nn<-ncol(x)
+  nitem<-nn-1
+  y<-x[,-nn]
+  fr<-x[,nn]
+  tot<-sum(fr)
+  mnvec<-apply(y*fr,2,sum)/tot ## means
+  vvec<-apply(y*y*fr,2,sum)/tot ## variances
+  vvec<-vvec-mnvec^2
+  cc<-matrix(0,nitem,nitem)
+  for(j in 1:(nitem-1))
+  { for(k in (j+1):nitem)
+    { ss<-sum(y[,j]*y[,k]*fr)/tot
+      den<-vvec[j]*vvec[k]
+      cc[j,k]<-(ss-mnvec[j]*mnvec[k])/sqrt(den)
+   	  cc[k,j]<-cc[j,k]
+    }
+  }
+  avcc<-apply(cc,2,sum)/(nitem-1)
+  #print(avcc)
+  bvec<-avcc*10; bvec[bvec>=1]=1; bvec[bvec<= -.2]= -.2
+  #print(bvec)
+  bvec
+}