Bookmark and Share

Notice: On April 23, 2014, Statalist moved from an email list to a forum, based at statalist.org.


[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]

st: R: RE: friedman test


From   Maria Letizia Bacchi Reggiani <[email protected]>
To   "[email protected]" <[email protected]>
Subject   st: R: RE: friedman test
Date   Thu, 19 Jul 2012 15:38:14 +0000

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: [email protected] [mailto:[email protected]] Per conto di Feiveson, Alan H. (JSC-SK311)
Inviato: giovedì 19 luglio 2012 15.37
A: [email protected]
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: [email protected] [mailto:[email protected]] On Behalf Of Maria Letizia Bacchi Reggiani
Sent: Wednesday, July 18, 2012 9:41 AM
To: [email protected]
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/


© Copyright 1996–2018 StataCorp LLC   |   Terms of use   |   Privacy   |   Contact us   |   Site index