Skip to content

Commit

Permalink
Various minor fixes for CRAN compliance
Browse files Browse the repository at this point in the history
  • Loading branch information
hawkrobe committed Sep 1, 2015
1 parent f439457 commit 467e10e
Show file tree
Hide file tree
Showing 4 changed files with 45 additions and 39 deletions.
5 changes: 3 additions & 2 deletions NAMESPACE
Original file line number Diff line number Diff line change
Expand Up @@ -8,5 +8,6 @@ export(GOF)
export(fit.grt)
export(mriTest)
export(riTest)
export(two_by_two_fit.grt)
export(two_by_two_plot.grt)
import(ellipse)
import(mnormt)
import(polycor)
65 changes: 35 additions & 30 deletions R/grt_base.R
Original file line number Diff line number Diff line change
@@ -1,3 +1,5 @@
#' @import ellipse mnormt polycor

# S3 grt object
#
# Constructor for a grt fit object, containing information about estimated parameters and likelihoods
Expand Down Expand Up @@ -247,7 +249,6 @@ mriTest <- function(x) {
}


#' @export
two_by_two_fit.grt <- function(freq, PS_x = FALSE, PS_y = FALSE, PI = 'none') {
if(!checkConfusionMatrix(freq)) return(FALSE); # Make sure confusion matrix valid
delta <- 1/10000; # Tolerance
Expand Down Expand Up @@ -365,7 +366,6 @@ two_by_two_fit.grt <- function(freq, PS_x = FALSE, PS_y = FALSE, PI = 'none') {



#' @export
two_by_two_plot.grt <- function(obj, xlab1, ylab1, level = .5) {
bin_width= .05; # determines smoothness of marginal plots
ex = .25 # determines relative size of main plot and marginal plots
Expand Down Expand Up @@ -461,9 +461,9 @@ isDefaultLabel <- function(label) {

#' Conduct goodness of fit tests
#'
#' Includes a number of common goodness of fit measures to compare different models.
#' Includes a number of common goodness of fit measures to compare different GRT models.
#'
#' @param bb a \code{grt} object
#' @param grtMod a \code{grt} object
#' @param teststat a string indicating which statistic to use in the test.
#' May be one of the following:
#' \itemize{
Expand All @@ -482,10 +482,10 @@ isDefaultLabel <- function(label) {
#' GOF(fit1, teststat = 'AIC')
#' GOF(fit2, teststat = 'AIC')
#' @export
GOF <- function(bb,teststat='X2',observed=NULL){
if (!identical(class(bb),'grt'))
GOF <- function(grtMod,teststat='X2',observed=NULL){
if (!identical(class(grtMod),'grt'))
stop('Argument must be object of class "grt"')
if (is.null(ff <- bb$fit) && is.null(observed))
if (is.null(ff <- grtMod$fit) && is.null(observed))
stop('Must have fitted model, observed frequencies, or both')
# added AIC, AICc, BIC 1.27.14
statlist <- c('X2','G2','AIC','AIC.c','BIC')
Expand All @@ -501,7 +501,7 @@ GOF <- function(bb,teststat='X2',observed=NULL){
nk <- apply(observed,3,sum)
ex <- array(dim=dim(observed))
for (k in 1:dim(observed)[3])
ex[,,k] <- bsd.freq(bb$rowcuts,bb$colcuts,bb$dists[k,],nk[k])
ex[,,k] <- bsd.freq(grtMod$rowcuts,grtMod$colcuts,grtMod$dists[k,],nk[k])
}
df <- length(observed) - dim(observed)[3] - df
if (test == 1){
Expand All @@ -514,31 +514,31 @@ GOF <- function(bb,teststat='X2',observed=NULL){
tstat <- 2*sum(observed*log(observed/ex))
}
if (test == 3){
map = bb$fit$map
map = grtMod$fit$map
k = 0
for(i in 1:ncol(map)){
k = k + sum(unique(map[,i])>0)
}
tstat <- 2*bb$fit$loglik + 2*k
tstat <- 2*grtMod$fit$loglik + 2*k
}
if (test == 4){
map = bb$fit$map
map = grtMod$fit$map
k = 0
for(i in 1:ncol(map)){
k = k + sum(unique(map[,i])>0)
}
n <- sum(observed)
aicStat <- 2*bb$fit$loglik + 2k
aicStat <- 2*grtMod$fit$loglik + 2*k
tstat <- aicStat + (2*k*(k+1))/(n-k-1)
}
if (test == 5){
map = bb$fit$map
map = grtMod$fit$map
k = 0
for(i in 1:ncol(map)){
k = k + sum(unique(map[,i])>0)
}
n <- sum(observed)
tstat <- 2*bb$fit$loglik + log(n)*k }
tstat <- 2*grtMod$fit$loglik + log(n)*k }
if (test < 3){
structure(tstat,names=teststat,df=df,class='bsdGOF')
}else{
Expand Down Expand Up @@ -599,16 +599,16 @@ n_by_n_fit.grt <- function (xx, pmap=NA, formula=x~., p0=NA, method=NA,
#print(pold)
#print(pmap)
#print(imap)
bb <- bsd.llike(pold,xx,pmap,imap,d.level=2,...)
#print(attr(bb,'ExpD2'))
#print(attr(bb,'gradient'))
dlt <- solve(attr(bb,'ExpD2'),attr(bb,'gradient'))
grtMod <- bsd.llike(pold,xx,pmap,imap,d.level=2,...)
#print(attr(grtMod,'ExpD2'))
#print(attr(grtMod,'gradient'))
dlt <- solve(attr(grtMod,'ExpD2'),attr(grtMod,'gradient'))
if (trace){
cat('Iteration number',iter,'\n')
cat('Row cutpoints', round(pold[imap$xi],4),'\n')
cat('Col cutpoints', round(pold[imap$eta],4),'\n')
print(round(bsd.map2array(pold,pmap,imap),4))
cat('Value',bb[1],'\n')
cat('Value',grtMod[1],'\n')
}
s <- stepsize
repeat{
Expand Down Expand Up @@ -670,7 +670,7 @@ n_by_n_fit.grt <- function (xx, pmap=NA, formula=x~., p0=NA, method=NA,
pp <- ffit$par
code <- ffit$convergence
iterations <- ffit$counts
bb <- bsd.llike(pp,xx,pmap,imap,d.level=2)
grtMod <- bsd.llike(pp,xx,pmap,imap,d.level=2)
found <- code < 4
}
if (!found) warning('Minimization routine encountered difficultites',
Expand All @@ -684,16 +684,16 @@ n_by_n_fit.grt <- function (xx, pmap=NA, formula=x~., p0=NA, method=NA,
muh <- array(dim=dim(xx),dimnames=dimnames(xx))
for (k in 1:KK) muh[,,k] <- bsd.freq(xi,eta,dists[k,],nk[k])
fit <- list(obs=xx,fitted=muh,estimate=pp,
expd2=attr(bb,'ExpD2'),map=pmap,
loglik=bb[1],code=code,iter=iterations)
expd2=attr(grtMod,'ExpD2'),map=pmap,
loglik=grtMod[1],code=code,iter=iterations)
if (verbose) {
if (found) cat('\nConvergence required',iterations,'iterations\n\n')
else cat (iterations, 'iterations used without reaching convergence\n\n')
cat('Parameter estimate vector\n'); # print(round(pp,3))
cat('Row cutpoints', round(xi,4),'\n')
cat('Col cutpoints', round(eta,4),'\n')
print(round(dists,4))
cat('Log likelihood =',bb[1],'\n')
cat('Log likelihood =',grtMod[1],'\n')
}
#print(dists);
output = grt(dists,fit=fit,rcuts = xi,ccuts = eta)
Expand Down Expand Up @@ -834,23 +834,28 @@ linear.hypothesis <- function(b,m,c=0,set='means'){

#' Compare nested GRT models
#'
#' Conducts a likelihood-ratio G-test on nested GRT models.
#' Conducts a likelihood-ratio G-test on nested GRT models. Currently only accepts pairs of nested models, not arbitrary sequences.
#'
#' @param model1 A fitted GRT model returned by fit.grt
#' @param model2 A larger GRT model, with model1 nested inside
#' @param object A fitted GRT model returned by fit.grt
#' @param ... A larger GRT model, with model1 nested inside
#' @export
anova.grt <- function(model1,model2){
anova.grt <- function(object, ...){
model1 <- object;
model2 <- list(...)[[1]];
g21 <- GOF(model1,teststat='G2')
df1 <- attr(g21,'df')
g22 <- GOF(model22,teststat='G2')
g22 <- GOF(model2,teststat='G2')
df2 <- attr(g22,'df')
# p.val added 1.24.14 -NHS
DG2 <- round(g21-g22,3)
ddf <- df1-df2
p.val <- round(pchisq(DG2,ddf,lower.tail=F),4)
table <- matrix(c(round(g21,3),round(g22,3),df1,df2,'',
DG2,'',ddf,'',p.val),2)
dimnames(table) <- list(c(substitute(model1),substitute(model2)),
# Want to get model var names passed in by user
arg <- deparse(substitute(object))
dots <- substitute(list(...))[-1]
modelNames <- c(arg, sapply(dots, deparse))
dimnames(table) <- list(modelNames,
c('G2', 'df', 'DG2', 'df','p-val'))
as.table(table)
}
Expand Down
6 changes: 3 additions & 3 deletions man/GOF.Rd
Original file line number Diff line number Diff line change
Expand Up @@ -4,10 +4,10 @@
\alias{GOF}
\title{Conduct goodness of fit tests}
\usage{
GOF(bb, teststat = "X2", observed = NULL)
GOF(grtMod, teststat = "X2", observed = NULL)
}
\arguments{
\item{bb}{a \code{grt} object}
\item{grtMod}{a \code{grt} object}

\item{teststat}{a string indicating which statistic to use in the test.
May be one of the following:
Expand All @@ -21,7 +21,7 @@ May be one of the following:
\item{observed}{optional, to provide a matrix of observed frequencies if no fit conducted.}
}
\description{
Includes a number of common goodness of fit measures to compare different models.
Includes a number of common goodness of fit measures to compare different GRT models.
}
\examples{
data(thomasA)
Expand Down
8 changes: 4 additions & 4 deletions man/anova.grt.Rd
Original file line number Diff line number Diff line change
Expand Up @@ -4,14 +4,14 @@
\alias{anova.grt}
\title{Compare nested GRT models}
\usage{
\method{anova}{grt}(model1, model2)
\method{anova}{grt}(object, ...)
}
\arguments{
\item{model1}{A fitted GRT model returned by fit.grt}
\item{object}{A fitted GRT model returned by fit.grt}

\item{model2}{A larger GRT model, with model1 nested inside}
\item{...}{A larger GRT model, with model1 nested inside}
}
\description{
Conducts a likelihood-ratio G-test on nested GRT models.
Conducts a likelihood-ratio G-test on nested GRT models. Currently only accepts pairs of nested models, not arbitrary sequences.
}

0 comments on commit 467e10e

Please sign in to comment.