Notice: On April 23, 2014, Statalist moved from an email list to a forum, based at statalist.org.
| From | Maria Letizia Bacchi Reggiani <maria.bacchireggiani@unibo.it> |
| To | "statalist@hsphsun2.harvard.edu" <statalist@hsphsun2.harvard.edu> |
| Subject | st: R: RE: R: RE: friedman test |
| Date | Fri, 20 Jul 2012 10:55:50 +0000 |
Hi Al,
again thanks a lot for your program that's working very well now.
Actually some code lines had invalid carriage return.
I really appreciate your help.
Best
letizia
-----Messaggio originale-----
Da: owner-statalist@hsphsun2.harvard.edu [mailto:owner-statalist@hsphsun2.harvard.edu] Per conto di Feiveson, Alan H. (JSC-SK311)
Inviato: giovedì 19 luglio 2012 18.05
A: statalist@hsphsun2.harvard.edu
Oggetto: st: RE: R: RE: friedman test
Letizia - Make sure that all the lines in my embedded code have a proper carriage return. I noticed that in my original message the lines
tab rank,matcell(A)
matrix D=vecdiag(A*A')
matrix T=D*A
were all stuck out to the right of a comment line - but then in your current response it looks like you fixed it.
If you like I can send it to you as an attachment outside of statalist.
Al
Here is the code again for fried.ado:
=============================
program define fried
/* Does Friedman's inference on two-way layout on ranks
as in Hollander and Wolfe CH. 7
For example, use rank_ex.dta (base-rounding data - p.141)
fried y1 y2 y3
*/
version 7.0
syntax varlist [if] [in]
tokenize `varlist'
local nv: word count `varlist'
local file "$S_FN"
cap keep `in'
cap keep `if'
forv j=1(1)`nv' {
local name`j': word `j' of `varlist'
qui drop if `name`j''==.
}
qui save temp,replace
/* First get ranks _r`j' */
local j 0
while `j' < `nv' {
local j = `j'+1
local name`j': word `j' of `varlist'
qui cap drop _r`j'
qui cap gen int _r`j'=.
}
qui cap drop _T
qui gen _T=.
local N = _N
local i 0
while `i' < `N' {
local i = `i'+1
preserve
qui {
keep in `i'
keep `varlist'
xpose,clear
gen k=_n
sort v1
gen m=_n
sort k
keep m v1 /* at this point m is integer rank */
gen int ord=_n
sort v1
egen rank=mean(m),by(v1)
sort ord
keep rank
/* ------- Calculate T = tie correction factor ---- */
tab rank,matcell(A)
matrix D=vecdiag(A*A')
matrix T=D*A
/* ------------------------------------------------ */
xpose,clear
mkmat v1-v`nv',matrix(M)
restore
replace _T=T[1,1]-`nv' in `i'
local j 0
while `j' < `nv' {
local j = `j'+1
replace _r`j'=M[1,`j'] in `i'
}
}
}
format _r* %5.1f
scalar rbar=0
local j 0
while `j' < `nv' {
local j = `j'+1
qui summ _r`j'
scalar rbar=scalar(rbar)+r(mean)
}
scalar rbar=rbar/`nv'
scalar nrbar=`N'*scalar(rbar)
scalar ss=0
local j 0
while `j' < `nv' {
local j = `j'+1
qui summ _r`j'
scalar ss=scalar(ss)+( r(sum)-scalar(nrbar))^2
}
qui summ _T
scalar bot=`N'*`nv'*(`nv'+1)-r(sum)/(`nv'-1)
scalar top=12*ss
scalar Sp=scalar(top)/scalar(bot)
scalar P = chiprob(`nv'-1,Sp)
di " "
di in green "Friedman's Two-way Rank Test - `N' observations, `nv' variables"
di in green "---------------------------------------------------------------"
di in yellow "S' = ",%8.2f scalar(Sp)," P = ",%8.5f scalar(P)
scalar Pall=scalar(P)
/* Multiple comparisons */
qui cap gen _u=.
qui cap gen _v=.
qui cap gen _P=.
di " "
di in green "Multiple Comparisons"
di in green "------------------------------------------------------"
di in yellow " Var1 Var2 Ru - Rv P (2-sided)"
local nuv=`nv'*(`nv'-1)/2
local N=_N
if `N' < `nuv' {
qui set obs `nuv'
}
local uv=0
local u 0
while `u' < `nv' {
local u = `u'+1
qui summ _r`u'
scalar U=r(sum)
local v `u'
while `v' < `nv' {
local v = `v'+1
qui summ _r`v'
scalar V=r(sum)
scalar D=U-V
scalar W=abs(U-V)
scalar sdw=sqrt(_N*`nv'*(`nv'+1)/12)
local w = W/sdw
frange 1000 `nv' `w'
scalar P=1-r(P)
local uv=`uv'+1
qui replace _u=`u' in `uv'
qui replace _v=`v' in `uv'
qui replace _P=scalar(P) in `uv'
noi di _col(2) "`name`u''",_col(14) "`name`v''",%11.2f scalar(D),_col(32) %8.4f scalar(P)
}
}
use `file',clear
end
-----Original Message-----
From: owner-statalist@hsphsun2.harvard.edu [mailto:owner-statalist@hsphsun2.harvard.edu] On Behalf Of Maria Letizia Bacchi Reggiani
Sent: Thursday, July 19, 2012 10:38 AM
To: statalist@hsphsun2.harvard.edu
Subject: st: R: RE: friedman test
Al - Thanks for your program, I think this could be the solution for me but I got an error message ' matrix option not allowed' .
Maybe it refers to
mkmat v1-v`nv',matrix(M)
even though I don't understand why
thanks anyway
letizia
-----Messaggio originale-----
Da: owner-statalist@hsphsun2.harvard.edu [mailto:owner-statalist@hsphsun2.harvard.edu] Per conto di Feiveson, Alan H. (JSC-SK311)
Inviato: giovedì 19 luglio 2012 15.37
A: statalist@hsphsun2.harvard.edu
Oggetto: st: RE: friedman test
Letizia - Here's what I got using my own friedman test program. I have tested this program on the data in Hollander and Wolfe CH. 7 (base-rounding data - p.141) and it agrees.
The program assumes the data is in wide form (as you have it):
It also assumes your data is saved before you run it.
I ran this on your data without and with Test 1:
Hope this helps -
Al Feiveson
. fried test2 test3 test4
Friedman's Two-way Rank Test - 26 observations, 3 variables
---------------------------------------------------------------
S' = 46.58 P = 0.00000
Multiple Comparisons
------------------------------------------------------
Var1 Var2 Ru - Rv P (2-sided)
test2 test3 -19.00 0.0229
test2 test4 -47.00 0.0000
test3 test4 -28.00 0.0003
. fried test1 test2 test3 test4
Friedman's Two-way Rank Test - 26 observations, 4 variables
---------------------------------------------------------------
S' = 68.23 P = 0.00000
Multiple Comparisons
------------------------------------------------------
Var1 Var2 Ru - Rv P (2-sided)
test1 test2 -14.50 0.4031
test1 test3 -36.00 0.0006
test1 test4 -67.50 0.0000
test2 test3 -21.50 0.0958
test2 test4 -53.00 0.0000
test3 test4 -31.50 0.0040
==================================================================================
program define fried
/* Does Friedman's inference on two-way layout on ranks
as in Hollander and Wolfe CH. 7
For example, use rank_ex.dta (base-rounding data - p.141)
fried y1 y2 y3
*/
version 7.0
syntax varlist [if] [in]
tokenize `varlist'
local nv: word count `varlist'
local file "$S_FN"
cap keep `in'
cap keep `if'
forv j=1(1)`nv' {
local name`j': word `j' of `varlist'
qui drop if `name`j''==.
}
qui save temp,replace
/* First get ranks _r`j' */
local j 0
while `j' < `nv' {
local j = `j'+1
local name`j': word `j' of `varlist'
qui cap drop _r`j'
qui cap gen int _r`j'=.
}
qui cap drop _T
qui gen _T=.
local N = _N
local i 0
while `i' < `N' {
local i = `i'+1
preserve
qui {
keep in `i'
keep `varlist'
xpose,clear
gen k=_n
sort v1
gen m=_n
sort k
keep m v1 /* at this point m is integer rank */
gen int ord=_n
sort v1
egen rank=mean(m),by(v1)
sort ord
keep rank
/* ------- Calculate T = tie correction factor ---- */
tab rank,matcell(A)
matrix D=vecdiag(A*A')
matrix T=D*A
/* ------------------------------------------------ */
xpose,clear
mkmat v1-v`nv',matrix(M)
restore
replace _T=T[1,1]-`nv' in `i'
local j 0
while `j' < `nv' {
local j = `j'+1
replace _r`j'=M[1,`j'] in `i'
}
}
}
format _r* %5.1f
scalar rbar=0
local j 0
while `j' < `nv' {
local j = `j'+1
qui summ _r`j'
scalar rbar=scalar(rbar)+r(mean)
}
scalar rbar=rbar/`nv'
scalar nrbar=`N'*scalar(rbar)
scalar ss=0
local j 0
while `j' < `nv' {
local j = `j'+1
qui summ _r`j'
scalar ss=scalar(ss)+( r(sum)-scalar(nrbar))^2
}
qui summ _T
scalar bot=`N'*`nv'*(`nv'+1)-r(sum)/(`nv'-1)
scalar top=12*ss
scalar Sp=scalar(top)/scalar(bot)
scalar P = chiprob(`nv'-1,Sp)
di " "
di in green "Friedman's Two-way Rank Test - `N' observations, `nv' variables"
di in green "---------------------------------------------------------------"
di in yellow "S' = ",%8.2f scalar(Sp)," P = ",%8.5f scalar(P)
scalar Pall=scalar(P)
/* Multiple comparisons */
qui cap gen _u=.
qui cap gen _v=.
qui cap gen _P=.
di " "
di in green "Multiple Comparisons"
di in green "------------------------------------------------------"
di in yellow " Var1 Var2 Ru - Rv P (2-sided)"
local nuv=`nv'*(`nv'-1)/2
local N=_N
if `N' < `nuv' {
qui set obs `nuv'
}
local uv=0
local u 0
while `u' < `nv' {
local u = `u'+1
qui summ _r`u'
scalar U=r(sum)
local v `u'
while `v' < `nv' {
local v = `v'+1
qui summ _r`v'
scalar V=r(sum)
scalar D=U-V
scalar W=abs(U-V)
scalar sdw=sqrt(_N*`nv'*(`nv'+1)/12)
local w = W/sdw
frange 1000 `nv' `w'
scalar P=1-r(P)
local uv=`uv'+1
qui replace _u=`u' in `uv'
qui replace _v=`v' in `uv'
qui replace _P=scalar(P) in `uv'
noi di _col(2) "`name`u''",_col(14) "`name`v''",%11.2f scalar(D),_col(32) %8.4f scalar(P)
}
}
use `file',clear
end
================================================================================================
program define frange,rclass
version 7.0
/* frange.do - evaluates CDF of range of k N(0,1) distribution*/
args NIT k w
preserve
drop _all
qui {
range x -6 6 `NIT'
gen F=normprob(x)
gen Fw=normprob(`w'+x)
gen fx=normd(x)
gen h=fx*`k'*(Fw-F)^(`k'-1)
integ h x
}
restore
return scalar P=r(integral)
end
=======================================================================================================
*/
-----Original Message-----
From: owner-statalist@hsphsun2.harvard.edu [mailto:owner-statalist@hsphsun2.harvard.edu] On Behalf Of Maria Letizia Bacchi Reggiani
Sent: Wednesday, July 18, 2012 9:41 AM
To: statalist@hsphsun2.harvard.edu
Subject: st: friedman test
I ‘ve some problems performing Friedman test on STATA 11.
The STATA command help suggests to xpose, clear data and then to write friedman v?
In this way I obtain :
friedman v?
Friedman = 22.7333
Kendall = 0.8420
P-value = 0.0000
With the form v1-v26 I obtain as result:
friedman v1-v26
Friedman = 59.5731
Kendall = 0.7638
P-value = 0.0000
Performing the same Friedman test with R and with Primer of Biostatistics I obtain a third result:
Friedman = 68.233
P-value = 0.0000
I’m a bit confused, could anybody help me?
Thanks
Letizia Bacchi
PS here you find the data
test1 test2 test3 test4
0 8 72 900
0 2 48 155
0 25 45 350
0 15 35 900
0 170 720 900
100 200 250 2000
0 0 0 750
0 36 150 2300
0 0 0 550
0 30 300 950
0 60 0 950
0 0 36 170
0 0 0 400
0 0 0 850
0 0 0 60
0 50 200 200
0 87 240 240
0 20 500 800
0 0 10 550
0 200 250 450
0 0 400 850
0 0 0 750
0 0 100 750
0 0 350 600
0 0 15 850
0 40 400 480
LA RICERCA C’È E SI VEDE:
5 per mille all'Università di Bologna - C.F.: 80007010376
http://www.unibo.it/5permille
Questa informativa è inserita in automatico dal sistema al fine esclusivo della realizzazione dei fini istituzionali dell’ente.
*
* For searches and help try:
* http://www.stata.com/help.cgi?search
* http://www.stata.com/support/statalist/faq
* http://www.ats.ucla.edu/stat/stata/
*
* For searches and help try:
* http://www.stata.com/help.cgi?search
* http://www.stata.com/support/statalist/faq
* http://www.ats.ucla.edu/stat/stata/
LA RICERCA C’È E SI VEDE:
5 per mille all'Università di Bologna - C.F.: 80007010376
http://www.unibo.it/5permille
Questa informativa è inserita in automatico dal sistema al fine esclusivo della realizzazione dei fini istituzionali dell’ente.
*
* For searches and help try:
* http://www.stata.com/help.cgi?search
* http://www.stata.com/support/statalist/faq
* http://www.ats.ucla.edu/stat/stata/
*
* For searches and help try:
* http://www.stata.com/help.cgi?search
* http://www.stata.com/support/statalist/faq
* http://www.ats.ucla.edu/stat/stata/
LA RICERCA C’È E SI VEDE:
5 per mille all'Università di Bologna - C.F.: 80007010376
http://www.unibo.it/5permille
Questa informativa è inserita in automatico dal sistema al fine esclusivo della realizzazione dei fini istituzionali dell’ente.
*
* For searches and help try:
* http://www.stata.com/help.cgi?search
* http://www.stata.com/support/statalist/faq
* http://www.ats.ucla.edu/stat/stata/