Stata The Stata listserver
[Date Prev][Date Next][Thread Prev][Thread Next][Date index][Thread index]

RE: st: Kwallis in a loop [now p-value of the Kwallis]


From   "Herve STOLOWY" <stolowy@hec.fr>
To   <statalist@hsphsun2.harvard.edu>
Subject   RE: st: Kwallis in a loop [now p-value of the Kwallis]
Date   Sat, 12 Nov 2005 18:41:25 +0100

Dear Nick:

Sorry for this late reply. I tried your program and it works perfectly. I thank you very much.

Although I perfectly read what you wrote about -makematrix-, I still like this command and my level of programming is not good enough to write the program you designed for me.

So, I would have an additional question, hopefully more simple.

With a ranksum test, I learnt from prior postings how to get a matrix:

makematrix m1, from(r(z) 2*normprob(-abs(r(z)))) label : ranksum prov_ass, by( sodas2cl)

I woud like to do the same with -kwallis-. In your program, I found the following formula which seems to be the computation of the p-value of the kwallis: 

max(chiprob(r(df), r(chi2_adj) + 1e-20),.0001) 

Then, I wrote the following command with -makematrix-: 

makematrix m2, from(r(chi2_adj)   max(chiprob(r(df), r(chi2_adj) + 1e-20),.0001) ) label : kwallis prov_ass, by(sodas3cl)

Unfortunately, it does not work. I get the following error message:

invalid '.0001' 
r(198);

I obviously did something wrong.

Best regards

Hervé


***********************************************************
Professeur/Professor
Coordinateur du Département/Head of Department
HEC Paris
Département Comptabilité Contrôle de gestion / Dept of Accounting and Management Control
1, rue de la Liberation
78351 - Jouy-en-Josas
France
Tel: +33 1 39 67 94 42 - Fax: +33 1 39 67 70 86
stolowy@hec.fr
http://campus.hec.fr/profs/stolowy/perso/home.htm
>>> n.j.cox@durham.ac.uk 11/06/05 11:39 PM >>>
This all looks more complicated than needed, 
especially given the alternative of 
a few seconds' copy and paste. 

-tabstat- is not really designed for matrix 
output, as is illustrated by the fact that you 
have to massage things with -tabstatmat-, which 
I rather regret writing, as it tempts people into 
bad habits. Similarly, -makematrix- was a kind of 
experiment, which was interesting at least for me, 
but it is disconcerting whenever people use it 
and there are better tools available for what 
they are doing. It doesn't support -by:-, 
partly because other things do. In your case, 
you intend using it to put a scalar in a 1 X 1 
matrix, which is not necessary. -statsmat-, on 
the other hand, does seem closer to what you 
want to do here. 

The main part of that is wanting a matrix, 
because you want to use -mat2txt-, because you want 
a tab-delimited file for output. Or so 
I understand. 

Let's generalise to (1) several 
variables and keep (2) a single grouping 
variable, and write a program on the 
grounds that we can write something 
a little more widely applicable with 
not much more effort. This produces 
a matrix with Ns, means and sums 
of ranks and the Kruskal-Wallis P-value. 
I rather gave up on where the variable 
names should go. 

program pourherve 
	version 9 
	syntax varlist [if] [in], by(varname) matrix(str) 

	marksample touse 
	markout `touse' `by', strok 
	qui count if `touse' 
	if r(N) == 0 error 2000 

	tempname col work  
	qui tab `by' if `touse' 
	matrix `col' = J(`r(r)',1,.) 
	matrix colnames `col' = "P-value" 
	
	qui foreach v of local varlist { 
		tempvar rank 
		egen `rank' = rank(`v') if `touse' 
		statsmat `rank', by(`by') s(N sum mean) matrix(`work') 
		mat `work' = `work' , `col' 
		kwallis `v', by(`by') 
		matrix `work'[1,4] = ///
		max(chiprob(r(df), r(chi2_adj) + 1e-20),.0001) 
		matrix `matrix' = nullmat(`matrix') \ `work' 
		drop `rank' 
	}
end 

With your data 

pourherve index2 prov_ass, by(sodas3cl) matrix(m) 
mat2txt, matrix(m) saving(table5)

may help you along. Given other variables just
change the varlist. 

Nick 
n.j.cox@durham.ac.uk 

Herve STOLOWY
 
> Thanks again for your help. I applied your suggestion and it 
> works. The principle of the loop is fine.
> 
> I have a new question concerning how to get the p-value of 
> the kwallis test.
> 
> Here is my loop:
> 
> local action replace
> foreach var of varlist index2 prov_ass {
> egen rank_`var' = rank(`var') if sodas3cl < . 
> bysort sodas3cl : egen ranksum_`var' = sum(rank_`var') 
> by sodas3cl : egen rankmean_`var' = mean(rank_`var') 
> 
> statsmat `var', s(N) by(sodas3cl) mat(m1)
> tabstat rankmean_`var', by(sodas3cl) save
> tabstatmat m2
> matrix m3 = m1,m2
> statsmat rankmean_`var', s(N 0) mat(m4)
> matrix m5 = m3\m4
> makematrix m6, from(r(chi2_adj) formula to add) label : 
> kwallis `var', by(sodas3cl)
> matrix m7 = m5\m6
> mat2txt, matrix(m7) saving(table5) `action'
> local action append
> }
> 
> As the p-value of the test is not returned as a scalar, I 
> learnt from prior postings to the list that I should add the 
> relevant formula. For the Mann-Whitney U test, I know that 
> the formula is: 2*normprob(-abs(r(z))).
> 
> Unfortunately, I don't know the right formula for the Kwallis 
> test. I searched in the kwallis ado file but don't find it. I 
> certainly miss it.

*
*   For searches and help try:
*   http://www.stata.com/support/faqs/res/findit.html
*   http://www.stata.com/support/statalist/faq
*   http://www.ats.ucla.edu/stat/stata/


*
*   For searches and help try:
*   http://www.stata.com/support/faqs/res/findit.html
*   http://www.stata.com/support/statalist/faq
*   http://www.ats.ucla.edu/stat/stata/



© Copyright 1996–2014 StataCorp LP   |   Terms of use   |   Privacy   |   Contact us   |   What's new   |   Site index