Showing the correlation of two variables using a plot

Question

I have a list of $\{x_i,y_i\}$ pairs and I want to show that $x_i$'s and $y_i$'s have positive correlation. The result of Correlation[Transpose[data][[1]], Transpose[data][[2]]] is 0.856906 and ListPlot[data, PlotRange -> All, AxesLabel -> {X, Y}] results in the following plot:

I just wonder if there is any other (better) type of plot in Mathematica to show the correlation between the variables?

Answer 1

You can also calculate the Coefficient of Determination, R Squared.

This is the same as the correlation squared, but by making use of LinearModelFit you can create some additional graphics.

To make a sample distribution you can use this:

CreateDistribution[] := DynamicModule[{savepts = {{-1, -1}}},
  Dynamic[
   EventHandler[
    ListPlot[pts, AxesOrigin -> {0, 0}, 
     PlotRange -> {{0, 7}, {0, 5}}], 
    "MouseDown" :> (savepts = 
       pts = DeleteCases[
         Append[pts, MousePosition["Graphics"]], {-1, -1}])],
   Initialization :> (pts = savepts)]]

CreateDistribution[]

Just click to add some points. The data is collected in the variable pts.

Then calculate R Squared:

lm = LinearModelFit[Sort@pts, a, a]; r2 = lm["RSquared"];
Show[Plot[lm[x], {x, 0, 7}], ListPlot[pts], AxesOrigin -> {0, 0}, 
 PlotRange -> {{0, 7}, {0, 5}}, 
 PlotLabel -> 
  "The correlation is " <> 
   If[D[lm["BestFit"], a] < 0, "negative", "positive"], 
 Epilog -> 
  Inset[Style[
    "\!\(\*SuperscriptBox[\"R\", \"2\"]\) = " <> ToString[r2], 
    11], {1.5, 4.5}]]

Whether the correlation is positive or negative is obtained here from the derivative of the BestFit.

You can add standard deviation bands for sigma = 1, 2 & 3 like this.

lm = LinearModelFit[Sort@pts, {1, x (*, x^2 *)}, x];
{bands68[x_], bands95[x_], bands99[x_]} = 
  Table[lm["SinglePredictionBands", 
    ConfidenceLevel -> cl], {cl, {0.6827, 0.9545, 0.9973}}];
Show[ListPlot[Sort@pts], 
 Plot[{lm[x], bands68[x], bands95[x], bands99[x]}, {x, -0.15, 7.2}, 
  Filling -> {2 -> {1}, 3 -> {2}, 4 -> {3}}], AxesOrigin -> {0, 0}, 
 PlotRange -> {{0, 7}, {0, 5}}, ImageSize -> 480, Frame -> True]

Uncomment x^2 for a quadratic fit.

Added note

The linear band calculation steps are detailed here: How to derive SinglePredictionBands

Answer 2

Basic, but useful options

Something that certainly helps in representing your data is to set the range of x and y to be equal:

ListPlot[data, PlotRange -> {{0,3},{0,3}}, AxesLabel -> {X, Y}]

Moreover, if the slope is also important you should set the aspect ratio to 1:

ListPlot[data, PlotRange -> {{0,3},{0,3}}, AxesLabel -> {X, Y}, AspectRatio -> 1]

---

Confidence ellipses

In some cases it is convenient to add a confidence ellipse.

Let's generate some data:

data = Table[
        RandomVariate[BinormalDistribution[{50, 50}, {5, 10}, .8]], {1000}];

Here we get the estimates of the generated data:

estDist = EstimatedDistribution[data, 
            BinormalDistribution[
               {μ1, μ2}, {σ1, σ2}, ρ]]

And, here there are two useful functions that I found somewhere some time ago. Unfortunately, I'm not able to refer to the original author.

Ellipse[{{r_, M_}, m_}, {x_, y_}] := 
 Show[ContourPlot[({x, y} - r).M.({x, y} - r) == m, 
       {x, 0, 100}, {y, 0, 100}, ContourStyle -> {Red, Thick}], 
      ListPlot[{r}]
 ];

showEllipse[r0_, s1_, s2_, rho_, perc_] := 
  Module[{dchi2, fi, fic, V, InvV}, 
    Clear[σ, ρ];
    dchi2 = Quantile[ChiSquareDistribution[2], perc]; 
    fi = Which[s1 == s2, π/4 Sign[rho], 
               s1 > s2, 1/2 ArcTan[(2 rho s1 s2)/(s1^2 - s2^2)], 
               s1 < s2, 1/2 ArcTan[(2 rho s1 s2)/(s1^2 - s2^2)] - π/2]; 
    fic = If[fi < 0, π + fi, fi]; 
    V = {{s1^2, rho s1 s2}, {rho s1 s2, s2^2}}; 
    InvV = Inverse[V]; 
    Ellipse[{{r0, InvV}, dchi2}, {x, y}]
  ]

Now we can plot the data and the 95% confidence ellipse:

Show[{
 ListPlot[data, PlotRange -> {{0, 100}, {0, 100}}, AspectRatio ->1],
 showEllipse[estDist[[1]], estDist[[2, 1]], estDist[[2, 2]], estDist[[3]], .95]
}]

Or, even multiple confidence ellipses (95% and 99%):

Show[{
 ListPlot[data, PlotRange -> {{0, 100}, {0, 100}}, AspectRatio ->1],
 showEllipse[estDist[[1]], estDist[[2, 1]], estDist[[2, 2]], estDist[[3]], .99],
 showEllipse[estDist[[1]], estDist[[2, 1]], estDist[[2, 2]], estDist[[3]], .95]
}]

Answer 3

Another approach could be also based on the smooth kernel distribution.

If this are our data:

data = Table[RandomVariate[BinormalDistribution[{50, 50}, {5, 7}, .8]], {1000}];

We can obtain a representation of a multivariate smooth kernel distribution based on the data values:

skd = SmoothKernelDistribution[data];

And here we get the plot:

ContourPlot[
    Evaluate@PDF[skd, {x, y}], {x, 30, 70}, {y, 30, 70},
    PlotRange -> All, PlotPoints -> 100]

This approach can be helpful to get a glimpse on how the data are distributed. If you consider a case in which the data originate from two different distributions, plotting them in this way gives you a clear picture.

data1 = Table[RandomVariate[BinormalDistribution[{50, 50}, {5, 7}, .8]], {500}];
data2 = Table[RandomVariate[BinormalDistribution[{45, 55}, {7, 5}, .8]], {500}];
data = Partition[Flatten[Join[{data1, data2}]], 2];
skd2 = SmoothKernelDistribution[data];
ContourPlot[Evaluate@PDF[skd2, {x, y}], {x, 30, 70}, {y, 30, 70},
    PlotRange -> All, PlotPoints -> 100]

Answer 4

I tried to use the above mentioned codes for plotting error ellipses for 2D-Data. However, I did not get the anticipated results, because an error occured when Mathematica tried to solve the equality in the function Counterplot for my data.

I found another solution based on the explicit calculation of the ellipse by means of covariance analysis. The ellipse is then rotated to coincide with the principal axis of the correlated data. The original explanation can be found at:

http://www.visiondummy.com/2014/04/draw-error-ellipse-representing-covariance-matrix/

Since I did not find an equivalent solution at stackexchange, I adapted the provided matlab code given by Vincent Spruyt to mathematica and answered to this 'historic' post.

define some random correlated data:

(*Random correlated data generation*)
s = 3;
var = 0.3;
mus= Table[RandomReal[], {i, 10000}];

x = RandomVariate[NormalDistribution[#, var]] & /@ (+s mus);
y = RandomVariate[NormalDistribution[#, var]] & /@ (-s mus);
data = {x, y}\[Transpose];

Perform a Covariance analysis and create a 2D-list, that contains the ellipse as points (percLevel is the percentage covering a given confidence interval)

(*define error Ellipse*)
ErrorEllipse[data_, percLevel_] :=
  Module[
   {eVa, eVec, coors, dchi, thetaGrid, ellipse, rEllipse},
   {eVa, eVec} = Eigensystem@Covariance[data];

   (* Get the coordinates of the data mean*)
   coors = Mean[data];

   (*Get the perc [%] confidence interval error ellipse*)
   dchi = √Quantile[ChiSquareDistribution[2], percLevel/100];

   (* define error ellipse in x and y coordinates*)
  thetaGrid = Table[i, {i, 0, 2 π, 2 π/99}];
   ellipse = {dchi √eVa[[1]] Cos[thetaGrid], 
     dchi √eVa[[2]] Sin[thetaGrid]};

   (* rotate the ellipse and center ellipse at coors*)
   rEllipse = coors + # & /@ (ellipse\[Transpose].eVec)
   ];

Use LineListPlot to show the error ellipses:

(* visualize results*)
percentOneSigma = 66.3;
percentTwoSigma = 95.4;

Show[
 ListPlot[data],
 ListLinePlot[ErrorEllipse[data, percentOneSigma],PlotStyle -> {Thick, Red}],
 ListLinePlot[ErrorEllipse[data, percentTwoSigma],PlotStyle -> {Thick, Blue}],
 ListPlot[{Mean[data]}, PlotStyle -> {Thick, Red}]
 ]

Hope, this might be useful.

Cheers

Answer 5

Needs["MultivariateStatistics`"]

data = Table[RandomVariate[BinormalDistribution[{50, 50}, {5, 7}, .8]], {1000}];

Show[
 ListPlot[data,
  Epilog -> {
    Text[Framed["1.00", Background -> LightBlue], {45, 25}],
    Text[Framed["0.98", Background -> LightBlue], {50, 33}],
    Text[Framed["0.90", Background -> LightBlue], {55, 43}]
    },
  PlotMarkers -> {Automatic, Tiny},
  PlotTheme -> "Detailed"],
 Graphics[EllipsoidQuantile[data, {0.9, 0.98, 1.0}]],
 ImageSize -> 500,
 PlotRange -> All]

Update

data = Table[
   RandomVariate[BinormalDistribution[{5, 5}, {5, 5}, .5]], {500}];

c = {Purple, Red, Green, Blue, Black};
q = {0.5, 0.7, 0.9, 0.98, 1.0};

Show[
 ListPlot[data,
  PlotLegends -> PointLegend[c, q],
  PlotMarkers -> {Automatic, Tiny},
  PlotTheme -> "Detailed"],
 Graphics[Transpose[{c, EllipsoidQuantile[data, q]}]],
 ImageSize -> 500,
 PlotRange -> All]