CODEWITHSUNDEEP
rfunctionargumentsstata
user25494
user25494

Reputation: 1371

Running existing function with non-default option

The code pasted below from ResourceSelection::hoslem.test performs a Hosmer and Lemeshow goodness of fit test. While investigating why the output that does not agree exactly with that performed by another software (Stata), I found that the difference relates to use of default R argument for the quantile function (type=7). I would like to use this function with a different default for calculation of quantiles (type=6).

FWIW, the reference to the 9 possible methods used by R can be found at:

https://www.amherst.edu/media/view/129116/original/Sample+Quantiles.pdf

The Stata manual for pctile refers to a default method and an 'altdef' method. I found it difficult to map these two methods to corresponding R types.

However, 

hoslem.test(yhat, y, type=6)

Produces:

> hl <- hoslem.test(y, yhat, type=6)
Error in hoslem.test(y, yhat, type = 6) : unused argument (type = 6)

Is there a way to run the function below with a non-default argument for the quantile function?

Ie. allows the following line adding ', type=6':

qq <- unique(quantile(yhat, probs = seq(0, 1, 1/g), type=6))

The function in question is:

> ResourceSelection::hoslem.test
function (x, y, g = 10) 
{
    DNAME <- paste(deparse(substitute(x)), deparse(substitute(y)), 
        sep = ", ")
    METHOD <- "Hosmer and Lemeshow goodness of fit (GOF) test"
    yhat <- y
    y <- x
    qq <- unique(quantile(yhat, probs = seq(0, 1, 1/g)))
    cutyhat <- cut(yhat, breaks = qq, include.lowest = TRUE)
    observed <- xtabs(cbind(y0 = 1 - y, y1 = y) ~ cutyhat)
    expected <- xtabs(cbind(yhat0 = 1 - yhat, yhat1 = yhat) ~ 
        cutyhat)
    chisq <- sum((observed - expected)^2/expected)
    PVAL = 1 - pchisq(chisq, g - 2)
    PARAMETER <- g - 2
    names(chisq) <- "X-squared"
    names(PARAMETER) <- "df"
    structure(list(statistic = chisq, parameter = PARAMETER, 
        p.value = PVAL, method = METHOD, data.name = DNAME, observed = observed, 
        expected = expected), class = "htest")
}

Upvotes: 4

Views: 180

Answers (3)

Vlo
Vlo

Reputation: 3188

We can modify pieces of functions. Look at the body of the function

as.list(body(hoslem.test))

See that the element we want to modify is the 6th element in the body

[[1]]
`{`

[[2]]
DNAME <- paste(deparse(substitute(x)), deparse(substitute(y)), 
    sep = ", ")

[[3]]
METHOD <- "Hosmer and Lemeshow goodness of fit (GOF) test"

[[4]]
yhat <- y

[[5]]
y <- x

[[6]]
qq <- unique(quantile(yhat, probs = seq(0, 1, 1/g)))

Modify the 6th element to what you want

body(hoslem.test)[[6]] = substitute(qq <- unique(quantile(yhat,
                                    probs = seq(0, 1, 1/g), type = 6)))

Upvotes: 2

user25494
user25494

Reputation: 1371

The two answers suggest a wrapper function to flexibly modify hoslem.test

myhoslem.test<-function(x, y, g = 10, mytype = 6){
  body(hoslem.test)[[6]] = substitute(qq <- unique(quantile(yhat,
                                              probs = seq(0, 1, 1/g), type = mytype))) 
  hoslem.test(x,y, g=10)
}

Upvotes: 0

jeremycg
jeremycg

Reputation: 24945

The easiest way would be to reenter the function as your own:

myhoslem.test<-function(x, y, g = 10, mytype = 6){
    DNAME <- paste(deparse(substitute(x)), deparse(substitute(y)), 
        sep = ", ")
    METHOD <- "Hosmer and Lemeshow goodness of fit (GOF) test"
    yhat <- y
    y <- x
    qq <- unique(quantile(yhat, probs = seq(0, 1, 1/g), type = mytype))
    cutyhat <- cut(yhat, breaks = qq, include.lowest = TRUE)
    observed <- xtabs(cbind(y0 = 1 - y, y1 = y) ~ cutyhat)
    expected <- xtabs(cbind(yhat0 = 1 - yhat, yhat1 = yhat) ~ 
        cutyhat)
    chisq <- sum((observed - expected)^2/expected)
    PVAL = 1 - pchisq(chisq, g - 2)
    PARAMETER <- g - 2
    names(chisq) <- "X-squared"
    names(PARAMETER) <- "df"
    structure(list(statistic = chisq, parameter = PARAMETER, 
        p.value = PVAL, method = METHOD, data.name = DNAME, observed = observed, 
        expected = expected), class = "htest")
}

The key change here is :

qq <- unique(quantile(yhat, probs = seq(0, 1, 1/g), type = mytype))

and allowing mytype as a argument to the function with default as 6

Upvotes: 1

Related Questions