423684_zh-CN

liliconcom

从二元图中提取因子范围

如何有效地提取该曲线及其 95% CI 与任意 Y 值的交点。

dt = Open( "$sample_data/Growth.jmp" );
biv = Bivariate(
 Y( :age ),
 X( :ratio ),
 Fit Line( {Line Color( {204, 121, 41} )}, save studentized residuals ),
 SendToReport(
  Dispatch(
   {},
   "Bivar Plot",
   FrameBox,
   {Grid Line Order( 3 ), Reference Line Order( 4 )}
  )
 )
);
biv << (Curve[1] << Save Studentized Residuals);

biv2 = Bivariate(
 Y( :Studentized Residuals age ),
 X( :ratio ),
 Fit Polynomial( 2, {Confid Shaded Fit( 1 ), Line Color( {212, 73, 88} )} ),
 SendToReport(
  Dispatch(
   {},
   "ratio",
   ScaleBox,
   {Min( 0.2796875 ), Max( 1.2666259765625 ), Inc( 0.2 ), Minor Ticks( 1 )}
  ),
  Dispatch(
   {},
   "Studentized Residuals age",
   ScaleBox,
   {Add Ref Line( 1.96, "Solid", "Black", "", 1 ),
   Add Ref Line( -1.96, "Solid", "Black", "", 1 )}
  )
 )
);

所以我想获得曲线及其 95%CI 交叉 +/- 1.96 的比率值（总共 6 个值）。 有没有比我能想到的蛮力更好的解决方案：将预测值和平均置信极限公式转储到表中，并将带有值的行添加到表中以查看它何时超过限制？

高级统计建模数据存取数据探索和可视化大规模定制回复：从二元图中提取因子范围

扩展于@Craige_Hales和@马克贝利的代码，您可以执行如下操作来找到总体最小值，然后将其用作每侧后续搜索的界限，这样您就不会那么依赖初始值。另外，如果您使用平方而不是绝对值，则可以消除角点。 但在实践中，我很幸运地使用了绝对值和其他函数，并在分析器的优化器中进行了步骤更改。

 

Names Default To Here( 1 );

// example data set
dt = Open( "$sample_data/Growth.jmp" );

// fit first-order polynomial model to response
biv = Bivariate(
 Y( :age ),
 X( :ratio ),
 Fit Line
);

// save residuals for next step
biv << (Curve[1] << Save Studentized Residuals);
biv << Close Window();

// fit second-order polynomial model to residuals
biv = Bivariate(
 Y( :Studentized Residuals age ),
 X( :ratio ),
 Fit Polynomial( 2, {Confid Shaded Fit( 1 ), Line Color( {212, 73, 88} )} )
);

// save formulas for lower / upper 95% confidence bound
biv << (Curve[1] << Mean Confidence Limit Formula);

// get formulas
fUpper = Column( dt, N Cols( dt ) ) << Get Formula;
fLower = Column( dt, N Cols( dt ) - 1 ) << Get Formula;

// replace predictor name with x
fsUpper = Substitute( Name Expr( fUpper ), Expr( :ratio ), Expr( x ) );
fsLower = Substitute( Name Expr( fLower ), Expr( :ratio ), Expr( x ) );

// initialize x, could get x that predicts mean residual = 0
x = 0.5;

// find the minimum of the whole expression, one solution will fall above and one below this value
fOverallMin = Name Expr( fsLower );
Minimize( fOverallMin, { x } );
xOverallMin = x;

// find where the value cross the x axis on each side of the overall min
x = xOverallMin - 1;
fMinCrossingTarget = Insert( Expr( Power() ), Insert( Expr( add(-1.96) ), Name Expr( fsLower ) ) );
Minimize( fMinCrossingTarget, { x(-1000, xOverallMin) } );
xMinLower = x;

x = xOverallMin + 1;
Minimize( fMinCrossingTarget, { x(xOverallMin, 1000) } );
xMinUpper = x;

x = xOverallMin - 1;
fMaxCrossingTarget = Insert( Expr( Power() ), Insert( Expr( add(-1.96) ), Name Expr( fsUpper ) ) );
Minimize( fMaxCrossingTarget, { x(-1000, xOverallMin) } );
xMaxLower = x;

x = xOverallMin + 1;
Minimize( fMaxCrossingTarget, { x(xOverallMin, 1000) } );
xMaxUpper = x;

show(xOverallMin, xMinLower, xMinUpper, xMaxLower, xMaxUpper);

 

回复：从二元图中提取因子范围

您可以调整因子网格（从红色三角形重置因子网格）以仅获得一个限制或另一个限制。

回复：从二元图中提取因子范围

我认为问题在于绝对值函数制作的曲线没有很好的导数，使得最小化可能失败。 我希望@ih这个想法会迫使最小化留在正确的邻域，但也许不是。

我现在没有时间去思考它；对于这个问题，似乎有一个标准的替代方案来解决这个问题......产生更平滑的功能。 （也许是平方？类似的东西。）

 

编辑：你的结果与我的一致，这就是为什么我有关于目视检查的评论。 abs 函数在函数中形成一个角以最小化，最小化函数使用该角附近的斜率，并且它期望答案的每一侧都有一个连续（二阶？）导数。 由于坡度即刻翻转真实答案的每一面，令人惊讶的是有时会得到接近的答案。 我认为日志中的消息可能是因为这个。

回复：从二元图中提取因子范围

我拼命地尝试想出一个起始点的算法，例如，我使用二次公式来计算与 1.96 交叉的预测曲线，即使你将该数字代入起始条件，你也会得到相当不可预测的行为：

 

因此，感谢您展示如何使用它，我确实从这个示例中学到了很多东西，但它还不够强大。

 

回复：从二元图中提取因子范围

@Craige_Hales ，为了解决初始值问题，如果您首先最小化整个函数，您可以使用该比率作为最小化搜索的限制，并保证找到每一侧的值。 对于这个问题，可能是 0.8：

 

Minimize( fa, { x([very small number],[min ratio - found from minimizing the function]) } );

//For this example -10 is arbitrary (could use -1E6 or something), and 0.8 is
//roughly the overall function minimum
Minimize( fa, { x(-10,0.8) } );

 

回复：从二元图中提取因子范围

微小变化

offset = 1.96;

// (I do not know how to pick these without a visual examination of the result.)
startingPoint = 1.1; // or 0.5 (trial and error. pick something close and check.)

// example data set
dt = Open( "$sample_data/Growth.jmp" );

// fit first-order polynomial model to response
biv = Bivariate( Y( :age ), X( :ratio ), Fit Line );

// save residuals for next step
biv << (Curve[1] << Save Studentized Residuals);
biv << Close Window();

// fit second-order polynomial model to residuals
biv = Bivariate( Y( :Studentized Residuals age ), X( :ratio ), Fit Polynomial( 2, {Confid Shaded Fit( 1 ), Line Color( {212, 73, 88} )} ) );

// save formulas for lower / upper 95% confidence bound
biv << (Curve[1] << save predicteds);
biv << (Curve[1] << Mean Confidence Limit Formula);
// get formulas, using order columns were saved
fupper = Column( dt, N Cols( dt ) ) << Get Formula;
flower = Column( dt, N Cols( dt ) - 1 ) << Get Formula;
fpredicted = Column( dt, N Cols( dt ) - 2 ) << Get Formula;

// replace predictor name with x
fsupper = Substitute( Name Expr( fupper ), Expr( :ratio ), Expr( x ) );
fslower = Substitute( Name Expr( flower ), Expr( :ratio ), Expr( x ) );
fspredicted = Substitute( Name Expr( fpredicted ), Expr( :ratio ), Expr( x ) );

fsupper2 = Substitute( Expr( a - offset ), Expr( a ), Name Expr( fsupper ) ); // find zero at offset 
fslower2 = Substitute( Expr( a - offset ), Expr( a ), Name Expr( fslower ) ); // find zero at offset 
fspredicted2 = Substitute( Expr( a - offset ), Expr( a ), Name Expr( fspredicted ) ); // find zero at offset 

// initialize x, could get x that predicts mean residual = 0
x = startingPoint;
// minimize expression that is Abs( residual ), or mean = 0
faupper = Insert( Expr( Abs() ), Name Expr( fsupper2 ) );
Minimize( faupper, {x} );
xupper = x;
Show( xupper, Eval( faupper ) );

// initialize x, could get x that predicts mean residual = 0
x = startingPoint;
// minimize expression that is Abs( residual ), or mean = 0
falower = Insert( Expr( Abs() ), Name Expr( fslower2 ) );
Minimize( falower, {x} );
xlower = x;
Show( xlower, Eval( falower ) );


// initialize x, could get x that predicts mean residual = 0
x = startingPoint;
// minimize expression that is Abs( residual ), or mean = 0
fapredicted = Insert( Expr( Abs() ), Name Expr( fspredicted2 ) );
Minimize( fapredicted, {x} );
xpredicted = x;
Show( xpredicted, Eval( fapredicted ) );

// add lines on plot
Report( biv )[axisbox( 2 )] << {Add Ref Line( xpredicted, "Solid", "Black", "P", 1 ), Add Ref Line( xupper, "Solid", "Black", "U", 1 ),
Add Ref Line( xlower, "Solid", "Black", "L", 1 )};
Report( biv )[axisbox( 1 )] << {Add Ref Line( offset, "Solid", "Black", "", 1 )};
Report( biv )[framebox( 1 )] << Frame Size( 752, 513 );

Low sideLow side低侧High sideHigh side高端

 

回复：从二元图中提取因子范围

您可以使用分析器中内置的优化器来执行此操作。设置响应限制以匹配目标，然后要查找上限和下限，只需使用重置因子网格来定义要搜索的值范围（高于或低于最小值）。

 

Profiler(
 Y( :Lower 95% Mean Studentized Residuals age ),
 Profiler(
  1,
  Desirability Functions( 1 ),
  Lower 95% Mean Studentized Residuals age <<
  Response Limits(
   {Lower( -0.04, 0.0183 ), Middle( 1.96, 1 ), Upper( 3.96, 0.0183 ),
   Goal( "Match Target" ), Importance( 1 )}
  ),
  Term Value( ratio( 0.8, Min( 0 ), Max( 0.8 ), Lock( 0 ), Show( 1 ) ) )
 ),
 SendToReport(
  Dispatch(
   {"Prediction Profiler"},
   "1",
   ScaleBox,
   {Min( 0 ), Max( 0.8 ), Inc( 0.2 ), Minor Ticks( 0 )}
  )
 )
) << Maximize Desirability;

Profiler(
 Y( :Lower 95% Mean Studentized Residuals age ),
 Profiler(
  1,
  Desirability Functions( 1 ),
  Lower 95% Mean Studentized Residuals age <<
  Response Limits(
   {Lower( -0.04, 0.0183 ), Middle( 1.96, 1 ), Upper( 3.96, 0.0183 ),
   Goal( "Match Target" ), Importance( 1 )}
  ),
  Term Value( ratio( 1.52, Min( 0.8 ), Max( 2 ), Lock( 0 ), Show( 1 ) ) )
 ),
 SendToReport(
  Dispatch(
   {"Prediction Profiler"},
   "1",
   ScaleBox,
   {Min( 0.8 ), Max( 2 ), Inc( 0.2 ), Minor Ticks( 0 )}
  )
 )
) << Maximize Desirability;

您还应该能够修改发送到最小化函数的表达式，只需减去 1.96 并最小化残差即可获得相同的结果。由于某种原因，即使在调整公差之后，我得到的值也略有不同（在探查器 v 1.03 中为 1.07，最小化）。

fa = Insert( Expr( Abs() ), parse( "abs( " || char( Name Expr( fs ) ) || " - 1.96 )" ) );

// find lower value
x = 0.2;
Minimize( fa, { x(0.2,0.8) } );
Show( x );

// find upper value
x = 1;
Minimize( fa, { x(0.8,1.5) } );
Show( x );

 

回复：从二元图中提取因子范围

是的，这确实适用于找到局部最小值，但仍然很难从那里到达下一个 1.96 的交叉点。 您仍然需要沿着曲线向上走，直到值超过 1.96。

 

为了暂时解决这个问题，我为其构建了一个函数，该函数只需要很少的步骤并在两个方向上输出，直到该值超过限制（或最大步骤数），但我预计 JMP 中有一个求解器函数可以实现更好的算法（就像非线性平台中内置的算法）。 我想现在我已经有了我需要的东西，但我总是愿意寻求更好的方法来做到这一点。

回复：从二元图中提取因子范围

是的，这个问题的问题在于，通常与目标值有多个交叉（一个在左侧，一个在右侧），并且使用意愿最大化器，您会随机获得两者之一。

我尝试了拟合模型下的逆预测，但不支持二次项。

回复：从二元图中提取因子范围

你的代码有点乱。 我清理了它并完成了这个例子，尽管我可能仍然在追求错误的目标。

 

Names Default To Here( 1 );

// example data set
dt = Open( "$sample_data/Growth.jmp" );

// fit first-order polynomial model to response
biv = Bivariate(
 Y( :age ),
 X( :ratio ),
 Fit Line
);

// save residuals for next step
biv << (Curve[1] << Save Studentized Residuals);
biv << Close Window();

// fit second-order polynomial model to residuals
biv = Bivariate(
 Y( :Studentized Residuals age ),
 X( :ratio ),
 Fit Polynomial( 2, {Confid Shaded Fit( 1 ), Line Color( {212, 73, 88} )} )
);

// save formulas for lower / upper 95% confidence bound
biv << (Curve[1] << Mean Confidence Limit Formula);

// get upper formula
f = Column( dt, N Cols( dt ) ) << Get Formula;

// replace predictor name with x
fs = Substitute( Name Expr( f ), Expr( :ratio ), Expr( x ) );

// initialize x, could get x that predicts mean residual = 0
x = 0.5;

// minimize expression that is Abs( residual ), or mean = 0
fa = Insert( Expr( Abs() ), Name Expr( fs ) );

Minimize( fa, { x } );
Show( x );

回复：从二元图中提取因子范围

我想知道另一种方法。 如果您将预测刻画器和匹配目标意愿函数与目标响应结合使用会怎样？ 然后您可以选择预测平均值和平均值的 95% 置信区间估计值。

 

我可能已经偏离了你最初的目标……

回复：从二元图中提取因子范围

@马克贝利谢谢，我尝试了一下，但我不清楚语法。

这是我所得到的，但现在我需要弄清楚如何为其提供正确的优化约束：

Names Default To Here( 1 );

dt = Open( "$sample_data/Growth.jmp" );
biv = Bivariate(
 Y( :age ),
 X( :ratio ),
 Fit Line( {Line Color( {204, 121, 41} )}, save studentized residuals ),
 SendToReport(
  Dispatch(
   {},
   "Bivar Plot",
   FrameBox,
   {Grid Line Order( 3 ), Reference Line Order( 4 )}
  )
 )
);
biv << (Curve[1] << Save Studentized Residuals);
biv << close window();

biv2 = Bivariate(
 Y( :Studentized Residuals age ),
 X( :ratio ),
 Fit Polynomial( 2, {Confid Shaded Fit( 1 ), Line Color( {212, 73, 88} )} ),
 SendToReport(
  Dispatch(
   {},
   "ratio",
   ScaleBox,
   {Min( 0.2796875 ), Max( 1.2666259765625 ), Inc( 0.2 ), Minor Ticks( 1 )}
  ),
  Dispatch(
   {},
   "Studentized Residuals age",
   ScaleBox,
   {Add Ref Line( 1.96, "Solid", "Black", "", 1 ),
   Add Ref Line( -1.96, "Solid", "Black", "", 1 )}
  )
 )
);
biv << (Curve[1] << Mean Confidence Limit Formula);
f =column(dt, n cols(dt)) << Get Formula;
fs = substitute(name expr(f), expr(:ratio), expr(x1));
/*A = [1 5];
b = [-1.96, 1.96];
minFun = Constrained Maximize(
 f,
 {x1( 0, 5 )},
 <<LESSTHANEQ></LESSTHANEQ>
minimize(name expr(f), {})

最后一部分不起作用。 感谢您的帮助！

回复：从二元图中提取因子范围

您也许可以使用 Minimize() 函数和置信边界表达式。 请参阅帮助优化方法。

回复：从图表中提取范围

我不知道直接的方法，但这里有一些想法。有些需要编写脚本，而且我怀疑对于小问题来说，没有一个会更快。

 

• 通过从双变量报告中的表中提取参数和标准误差，然后使用适当的公式创建列，将每行的公式保存到数据表中。然后使用分析器中内置的优化器来解决每个问题。
• 使用蛮力法保存相对较少的点，然后使用拟合模型对它们拟合一条新曲线，并且然后保存预测公式并使用优化器求解。
• 使用拟合模型进行拟合，您应该能够直接提取预测公式，但可能不能提取 CI。
• 使用强力方法并在最近的几个点之间进行插值，或者在最后两个点之间进行线性插值，或者对于该问题在最近的三个点之间进行二次插值。