create partial dependence plot (pdp) and individual conditional expectation (ice) plots -凯发k8网页登录

train a regression tree using the carsmall data set, and create a pdp that shows the relationship between a feature and the predicted responses in the trained regression tree.

load the carsmall data set.

load carsmall

specify weight, cylinders, and horsepower as the predictor variables (x), and mpg as the response variable (y).

x = [weight,cylinders,horsepower];
y = mpg;

train a regression tree using x and y.

mdl = fitrtree(x,y);

view a graphical display of the trained regression tree.

view(mdl,"mode","graph")

create a pdp of the first predictor variable, weight.

plotpartialdependence(mdl,1)

the plotted line represents averaged partial relationships between weight (labeled as x1) and mpg (labeled as y) in the trained regression tree mdl. the x-axis minor ticks represent the unique values in x1.

the regression tree viewer shows that the first decision is whether x1 is smaller than 3085.5. the pdp also shows a large change near x1 = 3085.5. the tree viewer visualizes each decision at each node based on predictor variables. you can find several nodes split based on the values of x1, but determining the dependence of y on x1 is not easy. however, the plotpartialdependence plots average predicted responses against x1, so you can clearly see the partial dependence of y on x1.

the labels x1 and y are the default values of the predictor names and the response name. you can modify these names by specifying the name-value arguments predictornames and responsename when you train mdl using fitrtree. you can also modify axis labels by using the and functions.

train a naive bayes classification model with the fisheriris data set, and create a pdp that shows the relationship between the predictor variable and the predicted scores (posterior probabilities) for multiple classes.

load the fisheriris data set, which contains species (species) and measurements (meas) on sepal length, sepal width, petal length, and petal width for 150 iris specimens. the data set contains 50 specimens from each of three species: setosa, versicolor, and virginica.

load fisheriris

train a naive bayes classification model with species as the response and meas as predictors.

mdl = fitcnb(meas,species);

create a pdp of the scores predicted by mdl for all three classes of species against the third predictor variable x3. specify the class labels by using the classnames property of mdl.

plotpartialdependence(mdl,3,mdl.classnames);

according to this model, the probability of virginica increases with x3. the probability of setosa is about 0.33, from where x3 is 0 to around 2.5, and then the probability drops to almost 0.

train a gaussian process regression model using generated sample data where a response variable includes interactions between predictor variables. then, create ice plots that show the relationship between a feature and the predicted responses for each observation.

generate sample predictor data x1 and x2.

rng("default") % for reproducibility
n = 200;
x1 = rand(n,1)*2-1;
x2 = rand(n,1)*2-1;

generate response values that include interactions between x1 and x2.

y = x1-2*x1.*(x2>0) 0.1*rand(n,1);

create a gaussian process regression model using [x1 x2] and y.

mdl = fitrgp([x1 x2],y);

create a figure including a pdp (red line) for the first predictor x1, a scatter plot (circle markers) of x1 and predicted responses, and a set of ice plots (gray lines) by specifying conditional as "centered".

plotpartialdependence(mdl,1,"conditional","centered")

when conditional is "centered", plotpartialdependence offsets plots so that all plots start from zero, which is helpful in examining the cumulative effect of the selected feature.

a pdp finds averaged relationships, so it does not reveal hidden dependencies especially when responses include interactions between features. however, the ice plots clearly show two different dependencies of responses on x1.

train an ensemble of classification models and create two pdps, one using the training data set and the other using a new data set.

load the census1994 data set, which contains us yearly salary data, categorized as <=50k or >50k, and several demographic variables.

load census1994

extract a subset of variables to analyze from the tables adultdata and adulttest.

x = adultdata(:,["age","workclass","education_num","marital_status","race", ...
   "sex","capital_gain","capital_loss","hours_per_week","salary"]);
xnew = adulttest(:,["age","workclass","education_num","marital_status","race", ...
   "sex","capital_gain","capital_loss","hours_per_week","salary"]);

train an ensemble of classifiers with salary as the response and the remaining variables as predictors by using the function fitcensemble. for binary classification, fitcensemble aggregates 100 classification trees using the logitboost method.

mdl = fitcensemble(x,"salary");

inspect the class names in mdl.

mdl.classnames

ans = 2x1 categorical
     <=50k 
     >50k 

create a partial dependence plot of the scores predicted by mdl for the second class of salary (>50k) against the predictor age using the training data.

plotpartialdependence(mdl,"age",mdl.classnames(2))

create a pdp of the scores for class >50k against age using new predictor data from the table xnew.

plotpartialdependence(mdl,"age",mdl.classnames(2),xnew)

the two plots show similar shapes for the partial dependence of the predicted score of high salary (>50k) on age. both plots indicate that the predicted score of high salary rises fast until the age of 30, then stays almost flat until the age of 60, and then drops fast. however, the plot based on the new data produces slightly higher scores for ages over 65.

create a pdp to analyze relationships between predictors and anomaly scores for an isolationforest object. you cannot pass an isolationforest object directly to the plotpartialdependence function. instead, define a custom function that returns anomaly scores for the object, and then pass the function to plotpartialdependence.

load the 1994 census data stored in census1994.mat. the data set consists of demographic data from the us census bureau.

load census1994

census1994 contains the two data sets adultdata and adulttest.

train an isolation forest model for adulttest. the function iforest returns an isolationforest object.

rng("default") % for reproducibility
mdl = iforest(adulttest);

define the custom function myanomalyscores, which returns anomaly scores computed by the isanomaly function of isolationforest; the custom function definition appears at .

create a pdp of the anomaly scores against the variable age for the adulttest data set. plotpartialdependence accepts a custom model in the form of a function handle. the function represented by the function handle must accept predictor data and return a column vector or matrix with one row for each observation. specify the custom model as @(tbl)myanomalyscores(mdl,tbl) so that the custom function uses the trained model mdl and accepts predictor data.

plotpartialdependence(@(tbl)myanomalyscores(mdl,tbl),"age",adulttest)
xlabel("age")
ylabel("anomaly score")

function scores = myanomalyscores(mdl,tbl)
[~,scores] = isanomaly(mdl,tbl);
end

train a regression ensemble using the carsmall data set, and create a pdp plot and ice plots for each predictor variable using a new data set, carbig. then, compare the figures to analyze the importance of predictor variables. also, compare the results with the estimates of predictor importance returned by the function.

load the carsmall data set.

load carsmall

specify weight, cylinders, horsepower, and model_year as the predictor variables (x), and mpg as the response variable (y).

x = [weight,cylinders,horsepower,model_year];
y = mpg;

train a regression ensemble using x and y.

mdl = fitrensemble(x,y, ...
    "predictornames",["weight","cylinders","horsepower","model year"], ...
    "responsename","mpg");

create the importance of predictor variables by using the plotpartialdependence and functions. the plotpartialdependence function visualizes the relationships between a selected predictor and predicted responses. predictorimportance summarizes the importance of a predictor with a single value.

create a figure including a pdp plot (red line) and ice plots (gray lines) for each predictor by using plotpartialdependence and specifying "conditional","absolute". each figure also includes a scatter plot (circle markers) of the selected predictor and predicted responses. also, load the carbig data set and use it as new predictor data, xnew. when you provide xnew, the plotpartialdependence function uses xnew instead of the predictor data in mdl.

load carbig
xnew = [weight,cylinders,horsepower,model_year];
figure
t = tiledlayout(2,2,"tilespacing","compact");
title(t,"individual conditional expectation plots")
for i = 1 : 4
    nexttile
    plotpartialdependence(mdl,i,xnew,"conditional","absolute")
    title("")
end

compute estimates of predictor importance by using predictorimportance. this function sums changes in the mean squared error (mse) due to splits on every predictor, and then divides the sum by the number of branch nodes.

imp = predictorimportance(mdl);
figure
bar(imp)
title("predictor importance estimates")
ylabel("estimates")
xlabel("predictors")
ax = gca;
ax.xticklabel = mdl.predictornames;

the variable weight has the most impact on mpg according to predictor importance. the pdp of weight also shows that mpg has high partial dependence on weight. the variable cylinders has the least impact on mpg according to predictor importance. the pdp of cylinders also shows that mpg does not change much depending on cylinders.

train a generalized additive model (gam) with both linear and interaction terms for predictors. then, create a pdp with both linear and interaction terms and a pdp with only linear terms. specify whether to include interaction terms when creating the pdps.

load the ionosphere data set. this data set has 34 predictors and 351 binary responses for radar returns, either bad ('b') or good ('g').

load ionosphere

train a gam using the predictors x and class labels y. a recommended practice is to specify the class names. specify to include the 10 most important interaction terms.

mdl = fitcgam(x,y,"classnames",{'b','g'},"interactions",10);

mdl is a classificationgam model object.

list the interaction terms in mdl.

mdl.interactions

ans = 10×2
     1     5
     7     8
     6     7
     5     6
     5     7
     5     8
     3     5
     4     7
     1     7
     4     5

each row of interactions represents one interaction term and contains the column indexes of the predictor variables for the interaction term.

find the most frequent predictor in the interaction terms.

mode(mdl.interactions,"all")

ans = 5

the most frequent predictor in the interaction terms is the 5th predictor (x5). create pdps for the 5th predictor. to exclude interaction terms from the computation, specify "includeinteractions",false for the second pdp.

plotpartialdependence(mdl,5,mdl.classnames(1))
hold on
plotpartialdependence(mdl,5,mdl.classnames(1),"includeinteractions",false)
grid on
legend("linear and interaction terms","linear terms only")
title("pdps of posterior probabilities for 5th predictor")
hold off

the plot shows that the partial dependence of the scores (posterior probabilities) on x5 varies depending on whether the model includes the interaction terms, especially where x5 is between 0.2 and 0.45.

train a support vector machine (svm) regression model using the carsmall data set, and create a pdp for two predictor variables. then, extract partial dependence estimates from the output of plotpartialdependence. alternatively, you can get the partial dependence values by using the partialdependence function.

load the carsmall data set.

load carsmall

specify weight, cylinders, displacement, and horsepower as the predictor variables (tbl).

tbl = table(weight,cylinders,displacement,horsepower);

construct an svm regression model using tbl and the response variable mpg. use a gaussian kernel function with an automatic kernel scale.

mdl = fitrsvm(tbl,mpg,"responsename","mpg", ...
    "categoricalpredictors","cylinders","standardize",true, ...
    "kernelfunction","gaussian","kernelscale","auto");

create a pdp that visualizes partial dependence of predicted responses (mpg) on the predictor variables weight and cylinders. specify query points to compute the partial dependence for weight by using the querypoints name-value argument. you cannot specify the querypoints value for cylinders because it is a categorical variable. plotpartialdependence uses all categorical values.

pt = linspace(min(weight),max(weight),50)';
ax = plotpartialdependence(mdl,["weight","cylinders"],"querypoints",{pt,[]});
view(140,30) % modify the viewing angle

the pdp shows an interaction effect between weight and cylinders. the partial dependence of mpg on weight changes depending on the value of cylinders.

extract the estimated partial dependence of mpg on weight and cylinders. the xdata, ydata, and zdata values of ax.children are x-axis values (the first selected predictor values), y-axis values (the second selected predictor values), and z-axis values (the corresponding partial dependence values), respectively.

xval = ax.children.xdata;
yval = ax.children.ydata;
zval = ax.children.zdata;

alternatively, you can get the partial dependence values by using the partialdependence function.

[pd,x,y] = partialdependence(mdl,["weight","cylinders"],"querypoints",{pt,[]});

pd contains the partial dependence values for the query points x and y.

if you specify conditional as "absolute", plotpartialdependence creates a figure including a pdp, a scatter plot, and a set of ice plots. ax.children(1) and ax.children(2) correspond to the pdp and scatter plot, respectively. the remaining elements of ax.children correspond to the ice plots. the xdata and ydata values of ax.children(i) are x-axis values (the selected predictor values) and y-axis values (the corresponding partial dependence values), respectively.