help mfpi,help mfpi_plot-------------------------------------------------------------------------------

Title

mfpi-- Modelling interactions between categorical and continuous/categorical covariates

Syntax

mfpi[,mfpi_options]:regression_cmd[yvar1[yvar2]]xvarlist[if] [in] [weight] [,regression_cmd_options]

mfpi_plotvarname[if] [in] [,mfpi_plot_options]

mfpi_optionsDescription -------------------------------------------------------------------------adjust(xvarlist2)adjust linearly forxvarlist2allinclude out-of-sample observations in generated variablescenter(center_list)determine centering for FP-transformed variablesdetailgive additional details of fitted modelsdf(df_list)set up degrees of freedom for each predictor inxvarlistflex(#)define flexibility of main-effects and interaction modelsfp1(fp1_zvarlist)use fractional polynomials (FPs) of degree 1 for main effects and interactions withfp1_zvarlistfp2(fp2_zvarlist)use FPs of degree 2 for main effects and interactions withfp2_zvarlistgendiff(stubname)generate new variable(s) that contains the difference between the estimated functionsgenf(string)generate new variables that contain fitted functionslinear(linear_zvarlist)define list of variables whose linear interaction withtrt_varis to be investigatedmfpopts(mfp_options)additional options formfpnoscalingsuppress scaling of continuous variables in FP transformationsnocisuppress standard errors and confidence intervals forgendiff()andgenf()optionsselect(select_list)set significance levels for variable selection inxvarlistshowmodelshow variables and FP transformations selected fromxvarlistbymfptreatment(trt_var)required - defines the categorical 'treatment' variableregression_cmd_optionsoptions forregression_cmd-------------------------------------------------------------------------

regression_cmdmay beclogit,cnreg,glm,intreg,logistic,logit,mlogit,nbreg,ologit,oprobit,poisson,probit,qreg,regress,rreg,stcox,stcrreg,stpm2(if installed),streg,xtgeeand

varnameformfpi_plotmust be an existing variable. Typically,varnameis a member oflinear_zvarlist,fp1_zvarlistorfp2_zvarlist. See also the remark onWorking with composite variablesfor details of plotting with composite predictors, such as spline basis functions or fractional polynomials, which can be done using a special feature oflinear_zvarlists.

mfpi_plot_optionsDescription -------------------------------------------------------------------------stubname(stubname)required if thegendiff()option ofmfpiwas not specifiedvn(#)variable number inlinear(),fp1(), andfp2()level(#)select level oftrt_varto compare with base levelplot(plot)add other plots to generated graphgraph_optionsoptions forgraph twoway-------------------------------------------------------------------------

All weight types supported by

regression_cmdare allowed; see weight.

yvaris not allowed forstreg,stcoxandstpm2; for these commands, you must firststsetyour data.

Description

mfpiis designed to investigate the interaction of a categorical covariate (trt_var) with covariate(s) specified by any combination of thefp1(),fp2()andlinear()options. Whether or nottrt_varis specified as a factor variable, it is handled as one. Typically,trt_varis the treatment indicator variable in a randomized controlled trial of one or more treatments against a control arm, with the lowest value signifying the control arm.

mfpiadjusts interactions withtrt_varfor variables inxvarlist(which may be subjected to FP transformation selected bymfp) and/or linearly for variables inxvarlist2.

mfpi_plotproduces a treatment-effect plot derived from variables saved with thegendiff()option ofmfpi. A treatment-effect plot shows how the estimated treatment effect varies with a continuous covariate when an interaction between the two has been modelled.

Options formfpi

adjust(xvarlist2)includesxvarlist2as covariates with a linear functional form in all the fitted models.

allallows prediction (see thegendiff()andgenf()options) for all available cases, irrespective of exclusion byif,in, or weights.

center(center_list)determines the centering of FP-transformed variables. The most likely requirement is to suppress centering for all such variables, which is done withcenter(no). The default behavior is centering on the mean of each continuous predictor. For further details, see the description of thecenter()option offracpoly.

detailgives additional details of the fitted models.

df(df_list)sets up the degrees of freedom (df) for each predictor inxvarlist. The df (not counting the regression constant,_cons) are twice the degree of the FP, so, for example, a member ofxvarlistfitted as a second-degree FP (FP2) has 4 df. The first item indf_listcan be either#orvarlist:#. Subsequent items must bevarlist:#. Items are separated by commas, andvarlistis specified in the usual way for variables. With the first type of item, the df for all predictors are taken to be#. With the second type of item, all members ofvarlist(which must be a subset ofxvarlist) have#df.The default df for a predictor of type

varlistspecified inxvarlistbut not indf_listare assigned according to the number of distinct (unique) values of the predictor, as follows:--------------------------------------------------------------- No. of distinct values Default df --------------------------------------------------------------- 1 (invalid -- covariate has variance 0) 2-3 1 4-5 min(2,

dfdefault(#)) >=6dfdefault(#)---------------------------------------------------------------dfdefault(#)is an option ofmfp(see themfpopts()option andmfp); its default#is 4, meaning an FP2 function.Example:

df(4)All variables have 4 df.Example:

df(2, weight displ:4)weightanddisplhave 4 df, and all other variables have 2 df.Example:

df(weight displ:4, mpg:2)weightanddisplhave 4 df,mpghas 2 df, and all other variables have the default of 1 df.Example:

df(weight displ:4, 2)All variables have 2 df because the final 2 overrides the earlier 4.

flex(#)defines the flexibility of the main-effects and interaction models;#= 1 is the least flexible and#= 4 is the most flexible. The default isflex(1). SeeRemarksfor further details.

fp1(fp1_zvarlist)defines a list of continuous variables whose interactions withtrt_varare to be investigated by fitting FP functions of degree 1 (i.e., FP1 functions) to each member offp1_zvarlistin turn, at each level oftrt_var. Also seeflex().

fp2(fp2_zvarlist)defines a list of continuous variables whose interactions withtrt_varare to be investigated by fitting FP functions of degree 2 (i.e., FP2 functions) to each member offp2_zvarlistin turn, at each level oftrt_var. Also seeflex().

gendiff(stubname)generates a new variable(s) calledstubname#_jthat contains fj-f0, the difference between the estimated functions (at level j minus level 0 oftrt_var), for the#th member of the list composed oflinear_zvarlist,fp1_zvarlist, andfp2_zvarlist, in that order. The difference fj-f0 is an estimate of the covariate-specific effect of level j compared with level 0 (e.g., the covariate-specific treatment effect).gendiff()also creates new variables calledstubname#s_j, which contains the standard error of fj-f0, andstubname#lb_jandstubname#ub_j, the lower and upper 95% confidence limits, thus providing all the quantities necessary for a treatment-effect plot.

genf(stubname)generates new variables calledstubname#_0,stubname#_1, etc., that contain the fitted functions at levels 0, 1, etc., oftrt_var, respectively, for the#th member of the list composed oflinear_zvarlist,fp1_zvarlist, andfp2_zvarlist, in that order. For variables infp1_zvarlistandfp2_zvarlist, the same FP transformation is used at each level oftrt_var. The estimated function at level 0 oftrt_varis adjusted to have mean 0.

linear(linear_zvarlist)defines a list of variables whose linear interactions withtrt_varare to be investigated. If a categorical variable inlinear_zvarlisthas more than 2 levels, the necessary dummy variables must be created and placed between parentheses to indicate that they should be tested together. More properly,linear_zvarlistis a list of variables of which some may be dummies for categorical variables.For example,

linear((who2 who3))would bindwho2andwho3together to create a predictor with 2 df for its main effect.For further information, see the remarks in

Working with compositevariables.

mfpopts(mfp_options)suppliesmfpoptions tomfpifor the creation of the adjustment model fromxvarlist.

noscalingsuppresses scaling of all the continuous variables that are subject to FP transformation. The default is automatic application of scale factors; the behavior can be turned off for all such variables by usingnoscaling. For further details, see the description of the same option infracpolyandfracgen.

select(select_list)sets the nominal P-values (significance levels) for variable selection amongxvarlistby backward elimination. See theselect()option ofmfpfor further details. A typical usage isselect(0.05), which selects all variables inxvarlistthat are significant at the 5% level according to a backward stepwise algorithm.Example:

select(0.05)All variables have nominal p-value 5 percent.Example:

select(0.05, weight:1)All variables exceptweighthave nominal p-value 5 percent;weightis forced into the model.Example:

select(a (b c):0.05)All variables excepta,b, andcare forced into the model.bandcare tested jointly with 2 df at the 5 percent level, andais tested singly at the 5 percent level.

showmodelshows the variables selected fromxvarlistbymfp, together with their FP powers, where relevant. It also shows the interaction model fitted last.showmodelis a concise alternative todetail.

treatment(trt_var)is required and defines the factor or categorical variable whose interactions with variables specified inlinear(),fp1(), orfp2()are of interest.trt_varmust have at least two distinct, nonmissing values, but the codes are arbitrary because they are mapped to 0, 1, 2, etc., internally. The category corresponding to the lowest value is taken as the reference category (level 0).

regression_cmd_optionsare any options forregression_cmd.

Options formfpi_plot

stubname(stubname)is required if thegendiff()option ofmfpiwas not specified. If thegendiff()option was specified, thenstubnameis identified automatically. Thestubname()option allows you to produce a treatment-effect plot for a variable whosestubnamewas used in an earlier run ofmfpi.

vn(#)identifies the variable inlinear(),fp1(), andfp2(), in that order, whose associated treatment effect is to be plotted. Only one variable can be plotted in a given call tomfpi_plot. For example, iflinear(x1 x2)fp1(x1 x3)fp2(x3)was specified, the list of variables isx1 x2 x1 x3 x3, and the total number of variables is 5 (not 3). Thus#would be an integer between 1 and 5. The default isvn(1).

level(#)defines the level oftrt_varyou wish to compare with the base level. Levels oftrt_varare coded 0, 1, 2, etc., with the base (reference) level being 0. These coded levels may or may not coincide with the actual values oftrt_varin the data. For example,level(1)denotes the second lowest value oftrt_var, which could be any positive integer (not necessarily 1). The default#is 1, meaning compare the second level oftrt_varwith the first level. Specifying a non-existent level raises a "variable ... not found" type of error message.

plot(plot)provides a way to add other plots to the generated graph; see [G]addplot option.

graph_optionsare any options ofgraph twoway, such asxtitle()andytitle().

Remarks

MethodologyThe algorithm used in

mfpican be summarized as follows. To fix ideas, it is assumed thattrt_varis a binary treatment indicator variable in a randomized controlled clinical trial with two parallel groups.trt_vartherefore has two levels, arbitrarily labeled 0 and 1.We first describe the mfpi algorithm for the default flexibility setting,

flex(1), i.e., the least flexible case.Let Z be the continuous or binary covariate of immediate interest, and let T be the binary treatment variable with values of 0 and 1. Let X be a vector of potentially prognostic factors. In

mfpi, Z is a member oflinear_zvarlist,fp1_zvarlist, orfp2_zvarlist; T istrt_var; and X isxvarlist. The basic idea is to model the relationship between the outcome and Z at each level of T, adjusting for selected members of X, the latter variables transformed by FP if appropriate. A test for interaction is performed on regression coefficients at the final step. The procedure for a Z infp2_zvarlist(see thefp2()option) is as follows:1. Apply the MFP algorithm to X with a p-value threshold of a* for selecting variables and FP transformations. Let X* be the resulting covariate vector, which we will call the adjustment model. X* can include (transformed) variables in X selected by the MFP algorithm. If all variables in X are uninfluential, X* may even be null. No covariate adjustment is then performed. In some cases, X* can be formulated from subject-matter knowledge, avoiding data-driven searching.

2. Find by maximum likelihood the best-fitting FP2 powers (p1,p2) for Z, adjusting for X*.

3. For groups j = 0, 1 and powers p_i (i = 1, 2), define new predictors Z_ji = Z^p_i if T = j, Z_ji = 0 otherwise.

4. The test of T x Z interaction is a likelihood-ratio test between the nested models T, Z_01, Z_02, Z_11, Z_12, X* and T, Z^p_1, Z^p_2, X*. The difference in deviance is compared with chi-squared on 2 df.

5. If an interaction is not found, Z is regarded as a potential prognostic factor only. To investigate if an FP2 function is still needed for Z, the final model is chosen by repeating step 1 but including Z as a potential prognostic factor.

The procedure for a Z in

fp1_zvarlist, i.e., using FP1 rather than FP2 functions of Z, is in principle identical. The interaction is tested on 1 df.If Z is binary or is modeled as linear (i.e., is a member of

linear_zvarlist), the approach reduces to the usual procedure to estimate and to test for an interaction adjusted for covariates X*.Extension to cases where T has more than 2 levels follows immediately.

For further details of the MFP algorithm used to construct X*, see Royston and Ambler (1999), Sauerbrei and Royston (1999), and

mfp. Further information on the method and applications is given by Royston and Sauerbrei (2004).

Other flexibility valuesThe covariates X are treated in the same way for all values of

flex(). Only the main effect of Z and its interaction with T are dealt with differently, as follows.With

flex(2), the best-fitting FP powers for Z are found within the context of interaction models. The powers are constrained to be the same for all levels of T. The FP powers for the main effect of Z are taken to be the same as for the Z x T interaction.

flex(3)is the same asflex(2), except that the FP powers for the main effect of Z are determined independently of the interaction, exactly as forflex(1).

flex(4), the most flexible option, calls for FP powers to be estimated for the main effect of Z and also for all levels of T in the interaction context.

Hypothesis testingThe main-effects and interaction models for

flex(1)are nested, so standard hypothesis testing may be used. Simulation studies (as yet unpublished) have confirmed that the significance level offlex(1)is close to its nominal value.

flex(2)uses the same FP function for the interaction as for the main effect. At first glance, a standard hypothesis test appears appropriate. However, no allowance has been made for degrees of freedom used to estimate FP powers across the levels of the categorical variable in the interaction model, so the p-values from interaction tests should be regarded as indicative, not definitive.Similar comments apply to

flex(3), except that here the main-effects and interaction models are non-nested (because different FP models may be used for the main effects and interaction). P-values should again be taken as indicative rather than definitive.The extra flexibility provided by

flex(4)is reflected in its additional degrees of freedom. Its significance level is currently under investigation in simulation studies.

Working with composite variablesA useful feature of

mfpiis the ability to model, and if appropriate plot, functions of predictors represented by one or more individual but linked variables. For example, a restricted cubic regression spline with 4 knots, 2 interior and 2 boundary, is represented by 3 basis functions. Regression is performed with these basis functions as predictors. Such a set of predictors 'belongs' to just one original continuous predictor and must be handled appropriately. This is done by surrounding in parentheses (()) the linked variables in thelinear()option.For example, suppose that

age_1 age_2 age_3were restricted cubic spline basis variables derived from a predictor calledage, e.g. using the Stata commandmkspline. Analysis of interaction between treatment and the spline function forageis specified by the optionlinear((age_ age_2age_3)):

. mkspline age_ = age, cubic nknots(4)

. mfpi, select(0.05) treatment(treat) linear((age_1 age_2 age_3)): stcoxz1 z2 z3When

treatis binary, the interaction in the above example has 3 d.f.A nice feature of this approach is that the treatment-effect plot created by

mfpi_plotworks happily with the original variable. Continuing theageexample above, we specify thegendiff()option ofmfpito create the variables necessary for the treatment-effect plot and then tellmfi_plotto plot the age-related treatment difference and its 95% pointwise confidence interval againstageas follows:

. mfpi, select(0.05) treatment(treat) linear((age_1 age_2 age_3))gendiff(d): stcox z1 z2 z3

. mfpi_plot ageThe same approach works with categorical covariates with several levels, represented by dummy variables. Alternatively, such covariates can be defined as factor variables, for example,

linear(i.group). This avoids the need to create dummy variables. Factor variables are not allowed to be combined into groups using parentheses.

Note on the base level

mfpiassumes that the lowest value oftrt_varrepresents the base level (reference value). If you want a different value to be the base level, you have to recodetrt_varbefore runningmfpi. For example

. recode trt (1=3) (3=1)would reverse the ordering of the levels 1, 2, 3 of the variable

trt. Seerecodefor details.

WarningExcitement engendered by the discovery of interactions should be tempered by proper caution. There is a significant problem of multiplicity and of potential instability of the FP functions found. See Royston and Sauerbrei (2004) for a more detailed discussion of this issue. It is prudent to check whether the functions in each treatment group suggest a cutpoint(s) that can reproduce a difference in response; whether the interaction is present in univariate as well as multivariable models; and whether it is robust to the choice of adjustment model.

Examples

. mfpi, select(0.05) showmodel treatment(i.treat) linear(c1 (c2 c3))fp2(x1 x2): stcox x1 x2 c1 c2 c3

. mfpi, select(1) degree(1) dfdefault(2) treatment(i.treat) linear(c1 (c2c3)) fp1(x1 x2) fp2(x1 x2) genf(f) gendiff(d): logit y x1 x2 c1 c2 c3

. mfpi, treatment(i.treat) linear(i.group x): regress y

. mkspline age_ = age, cubic nknots(4)

. mfpi, select(0.05) treatment(i.treat) linear((age_1 age_2 age_3))gendiff(d): stcox z1 z2 z3

. mfpi_plot age

Saved resultsIn addition to what

regression_cmdsaves ine(),mfpisaves the some or all of the following inr():Scalars

r(nxvar)number of predictors inxvarlistr(nxf)number of predictors inxvarlistselected bymfpr(totdflin)total degrees of freedom for linear interaction modelr(totdffp1)total degrees of freedom for FP1 interaction modelr(totdffp2)total degrees of freedom for FP2 interaction modelr(devlin)deviance for linear interaction modelr(devfp1)deviance for FP1 interaction modelr(devfp2)deviance for FP2 interaction modelr(aiclin)AIC for linear interaction modelr(aicfp1)AIC for FP1 interaction modelr(aicfp2)AIC for FP2 interaction modelr(chi2lin)chi-squared value for testing linear interaction modelr(chi2fp1)chi-squared value for testing FP1 interaction modelr(chi2fp2)chi-squared value for testing FP2 interaction modelr(Plin)p-value for testing linear interaction modelr(Pfp1)p-value for testing FP1 interaction modelr(Pfp2)p-value for testing FP2 interaction modelMacros

r(varlist)xvarlistr(adjust)xvarlist2r(treatment)trt_varr(Linear)list of variables inlinear()r(Fp1)list of variables infp1()r(Fp2)list of variables infp2()r(z)z, the final continuous variable processedr(powmain)FP powers for z, main effectr(powint0)FP powers for z at level 0 oftrt_varr(powint`j')FP powers for z at level `j' oftrt_var(`j' > 0)

AuthorPatrick Royston MRC Clinical Trials Unit London, UK pr@ctu.mrc.ac.uk

ReferencesRoyston, P., and G. Ambler. 1999. sg112.1: Nonlinear regression models involving power or exponential functions of covariates: Update.

Stata Technical Bulletin50: 26. Reprinted inStata TechnicalBulletin Reprints, vol. 9, p. 180. College Station, TX: Stata Press.Royston, P., and W. Sauerbrei. 2004. A new approach to modelling interactions between treatment and continuous covariates in clinical trials by using fractional polynomials.

Statistics in Medicine23: 2509-2525.Royston, P., and W. Sauerbrei. 2009. Two techniques for investigating interactions between treatment and continuous covariates in clinical trials.

Stata Journal9 (2): 230-251.Sauerbrei, W., and P. Royston. 1999. Building multivariable prognostic and diagnostic models: Transformation of the predictors by using fractional polynomials.

Journal of the Royal Statistical Society,Series A162: 71-94.

Also seeArticle:

Stata Journal, volume 9, number 2: st0164Manual:

[R] fracpoly,[R] mfpOnline:

mfp