mw {ldsa}R Documentation

~~function to do ... ~~

Description

~~ A concise (1-5 lines) description of what the function does. ~~

Usage

mw(y, n = NULL, lat = NULL, freqcount = FALSE, niter = 100, err.mode = c("symhom", "asymhom", "symitem", "asymitem"), recover.df = TRUE, item.labels = colnames(y)[1:(NCOL(y) - freqcount)], tol = 1e-10, verbose = FALSE)

Arguments

y ~~Describe y here~~
n ~~Describe n here~~
lat ~~Describe lat here~~
freqcount ~~Describe freqcount here~~
niter ~~Describe niter here~~
err.mode ~~Describe err.mode here~~
recover.df ~~Describe recover.df here~~
item.labels ~~Describe item.labels here~~
tol ~~Describe tol here~~
verbose ~~Describe verbose here~~

Details

~~ If necessary, more details than the description above ~~

Value

~Describe the value returned If it is a LIST, use

comp1 Description of 'comp1'
comp2 Description of 'comp2'

...

Warning

....

Note

~~further notes~~

~Make other sections like Warning with section{Warning }{....} ~

Author(s)

~~who you are~~

References

~put references to the literature/web site here ~

See Also

~~objects to See Also as help, ~~~

Examples

##---- Should be DIRECTLY executable !! ----
##-- ==>  Define data, use random,
##--    or do  help(data=index)  for the standard data sets.

## The function is currently defined as
function(y,n=NULL,lat=NULL,freqcount=FALSE,niter=100,err.mode=c("symhom","asymhom","symitem","asymitem"),recover.df=TRUE,item.labels=colnames(y)[1:(NCOL(y)-freqcount)],tol=1e-10,verbose=FALSE){
  err.mode<-match.arg(err.mode)  #Clarify this parameter
  #Preprocess y, if need be
  if(!freqcount)
    y<-mw.enumerate(y)
  if(is.null(item.labels))
    item.labels<-1:(NCOL(y)-1)
  colnames(y)<-c(item.labels,"Count")
  #Make sure that we've got the input lattice
  if(is.null(lat)){  #If the lattice isn't given, extract it from the data
    if(is.null(n))
      stop("One of n or lat must be specified.\n")
    else
      lat<-mw.intclose(mw.freqsel(y,n,freqcount=TRUE))
  }
  #Get some useful stats
  nmb<-NCOL(y)-1      #Number of micro-beliefs
  nl<-NROW(lat)       #Number of lattice elements
  npat<-2^nmb         #Number of potential patterns
  n<-sum(y[,nmb+1])   #Total number of observations
  b1<-y[,1:nmb]       #Patterns - 1-coded
  b0<-1-b1            #Patterns - 0-coded
  #Create initial class probabilities (randomly)
  cp<-runif(nl)
  cp<-cp/sum(cp)
  names(cp)<-1:nl
  #Create initial error rates (again, randomized)
  err<-switch(err.mode,
    symhom=matrix(runif(1,0,0.5),nr=nmb,nc=2),
    asymhom=sapply(runif(2,0,0.5),rep,nmb),
    symitem=t(sapply(runif(nmb,0,0.5),rep,2)),
    asymitem=matrix(runif(nmb*2,0,0.5),nmb)
  )
  #Do some time-saving preprocessor work
  pas<-matrix(0,nl,npat)
  cs<-nl
  errfp<-array(0,dim=c(nmb,nl,npat))    #False positives and false negatives
  errfn<-array(0,dim=c(nmb,nl,npat))    #For each item, by pattern and class
  for(i in 1:nl){
    errfn[,i,]<-t(sweep(b0,2,lat[i,],"&"))
    errfp[,i,]<-t(sweep(b1,2,1-lat[i,],"&"))
  }  
  #Perform an EM loop, hopefully converging to the Truth(TM)
  for(i in 1:niter){
    #Calculate the pattern probability distribution
    p0<-sweep(lat,2,err[,1],"*")+sweep(1-lat,2,1-err[,2],"*")
    p1<-1-p0
    #Find the conditional probabilities for each lattice member (i.e. class)
    for(j in 1:nl){
      #Compute p(pattern|class)
      bprob<-sweep(b1,2,p1[j,],"*")+sweep(b0,2,p0[j,],"*") #Calculate response probs
      pas[j,]<-drop(exp(log(bprob)%*%rep(1,nmb)))  #Aggregate response probs
    }
    #Compute p(pattern|class)*p(class)
    pas<-sweep(pas,1,cp,"*")
    #Compute p(class|pattern)
    pp<-drop(rep(1,nl)%*%pas)
    pas<-t(t(pas)/pp)    
    #Compute the expected number of cases per pattern per class
    spas<-t(t(pas)*y[,nmb+1])
    cs<-drop(spas%*%rep(1,npat))
    cp<-cs/n   #Normalize by total number of cases to get class probs
    #Calculate expected error rates
    if(err.mode=="symhom"){
      temp<-sum(sweep(errfn,c(2,3),spas,"*"))+sum(sweep(errfp,c(2,3),spas,"*"))
      temp<-max(1e-20,min(0.5-1e-20,temp/(nmb*n)))
      err<-matrix(temp,nr=nmb,nc=2)
    }else if(err.mode=="asymhom"){
      temp<-c(sum(sweep(errfn,c(2,3),spas,"*")),sum(sweep(errfp,c(2,3),spas,"*")))
      #temp<-pmax(temp/(c(mean(lat),mean(1-lat))*n),1e-20)
      temp<-pmax(temp/c(sum(apply(lat,1,sum)*cp)*n,sum(apply(1-lat,1,sum)*cp)*n),1e-20)
      temp<-temp/max(1,sum(temp))
      err<-sapply(temp,rep,nmb)
    }else if(err.mode=="symitem"){
      temp<-apply(sweep(errfn,c(2,3),spas,"*"),1,sum)+apply(sweep(errfp,c(2,3),spas,"*"),1,sum)
      err[,1]<-pmax(pmin(temp/n,0.5-1e-20),1e-20)
      err[,2]<-err[,1]
    }else if(err.mode=="asymitem"){
      temp<-apply(sweep(errfn,c(2,3),spas,"*"),1,sum)
      err[,1]<-temp/(cp%*%lat*n)
      temp<-apply(sweep(errfp,c(2,3),spas,"*"),1,sum)
      err[,2]<-temp/(cp%*%(1-lat)*n)
      #Patch the errors, if need be
      err<-pmax(pmin(err,1-1e-20),1e-20)
      err<-sweep(err,1,pmax(apply(err,1,sum),1),"/")
    }
    #If desired, dump some diagnostics at this point
    if(verbose){
      cat("Iteration",i,"\n")
    }
  }
  #Calculate the final pattern probability distribution
  p0<-sweep(lat,2,err[,1],"*")+sweep(1-lat,2,1-err[,2],"*")
  p1<-1-p0
  #Find the conditional probabilities for each lattice member (i.e. class)
  for(j in 1:nl){
    #Compute final p(pattern|class)p(class)
    bprob<-sweep(b1,2,p1[j,],"*")+sweep(b0,2,p0[j,],"*") #Calculate response probs
    pas[j,]<-drop(exp(log(bprob)%*%rep(1,nmb)))*cp[j]  #Aggregate response probs
  }
  #If recovering degrees of freedom, do so now
  if(recover.df){
    df.recovered<-sum(cp<tol)+sum(err<tol)
  }else
    df.recovered<-0
  #Calculate the maximized likelihood
  ll<-sum(log(drop(rep(1,nl)%*%pas))*y[,nmb+1])
  yc<-y[y[,nmb+1]>0,nmb+1]
  ll0<-sum(yc*log(yc/n))
  #Fit the null model
  prob.null<-as.vector(t(y[,1:nmb])%*%y[,nmb+1])/n
  names(prob.null)<-item.labels
  
  lln<-sum((log(sweep(y[,1:nmb],2,prob.null,"*")+sweep(1-y[,1:nmb],2,1-prob.null,"*"))%*%rep(1,nmb))*y[,nmb+1])
  #Calculate the entropies of the saturated, null, and fitted distributions
  temp<-log2(drop(rep(1,nl)%*%pas))    #Base 2 logs for fitted
  etot<-sum(-temp*(2^temp),na.rm=T) 
  edat<-sum(-temp*(2^temp)*y[,nmb+1],na.rm=T)
  temp<-log2(y[,nmb+1]/n)  #Saturated model
  etot0<-sum(-temp*(2^temp),na.rm=T) 
  edat0<-sum(-temp*(2^temp)*y[,nmb+1],na.rm=T)
  temp<-log2(sweep(y[,1:nmb],2,prob.null,"*")+sweep(1-y[,1:nmb],2,1-prob.null,"*"))%*%rep(1,nmb)  #Null model
  etotn<-sum(-temp*(2^temp),na.rm=T) 
  edatn<-sum(-temp*(2^temp)*y[,nmb+1],na.rm=T)
  #Dump the final state, if desired
  if(verbose){
    cat("Log-likelihood:",ll,"Saturated Log-Likelihood:",ll0,"\n")
    cat("Likelihood Ratio Chi-square:",2*(ll0-ll),"BIC:",-2*ll+log(n)*(nl*(nmb+1)-1),"\n")
    cat("Class prob:",cp,"\n")
    cat("Err:\n")
    print(err)
  }
  #Compute final p(class|pattern)
  pp<-drop(rep(1,nl)%*%pas)
  pas<-t(t(pas)/pp)  
  #Prettify some of the output
  colnames(err)<-c("p(false neg)","p(false pos)")
  rownames(err)<-paste("Item",item.labels)
  D<-mw.D(lat)
  if(is.matrix(D)&&(NROW(D)>0))
    rownames(D)<-c(letters,1:(NROW(D)-26))[1:NROW(D)]
  #Prepare to return the results
  out<-list()
  out$y<-y[y[,nmb+1]>0,]
  out$lik<-ll
  out$slik<-ll0
  out$nlik<-lln
  out$lrchisq<-2*(ll0-ll)
  out$df.model<-switch(err.mode,
    symhom=nl-df.recovered,
    asymhom=nl+1-df.recovered,
    symitem=nl+nmb-1-df.recovered,
    asymitem=nl+2*nmb-1-df.recovered
  )
  out$df.sat<-npat-1
  out$df.residual<-out$df.sat-out$df.model
  out$df.recovered<-df.recovered
  out$df.null<-nmb
  out$AIC<--2*ll+2*out$df.model
  out$BIC<--2*ll+log(n)*out$df.model
  out$sAIC<--2*ll0+2*(npat-1)
  out$sBIC<--2*ll0+log(n)*(npat-1)
  out$nAIC<--2*lln+2*nmb
  out$nBIC<--2*lln+log(n)*nmb
  out$entropy.total<-etot
  out$entropy.data<-edat
  out$sentropy.total<-etot0
  out$sentropy.data<-edat0
  out$nentropy.total<-etotn
  out$nentropy.data<-edatn
  out$classprob<-cp
  out$memb<-max.col(t(pas))
  out$niter<-niter
  out$err.mode<-err.mode
  out$err<-switch(err.mode,
    symhom=err[1,1],
    asymhom=err[1,],
    symitem=err[,1],
    asymitem=err
  )
  out$n<-n
  out$npat<-npat
  out$obspat<-length(yc)
  out$nmb<-nmb
  out$nl<-nl
  out$lattice<-lat
  out$D<-D
  out$labels<-item.labels
  out$null.probs<-prob.null
  #Return the result
  class(out)<-"mw"
  out
  }

[Package ldsa version 0.1-2 Index]