[R] Count function calls

[Simon Zehnder]

Dear Giovanni,

apologize for this late reply! I was testing and reading a lot of stuff. I tried your suggestions and the problem of singularity in the regressor cross product vanishes when using the Group Mean function 'pgm' instead of 'pvcm'. Nevertheless, I found the collinearity in the regressors and 'pvcm' ran further, until I got another error. This time a little different but still directed towards the same area:

'Error in solve.default(crossprod(X[[i]])[!coefna[i, ], !coefna[i, ]])) :
	system is computationally singular: reciprocal condition number = 1.86131e-17'

Following the description of R's function 'rcond', the reciprocal condition number measures how close a matrix is to be rank deficient. So in this case one regressor crossproduct seems to be sufficiently close to singularity, such that inversion becomes impossible. 

I attached you a simple subsample. If you let run the following commands on it:

require(plm)
data <- read.csv(path, sep = ",", header = TRUE)
pdata <- plm.data(data, index = c("gvkey", "fyear"))
debug(plm:::pvcm)
model.pvcm <- pvcm(LOGDXSGA~LOGDSALE + DECRLOGDSALE + DECRLOGDSALEWDAVG, data = pdata, model = "random")
...... go until
Browse[2]> 
debug: ml <- split(data cond)

Then type:
A <- as.matrix(ml$'25062'[, 2:4])
XX <- t(A) %*% A
qr(XX)$rank
[1] 2                            Aha, the rank is only 2! 
rcond(XX) 
[1] 9.339817e-16    Alright, the matrix XX is pretty near to singularity, even it is not really singular!

Let us have a look at it:
A[, 2]/A[, 3]
.....                              A[, 2] is almost a multiple of A[, 3], if the latter is multiplied by -6.231979 or -6.231331 ..... at least it suffices to let the rank shrink to 2. 

The thing is, the same regression goes through in the Stata package. As I would say, I am an R fanatist, I would like to know, why it runs in Stata and not in R.....and that one can let it run in R as well. Different number formats/precision? Maybe it suffices in such a case to enforce full rank of the crossproduct by enforcing positive definitess, for instance via the function 'nearPD' in the R package 'Matrix':

Try 
library(Matrix)
qr(nearPD(XX)$mat)$rank
[1] 3                                                 Yey!..... :)

Tell me your opinion Giovanni! Give me a little help here please! 


Best 

Simon
 

P.S. Thank your for this fantastic package making panel data estimation possible in the R environment! 

-------------- next part --------------



On Feb 27, 2013, at 1:40 PM, Millo Giovanni <Giovanni_Millo at Generali.com> wrote:

> Hello.
> Another thing you may want to do depends on whether you are using
> model="within" (the default) or model="random".
> 
> In the first case, pvcm() estimates separate regressions, so you just
> need to loop lm() on individual indices to spot where it fails.
> 
> In the second case, what you may want to do is try a similar estimator:
> pmg(..., model="mg"), which is an unweighted version of Swamy's
> estimator in pvcm() and a simpler function to read and modify, possibly
> using the global assignment operator '<<-' as already suggested by
> William to output diagnostics.
> 
> Best,
> Giovanni
> 
>> -----Original Message-----
>> From: r-help-bounces at r-project.org
> [mailto:r-help-bounces at r-project.org] On Behalf
>> Of Simon Zehnder
>> Sent: Tuesday, February 26, 2013 2:53 AM
>> To: r-help at r-project.org help
>> Subject: [R] Count function calls
>> 
>> Dear R-users,
>> 
>> I have the following problem: I am running the function 'pvcm' from
> the 'plm' Panel Data
>> package. Inside this function 'solve' is called and gives for a
> certain individual data series
>> an exception because of singularity. I would like to know which
> individual data series
>> causes this error. I tried to debug it, but this is truly painful, as
> solve is called inside of
>> 'lapply' and there are over 5,000 individual data series in the panel.
>> 
>> Now, what I would like to do is to count the calls to 'solve' inside
> the function, so I can
>> see, where the function throws the exception. I tried to use 'trace'
> with a count variable,
>> but I have no clue how to define a global variable to be used by trace
> and updated at
>> every call.....is there another approach?
>> 
>> 
>> Best
>> 
>> Simon
>> 
>> ______________________________________________
> 
>  
> Ai sensi del D.Lgs. 196/2003 si precisa che le informazioni contenute in questo messaggio sono riservate ed a uso esclusivo del destinatario. Qualora il messaggio in parola Le fosse pervenuto per errore, La invitiamo ad eliminarlo senza copiarlo e a non inoltrarlo a terzi, dandocene gentilmente comunicazione. Grazie.
> 
> Pursuant to Legislative Decree No. 196/2003, you are hereby informed that this message contains confidential information intended only for the use of the addressee. If you are not the addressee, and have received this message by mistake, please delete it and immediately notify us. You may not copy or disseminate this message to anyone. Thank you.