CODEWITHSUNDEEP
rfunctionstatisticsarguments
Emma Wiik
Emma Wiik

Reputation: 67

Sensitivity analysis with correlated predictors: Package pse LHS() with correlation structure does not converge

I noticed that function LHS() in package pse provides the argument opts = list() which accepts two arguments: COR and eps, the former being a correlation matrix and the latter the tolerance for the deviation between the prescribed correlation matrix and the result. Both feed down into a function called LHScorcorr()

I can run the function successfully for a high eps (e.g. 0.5, 1), however the default is 0.005. At low values of eps (< 0.1 or so) the function does not converge (see error message at bottom).

I would like to allow it to converge at lower values, and looking at LHScorcorr.R it does take a maxIt argument (maximum iterations). It is set to maxIt = 2*sqrt(dim(vars)[1]), however, unlike eps, cannot seemingly be fed into LHS(). According to ?LHS, it also doesn't take ... arguments.

1) Am I missing something or is there really no way of feeding a maxIt into the function call? 2) Is that even the cause of the problem I am encountering? 3) What are the "acceptable" ranges of eps to feed?

Here is a subset 'subdf' of my real-data data frame (I've tested my code on it): https://drive.google.com/file/d/0BwjVzXSG-JMJX2lZSG44TjNjb28/view?usp=sharing

And here is the function I have (currently still poorly commented since it is a modified version of one which takes more options): https://github.com/ewiik/flux/blob/master/functions/gasExchangeSensitivity.R

This is what I run:

datacorr <- cor(subdf, method = "spearman",use = "complete.obs")

factors <- c("Temperature", "Conductivity", "pH", "meanWindMS", "SalCalc", "TICumol", 
             "Pressure", "pco2atm")

distro <- c("qweibull", "qnorm", "qlogis", "qnorm", "qweibull", "qnorm", "qnorm", "qunif") 

props <- list( list(shape=4.02, scale=18.65), list(mean=1013, sd=499), 
               list(location=8.84,scale=0.31), list(mean=4.98, sd=0.83),
               list(shape=2.13, scale=0.68), list(mean=3821, sd=1068),
               list(mean=94.6, sd=0.17), list(min=356.9, max=402.2))


latincorr <- LHS(gasExchangeSens, factors = factors, N = 200, q = distro, q.arg = props, 
                 nboot = 200, opts = list(COR = datacorr, eps = 0.1))

latincorr <- at the above eps setting gives:

In internal.LHScorcorr(vars, COR, l = l, eps = eps, it = it + 1,  :
  LHScorcorr: correlation does not converge after maximum iterations

Upvotes: 2

Views: 397

Answers (1)

andrechalom
andrechalom

Reputation: 737

Package author here :)

That is certainly a mistake on our part, the user SHOULD be able to set maxIt explicitly and it solves the problem in your case. I'll try to release a new version soon with the fix. (Update: pse 0.4.5 submitted to CRAN. Should be online soon)

However, it is important to note that there is another problem that might happen in cases such as this. The problem is that it may be impossible to generate samples from a multivariate distribution with fixed parameters and correlation matrix. In some situations, all correlation matrices are admissible (for instance, if all your q distributions are qnorm, you can specify the correlation terms freely). However, more complicated distributions sometimes do not allow for a free specification of correlation terms. This is the case, for your data, in the qweibull and qlogis distributions. So it might happen that the program will not converge, no matter how large maxIt is set.

We wrote a small section about this problem here: http://arxiv.org/abs/1210.6278

A better (but maybe denser) mathematical background can be found here: http://arxiv.org/abs/1010.6118

Upvotes: 1

Related Questions