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!
On Feb 27, 2013, at 1:40 PM, Millo Giovanni <[email protected]> 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: [email protected]
> [mailto:[email protected]] On Behalf
>> Of Simon Zehnder
>> Sent: Tuesday, February 26, 2013 2:53 AM
>> To: [email protected] 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.
______________________________________________
[email protected] mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
