6
$\begingroup$

I would like help with R code to estimate theta for the Inverse Hyperbolic Sine Transformation. This transformation is useful to transform skewed data that contain negative values or zeros.

There are a couple of related posts that discuss the IHS and suggest there is a maximum likelihood approach to estimating theta but I can't work out how to apply it. These related posts are:

Please see example code below. I've made skewed data by using the inverse of the IHS. Now, given the skewed data, and no prior knowledge of theta, how can I work out what theta should be? I would be most grateful for R code to undertake this analysis.

# Define the IHS transformation and its inverse
IHS <- function(x, theta){  # Inverse IHS transformation
  (1/theta)*asinh(theta * x)
}

Inv.IHS <- function(x, theta){  # IHS transformation
  (1/theta)*sinh(theta * x)
}

set.seed(1)
# generate some normal data
x <- rnorm(1000)
hist(x, breaks='FD')

# skew it by applying the Inverse of the IHS transformation
xt <- Inv.IHS(x, theta=2)
hist(xt, breaks='FD') # yep this is skewed.  How could we estimate theta?

usng inverse of the IHS to skew data

$\endgroup$
1
  • $\begingroup$ If you want to estimate it, maximum likelihood estimation, perhaps? $\endgroup$
    – Glen_b
    Dec 8, 2013 at 9:52

1 Answer 1

4
$\begingroup$

After a couple more days of thinking about the problem, I have two tentative answers.

  1. Select theta so that the transformed data is close to normal as measured by goodness of fit. For example choose theta to maximize the p-value of the Shapiro-Wik test.

    set.seed(1)
    x <- rnorm(1000)
    xt <- Inv.IHS(x, theta=2)
    
    Shapiro.test.pvalue <- function(theta, x){ 
      x <- IHS(x, theta) 
      shapiro.test(x)$p.value 
    }
    
    optimise(Shapiro.test.pvalue, lower=0.001, upper=50, x=xt, maximum=TRUE)  # 2.069838
    
  2. Maximum likelihood estimation of theta.

    Looking at the paper by Burbidge et al. I think the likelihood function for a single variable can be expressed as follows.

    IHS.loglik <- function(theta,x){
    
      IHS <- function(x, theta){  # function to IHS transform
      asinh(theta * x)/theta
    }
    
    n <- length(x)
    xt <- IHS(x, theta)
    
    log.lik <- -n*log(sum((xt - mean(xt))^2))- sum(log(1+theta^2*x^2))
    return(log.lik)
    }     
    
    # try this on our data
    optimise(IHS.loglik, lower=0.001, upper=50, x=xt, maximum=TRUE) # 2.0407
    

In both cases we get close to the expected theta, which is encouraging. But when I try the maximum likelihood approach on my real data it doesn't seem to give reasonable answers.

$\endgroup$
2
  • 1
    $\begingroup$ +1 (on both Q and A). ... "when I try the maximum likelihood approach on my real data it doesn't seem to give reasonable answers" - if the other approach does give reasonable answers, it might suggest a convergence problem with the maximization routine. What happens if you start it at the other estimate? $\endgroup$
    – Glen_b
    Dec 9, 2013 at 22:32
  • $\begingroup$ Unfortunately the shapiro version just works for small samples. It would be interesting to see an stable answer using ML. $\endgroup$
    – marbel
    Nov 6, 2014 at 20:20

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.