Hola a todos, Quería proporcionar alguna resolución a este hilo cortesía de @sseligman en el soporte al cliente de SAS JMP (¡increíble!), que proporcionó un script .jsl que calcula los Intervalos de Confianza de Hodges-Lehmann fuera de la plataforma Fit Y by X.
Simplemente modificando el nivel alfa, uno puede calcular los límites de 1 cara correspondientes vs. JMPs por defecto los límites de 2 caras. Además, este script ofrece la ventaja de calcular estos límites no paramétricos tanto con y sin una llamada "corrección de continuidad" que JMP emplea de forma predeterminada, y SAS no. Esto podría ser particularmente útil para usted si, por ejemplo, tiene dos grupos, el Grupo 1 y el Grupo 2, y está calculando intervalos para la evaluación de riesgos a través de ellos con el fin de estimar qué cambio necesitaría aplicar a una especificación existente para el Grupo 1, de modo que sea igualmente conservador para el Grupo 2. Por supuesto, esto pone de manifiesto que existe una diferencia estadísticamente detectable en la diferencia de ubicación entre los Grupos (y que esa diferencia es de importancia práctica).
Para intervalos de 1 cara en general, reemplace "α / 2" con "α". Entonces, para el caso de 1 cara Confianza = 1 - α. Para el caso de 2 caras, Confianza = 1 - (α / 2).
Finalmente, si desea cambiar el orden de la diferencia tomada (y los cálculos subsiguientes sobre la diferencia de ubicación, incluidos los IC), simplemente debe cambiar "Ascendente" a "Descendente" en la siguiente sección del script:
// Pick up distinct values
sortx = Sort Ascending( x );
Aquí está el script completo, y también lo estoy incluyendo como un archivo adjunto.
// Main routine
MainNonParMCP=Function({},{default local},
// Column dialog to get data values
alpha=0.10;
{alpha,ng, xdistinct, x, y}=GetDataMatrix(alpha);
/* ---------------------------------------------------*/
/* for each pair */
OB1=NonparMCP(ng, xdistinct, x, y, alpha, method=1);
/* for all pairs (Steel-Dwass) */
OB2=NonparMCP(ng, xdistinct, x, y, alpha, method=2);
/* comparisons with the control */
OB3=NonparMCP(ng, xdistinct, x, y, alpha, method=3);
/* ---------------------------------------------------*/
/* Bonferroni for Joint Ranking */
/* for each pair */
OB4=BonferroniJointRanking(ng, xdistinct, x, y, alpha, method=1);
/* for all pairs */
OB5=BonferroniJointRanking(ng, xdistinct, x, y, alpha, method=2);
/* comparisons with the control */
OB6=BonferroniJointRanking(ng, xdistinct, x, y, alpha, method=3);
New Window("Result",
OutlineBox("Steel and Steel-Dwass",
OB1,OB2, OB3
),
OutlineBox("Bonferroni for Joint Ranking",
OB4,OB5, OB6
)
)
);
//
// Invoke column data dialog box, and
// get data values from a data table.
//
GetDataMatrix = Function(
{alpha},{default local},
dt = Current Data Table();
If( Is Empty( dt ),
dt = Open()
);
If( Is Scriptable( dt )==0,
Dialog( "Please open a data table before running this script" );
Throw();
);
// Column Dialog
dlg = Column Dialog(
"Select X and Y Column",
ycol = ColList( "Y, Continuous", DataType( numeric ), MinCol( 1 ), MaxCol( 1 ) ),
xcol = ColList( "X, Nominal", MinCol( 1 ), MaxCol( 1 ) ),
HList("Alpha : ",alpha =EditNumber(alpha))
);
// If canceled
If( Arg( dlg[N Arg( dlg )] ) != 1,
Dialog( "Canceled" );
Throw();
);
// Make Data
xx = Column( dt, dlg["xcol"] ) << GetAsMatrix;
yy = Column( dt, dlg["ycol"] ) << GetAsMatrix;
x = {};
y = [];
// Exclude missing values or excluded values
For( i = 1, i <= N Row( yy ), i++,
If( Excluded( Row State( i ) ) == 0 & Is Missing( xx[i] ) == 0 & Is Missing( yy[i] ) == 0,
Insert Into( x, xx[i] );
y |/= yy[i];
)
);
// Pick up distinct values
sortx = Sort Descending( x );
/*change "Descending" to "Ascending" in the above line of code to change the order of the groups*/
xdistinct = {};
Insert Into( xdistinct, sortx[1] );
For( i = 2, i <= N Items( sortx ), i++,
If( sortx[i] != sortx[i - 1],
Insert Into( xdistinct, sortx[i] )
)
);
ng = N Items( xdistinct );
If(IsMissing(dlg["alpha"])!=0,
alpha=dlg["alpha"]
);
Eval List( {alpha, ng, xdistinct, x, y} );
);
// Wilcoxon tests and Hodges-Lehman confidence interval for two samples
Wilcoxon_HL=Function({y1,y2,q},{default local},
n1=NRow(y1); /* Sample Size of Group 1 */
n2=NRow(y2); /* Sample Size of Group 2 */
n=n1+n2;
r=RankingTie(y1|/y2);
T1=Sum(r[1::n1]);
T2=n*(n+1)/2-T1;
diff=T2/n2-T1/n1;
stderr=Stddev(r)*Sqrt(1/n1+1/n2);
z=diff/stderr;
If(diff>0,
diffc=max(0,(T2-0.5)/n2-(T1+0.5)/n1),
diffc=min(0,(T2+0.5)/n2-(T1-0.5)/n1)
);
zc=diffc/stderr;
/* Hodges-Lehman: Point estimation and Confidence Interval */
dvec=J(n1*n2,1,.);
k=0;
For(i=1,i<=n1,i++,
For(j=1,j<=n2,j++,
k++;
dvec[k]=y2[j]-y1[i]
)
);
dvec=Sort Ascending(dvec);
/* Point Estimate */
If(Mod(n1*n2,2)==1,
HLEst= dvec[(n1*n2+1)/2] ,
HLEst= (dvec[(n1*n2)/2]+dvec[(n1*n2)/2+1])/2
);
/* Confidence interval without continuity correction */
C = Floor( n1*n2/2 - q*Sqrt(n1*n2*(n1+n2+1)/12));
If(C<1,
L=.;
H=.;
,
L=dvec[C];
H=dvec[n1*n2+1-C]
);
/* Confidence interval without continuity correction */
Cc = Floor( n1*n2/2 -0.5 - q*Sqrt(n1*n2*(n1+n2+1)/12));
If(Cc<1,
Lc=.;
Hc=.;
,
Lc=dvec[Cc];
Hc=dvec[n1*n2+1-Cc]
);
EvalList({diff,stderr,z,zc,HLEst,L,H,Lc,Hc});
);
//Nonparametric comparison
// Method==1 ... For each pair, without continuity correction
// Method==2 ... For all pairs, nonparametric Tukey based on pairwise ranking
// "Steel-Dwass test"
// Mehtod==3 ... Comparison with the control, nonparametric Dunnett based on pairwise ranking
// "Steel test"
NonparMCP=Function({ng, xdistinct, x, y, alpha, method},{default local},
Result=[];
g1Label={};
g2Label={};
n=J(NItems(xdistinct),1,0);
For(i=1,i<=NItems(x),i++,
For(j=1,j<=ng,j++,
If(x[i]==xdistinct[j],
n[j]++;
Break()
);
)
);
/* In order to perform, Steel test (nonparametric Dunnett),
we need to get sample size for each group*/
If(method==3,
ncont=n[1];
nother=Matrix(n[2::NRow(n)]);
rho=(1+ncont:/nother)^(-0.5);
);
/* Calculate Crit. value */
If(
method==1,
q=NormalQuantile(1-alpha/2),
method==2,
q=Tukey HSD Quantile(1-alpha,ng,.),
method==3,
q=DunnettQuantile(1-alpha,Nrow(nother),df=.,rho);
);
/* Loop Start */
For(i=1,i<=ng-1,i++,
If(method==3&i>=2,Break());
For(j=i+1,j<=ng,j++,
gval1=xdistinct[i];
gval2=xdistinct[j];
n1=n[i];
n2=n[j];
y1=[];
y2=[];
For(k=1,k<=NRow(y),k++,
If(
x[k]==gval1, y1|/=y[k],
x[k]==gval2, y2|/=y[k]
)
);
// Caclulate z-value and CI */
{diff,stderr,z,zc,HLEst,L,H,Lc,Hc}=Wilcoxon_HL(y1,y2,q);
// Calculate p-values from z
If(method==1,
p=2*NormalDistribution(-Abs(z));
pc=2*NormalDistribution(-Abs(zc));
,
method==2,
p=Tukey HSD P value(Abs(z),ng,.);
pc=Tukey HSD P Value(Abs(zc),ng,.);
,
method==3,
p=Dunnett P Value(Abs(z),ng-1,.,rho);
pc=Dunnett P Value(Abs(zc),ng-1,.,rho);
);
Insert Into(g1Label,Char(gval1));
Insert Into(g2Label,Char(gval2));
Result|/= (n1||n2||diff||stderr||z||p||zc||pc||HLEst||L||H||Lc||Hc);
)
);
If(
method==1, titleOB="Nonparametric Comparison for Each Pair",
method==2, titleOB="Nonparametric Comparison for All Pairs (Steel-Dwass)",
method==3, titleOB="Nonparametric Comparison for Each Pair (Steel)"
);
OB=OutlineBox(titleOB,
TableBox(
NumberColBox(Char(100*alpha)||"% Critical value",Matrix(q))
),
TableBox(
StringColBox("Group2",g2Label),
StringColBox("v.s. Group1",g1Label),
NumberColBox("n2", Matrix(Result[0,2])),
NumberColBox("n1", Matrix(Result[0,1])),
NumberColBox("Diff in Mean Rank", Matrix(Result[0,3])),
NumberColBox("Std Err", Matrix(Result[0,4])),
NumberColBox("Z", Matrix(Result[0,5])),
NumberColBox("p", Matrix(Result[0,6]),<<SET format="">=2,Break());
For(j=i+1,j<=ng,j++,
gval1=xdistinct[i];
gval2=xdistinct[j];
n1=n[i];
n2=n[j];
stderr=STotal*Sqrt(1/n1+1/n2);
diff=T[j]/n2-T[i]/n1;
z=diff/stderr;
If(diff>0,
diffc=Max(0,(T[j]-0.5)/n2-(T[i]+0.5)/n1),
diffc=Min(0,(T[j]+0.5)/n2-(T[i]-0.5)/n1)
);
zc=diffc/stderr;
If(method==1,
p=2*NormalDistribution(-Abs(z));
pc=2*NormalDistribution(-Abs(zc));
,
method==2,
p=Min(1,2*ng*(ng-1)/2*NormalDistribution(-Abs(z)));
pc=Min(1,2*ng*(ng-1)/2*NormalDistribution(-Abs(zc)));
,
method==3,
p=Min(1,2*(ng-1)*NormalDistribution(-Abs(z)));
pc=Min(1,2*(ng-1)*NormalDistribution(-Abs(zc)));
);
Insert Into(g1Label,Char(gval1));
Insert Into(g2Label,Char(gval2));
Result|/= (n1||n2||diff||stderr||z||p||zc||pc);
)
);
If(
method==1, titleOB="Nonparametric Comparison for Each Pair (without adjustment for Joint Ranking)",
method==2, titleOB="Nonparametric Comparison for All Pairs (Bonferonni for Joint Ranking)",
method==3, titleOB="Nonparametric Comparison for Each Pair (Bonferonni for Joint Ranking)"
);
OB=OutlineBox(titleOB,
TableBox(
StringColBox("Group2",g2Label),
StringColBox("v.s. Group1",g1Label),
NumberColBox("n2", Matrix(Result[0,2])),
NumberColBox("n1", Matrix(Result[0,1])),
NumberColBox("Diff in Mean Rank", Matrix(Result[0,3])),
NumberColBox("Std Err", Matrix(Result[0,4])),
NumberColBox("Z", Matrix(Result[0,5])),
NumberColBox("p", Matrix(Result[0,6]),<<SET format=""></SET></SET>
Entonces, los IC tampoco están asociados con la media. Están asociados con alguna medida de la mediana, según las fórmulas dadas por SAS:
Siguiendo a Hollander y Wolfe (1999), los límites de confianza superior e inferior asintóticos para el cambio de ubicación son
