0

I'm very newbie in R and i am creating a function to calculate the Furnival Index, trying to make it cleaner, where the user will only need to insert the adjusted model. I'm having some trouble to identify if the log transformation in the model, is a natural log or any other kind, because the index will change according to this information. So i thought to calculate this information using a=b^1/x, where "a" is the logarithm base, "x" and "b" are the formula information with/without log transformation respectively. But for doing that i need the original value from the model, because by using "model$model" i only get the logarithm value.

Here is what i have done so far:

furnival=function(object=NULL)
{
  w <- object$weights
  if(!is.null(object) && is.numeric(object))
    stop("'object' must be a formula")
  if(is.null(w <- object$weights)){
    w <- 1
  }else{
    w
  }
  if(length(grep("log", formula(object)))!=0){
    y <- as.numeric(as.matrix(object$model[1L]))
    modelValues <- object[Something to identify the original value]
    routine <- object$model == 1        
    if(any(routine))
       modelValues[!routine]
    modelValues <- sample(modelValues,1)
    a <- modelValues^(1/y)
    if(grep("log", formula(object))[1L]==2)
      y <- a^y
    if(a == exp(1)){ 
      df <- df.residual(object)
      MSE <- sum((residuals(object)^2)*w)
      index <- (exp(mean(log(y))))*(sqrt(MSE/df))
      return(index)
    }else{
      df <- df.residual(object)
      MSE <- sum((residuals(object)^2)*w)
      index <- (a^(mean(log(y,a))))*(sqrt(MSE/df))*(log(exp(1),a)^-1)
      return(index)
    }
  }
  else{
    df <- df.residual(object)
    MSE <- sum((residuals(object)^2)*w)
    index <- sqrt((MSE/df))
    return(index)
  }
}

If there is some way to do this, or even if there is a smarter way to make this function.

somoto
  • 27
  • 3
  • Can you also show how you are passing values to this function? Perhaps you could show sample models with different `log`s and show exactly what you want the desired behavior to be for each of those different cases. Also, it seems like some of those code might not really be necessary to make a *minimal* [reproducible example](http://stackoverflow.com/questions/5963269/how-to-make-a-great-r-reproducible-example). Try to focus on a very specific question/task and remove anything that's unnecessary to helping you solve your problem. – MrFlick Sep 07 '14 at 04:41
  • It's the first that i make a question here, i didn't know very well how to ask this. Sorry for the inconvenient. Some examples of what i am trying to do are like this: mod <- lm(log(z) ~ y) or even mod2 <- lm(log(z,10) ~ y) With that information i want to extract the log base. With the idea that i had, that is using a=b^1/x, the "a" will be the log base, the "b" will be the original value without the transformation, only "z" in the example and "x" the value transformed, in the example "log(z)". I achieved to obtain the "log(z)" value, but i didn't figure it out how to obtain only "z". – somoto Sep 07 '14 at 12:54
  • How about if you obtain predictions from the model at two different sets of predictor values (using `predict(model, newdata = data.frame(...))`, and compare them? By judiciously choosing those two predictor settings, it should be able to identify what base of logarithm was used. – Russ Lenth Sep 07 '14 at 13:28
  • @somoto Welcome to Stack Overflow. If people ask for more information, it's generally better to edit your original question to include the requested information or description. Then you can add a comment with their "@" username to notify them of the update. It's not a good idea to put important information in comments where it is hard to scan and format. – MrFlick Sep 07 '14 at 13:49

1 Answers1

0

If I isolate the part only where you try to determine the base of the log transformation of the response of your formula, this helper function should to that.

getresplogbase <- function(obj) {
    if(class(obj)=="lm") {
        obj = terms(obj)
    }
    stopifnot(is(obj,"formula"))
    rhs <- obj[[2]]
    if (is.recursive(rhs)) {
        if(rhs[[1]]==quote(log)) {
            if(length(rhs)==2) {
                return(exp(1))
            } else {
                return(eval(rhs[[3]], environment(obj)))
            }
        } else {
            stop("unable to parse:", deparse(rhs))
        }
    } else {
        NA
    }
}

For example, you can use it like

getresplogbase(y~x)
# [1] NA
getresplogbase(log(y)~x)
# [1] 2.718282
getresplogbase(log(y,10)~x)
# [1] 10
a<-2
getresplogbase(log(y,a)~x)
# [1] 2

You can also pass in lm() models

dd <- data.frame(y=runif(50,4,50)); dd$z=log(dd$y,2)+rnorm(50)
mod <- lm(log(z) ~ y, dd)
getresplogbase(mod)
# [1] 2.718282

All of this is done through a careful deparsing of the formula object that was used to fit the model.

MrFlick
  • 195,160
  • 17
  • 277
  • 295