<p><div><div><div><div data-width="100" data-show-faces="false" data-send="true" data-layout="button_count" data-height="21" data-href="http://www.r-bloggers.com/how-to-use-r-for-matching-samples-propensity-score/"></div></div><div><a data-url="http://www.r-bloggers.com/how-to-use-r-for-matching-samples-propensity-score/" data-counturl="http://www.r-bloggers.com/how-to-use-r-for-matching-samples-propensity-score/" data-count="horizontal" data-text="How to use R for matching samples (propensity score)" class="twitter-share-button" data-via="rbloggers" href="https://twitter.com/share"></a></div><div></div></div><div></div></div>
<div>
(This article was first published on <strong><a href="http://datascienceplus.com/how-to-use-r-for-matching-samples-propensity-score/"> DataScience+</a></strong>, and kindly contributed to <a rel="nofollow" href="http://www.r-bloggers.com/">R-bloggers)</a> </div>
<p>According to <a target="_blank" rel="nofollow" href="https://en.wikipedia.org/wiki/Propensity_score_matching">Wikipedia</a>, propensity score matching (PSM) is a “statistical matching technique that attempts to estimate the effect of a treatment, policy, or other intervention by accounting for the covariates that predict receiving the treatment”. In a broader sense, propensity score analysis assumes that an unbiased comparison between samples can only be made when the subjects of both samples have similar characteristics. Thus, PSM can not only be used as “an alternative method to estimate the effect of receiving treatment when random assignment of treatments to subjects is not feasible” (<a target="_blank" rel="nofollow" href="http://umanitoba.ca/faculties/health_sciences/medicine/units/community_health_sciences/departmental_units/mchp/protocol/media/propensity_score_matching.pdf">Thavaneswaran 2008</a>). It can also be used for the comparison of samples in epidemiological studies. Let’s give an example:</p>
<p>Health-related quality of life (HRQOL) is considered an important outcome in cancer therapy. One of the most frequently used instruments to measure HRQOL in cancer patients is the core quality-of-life questionnaire of the <em>European Organisation for Research and Treatment of Cancer</em>. The <a target="_blank" rel="nofollow" href="http://groups.eortc.be/qol/sites/default/files/img/slider/specimen_qlq-c30_english.pdf">EORTC QLQ-C30</a> is a 30-item instrument comprised of five functioning scales, nine symptom scales and one scale measuring Global quality of life. All scales have a score range between 0 and 100. While high scores of the symptom scales indicate a high burden of symptoms, high scores of the functioning scales and on the GQoL scale indicate better functioning resp. quality of life. </p>
<p>However, without having any reference point, it is difficult if not impossible to interpret the scores. Fortunately, the EORTC QLQ-C30 questionnaire was used in several general population surveys. Therefore, patient scores may be compared against scores of the general population. This makes it far easier to decide whether the burden of symptoms or functional impairments can be attributed to cancer (treatment) or not. PSM can be used to make both patient and population samples comparable by matching for relevant demographic characteristics like age and sex. </p>
<p>In this blog post, I show how to do PSM using R. A more comprehensive PSM guide can be found under: “<a target="_blank" rel="nofollow" href="http://pareonline.net/getvn.asp?v=19&n=18">A Step-by-Step Guide to Propensity Score Matching in R</a>“.</p>
<h2>Creating two random dataframes</h2>
<p>Since we don’t want to use real-world data in this blog post, we need to emulate the data. This can be easily done using the <code>Wakefield</code> package. </p>
<p>In a first step, we create a dataframe named <code>df.patients</code>. We want the dataframe to contain specifications of age and sex for 250 patients. The patients’ age shall be between 30 and 78 years. Furthermore, 70% of patients shall be male.</p>
<pre>
set.seed(1234)
df.patients &lt;- r_data_frame(n = 250, 
                            age(x = 30:78, 
                                name = 'Age'), 
                            sex(x = c(&quot;Male&quot;, &quot;Female&quot;), 
                                prob = c(0.70, 0.30), 
                                name = &quot;Sex&quot;))
df.patients$Sample &lt;- as.factor('Patients')
</pre>
<p>The <code>summary</code>-function returns some basic information about the dataframe created. As we can see, the mean age of the patient sample is 53.7 and roughly 70% of the patients are male (69.2%).</p>
<pre>summary(df.patients)<em>
## Age Sex Sample ## Min. :30.00 Male :173 Patients:250 ## 1st Qu.:42.00 Female: 77 ## Median :54.00 ## Mean :53.71 ## 3rd Qu.:66.00 ## Max. :78.00
</em></pre>
<p>In a second step, we create another dataframe named <code>df.population</code>. We want this dataframe to comprise the same variables as <code>df.patients</code> with different specifications. With 18 to 80, the age-range of the population shall be wider than in the patient sample and the proportion of female and male patients shall be the same.</p>
<pre>
set.seed(1234)
df.population &lt;- r_data_frame(n = 1000, 
                              age(x = 18:80, 
                                  name = 'Age'), 
                              sex(x = c(&quot;Male&quot;, &quot;Female&quot;), 
                                  prob = c(0.50, 0.50), 
                                  name = &quot;Sex&quot;))
df.population$Sample &lt;- as.factor('Population')
</pre>
<p>The following table shows the sample’s mean age (49.5 years) and the proportion of men (48.5%) and women (51.5%).</p>
<pre>summary(df.population)<em>
## Age Sex Sample ## Min. :18.00 Male :485 Population:1000 ## 1st Qu.:34.00 Female:515 ## Median :50.00 ## Mean :49.46 ## 3rd Qu.:65.00 ## Max. :80.00
</em></pre>
<h2>Merging the dataframes</h2>
<p>Before we match the samples, we need to merge both dataframes. Based on the variable <em>Sample</em>, we create a new variable named <em>Group</em> (type <em>logic</em>) and a further variable (<em>Distress</em>) containing information about the individuals’ level of distress. The <em>Distress</em> variable is created using the <code>age</code>-function of the <code>Wakefield</code> package. As we can see, women will have higher levels of distress.</p>
<pre>
mydata &lt;- rbind(df.patients, df.population)
mydata$Group &lt;- as.logical(mydata$Sample == 'Patients')
mydata$Distress &lt;- ifelse(mydata$Sex == 'Male', age(nrow(mydata), x = 0:42, name = 'Distress'),
                                                age(nrow(mydata), x = 15:42, name = 'Distress'))
</pre>
<p>When we compare the distribution of age and sex in both samples, we discover significant differences:</p>
<pre>
pacman::p_load(tableone)
table1 &lt;- CreateTableOne(vars = c('Age', 'Sex', 'Distress'), 
                         data = mydata, 
                         factorVars = 'Sex', 
                         strata = 'Sample')
table1 &lt;- print(table1, 
                printToggle = FALSE, 
                noSpaces = TRUE)
kable(table1[,1:3],  
      align = 'c', 
      caption = 'Table 1: Comparison of unmatched samples')
</pre>
<p>Table 1: Comparison of unmatched samples</p>
<table>
<thead>
<tr>
<th align="left"></th>
<th align="center">Patients</th>
<th align="center">Population</th>
<th align="center">p</th>
</tr>
</thead>
<tbody>
<tr>
<td align="left">n</td>
<td align="center">250</td>
<td align="center">1000</td>
<td align="center"></td>
</tr>
<tr>
<td align="left">Age (mean (sd))</td>
<td align="center">53.71 (13.88)</td>
<td align="center">49.46 (18.33)</td>
<td align="center">0.001</td>
</tr>
<tr>
<td align="left">Sex = Female (%)</td>
<td align="center">77 (30.8)</td>
<td align="center">515 (51.5)</td>
<td align="center">&lt;0.001</td>
</tr>
<tr>
<td align="left">Distress (mean (sd))</td>
<td align="center">22.86 (11.38)</td>
<td align="center">25.13 (11.11)</td>
<td align="center">0.004</td>
</tr>
</tbody>
</table>
<p>Furthermore, the level of distress seems to be significantly higher in the population sample.</p>
<h2>Matching the samples</h2>
<p>Now, that we have completed preparation and inspection of data, we are going to match the two samples using the <code>matchit</code>-function of the <code>MatchIt</code> package. The method command <code>method=&quot;nearest&quot;</code> specifies that the nearest neighbors method will be used. Other matching methods are <em>exact matching, subclassification, optimal matching, genetic matching</em>, and <em>full matching</em> (<code>method = c(&quot;exact&quot;, &quot;subclass&quot;, &quot;optimal&quot;, &quot;&quot;genetic&quot;, &quot;full&quot;)</code>). The ratio command <code>ratio = 1</code> indicates a one-to-one matching approach. With regard to our example, for each case in the patient sample exactly one case in the population sample will be matched. Please also note that the Group variable needs to be <em>logic</em> (TRUE vs. FALSE).</p>
<pre>
set.seed(1234)
match.it &lt;- matchit(Group ~ Age + Sex, data = mydata, method=&quot;nearest&quot;, ratio=1)
a &lt;- summary(match.it)
</pre>
<p>For further data presentation, we save the output of the <code>summary</code>-function into a variable named a.</p>
<p>After matching the samples, the size of the population sample was reduced to the size of the patient sample (n=250; see table 2).</p>
<pre>
kable(a$nn, digits = 2, align = 'c', 
      caption = 'Table 2: Sample sizes')
</pre>
<p>Table 2: Sample sizes</p>
<table>
<thead>
<tr>
<th align="left"></th>
<th align="center">Control</th>
<th align="center">Treated</th>
</tr>
</thead>
<tbody>
<tr>
<td align="left">All</td>
<td align="center">1000</td>
<td align="center">250</td>
</tr>
<tr>
<td align="left">Matched</td>
<td align="center">250</td>
<td align="center">250</td>
</tr>
<tr>
<td align="left">Unmatched</td>
<td align="center">750</td>
<td align="center">0</td>
</tr>
<tr>
<td align="left">Discarded</td>
<td align="center">0</td>
<td align="center">0</td>
</tr>
</tbody>
</table>
<p>The following output shows, that the distributions of the variables Age and Sex are nearly identical after matching.</p>
<pre>
kable(a$sum.matched[c(1,2,4)], digits = 2, align = 'c', 
      caption = 'Table 3: Summary of balance for matched data')
</pre>
<p>Table 3: Summary of balance for matched data</p>
<table>
<thead>
<tr>
<th align="left"></th>
<th align="center">Means Treated</th>
<th align="center">Means Control</th>
<th align="center">Mean Diff</th>
</tr>
</thead>
<tbody>
<tr>
<td align="left">distance</td>
<td align="center">0.23</td>
<td align="center">0.23</td>
<td align="center">0.00</td>
</tr>
<tr>
<td align="left">Age</td>
<td align="center">53.71</td>
<td align="center">53.65</td>
<td align="center">0.06</td>
</tr>
<tr>
<td align="left">SexMale</td>
<td align="center">0.69</td>
<td align="center">0.69</td>
<td align="center">0.00</td>
</tr>
<tr>
<td align="left">SexFemale</td>
<td align="center">0.31</td>
<td align="center">0.31</td>
<td align="center">0.00</td>
</tr>
</tbody>
</table>
<p>The distributions of propensity scores can be visualized using the <code>plot</code>-function which is part of the <code>MatchIt</code> package .</p>
<pre>
plot(match.it, type = 'jitter', interactive = FALSE)
</pre>
<p>Here is the plot:<br>
<a title="" data-rel="lightbox-0" target="_blank" href="http://i0.wp.com/datascienceplus.com/wp-content/uploads/2016/06/plot-19.png"><img sizes="(max-width: 490px) 100vw, 490px" srcset="http://i0.wp.com/datascienceplus.com/wp-content/uploads/2016/06/plot-19-490x327.png?resize=456%2C304 490w, http://datascienceplus.com/wp-content/uploads/2016/06/plot-19-768x512.png 768w, http://datascienceplus.com/wp-content/uploads/2016/06/plot-19-600x400.png 600w, http://datascienceplus.com/wp-content/uploads/2016/06/plot-19-300x200.png 300w, http://datascienceplus.com/wp-content/uploads/2016/06/plot-19.png 864w" alt="plot-19" class="wp-image-4188" src="http://i0.wp.com/datascienceplus.com/wp-content/uploads/2016/06/plot-19-490x327.png?resize=456%2C304" data-recalc-dims="1"></a></p>
<h2>Saving the matched samples</h2>
<p>Finally, the matched samples will be saved into a new dataframe named <code>df.match</code>.</p>
<pre>
df.match &lt;- match.data(match.it)[1:ncol(mydata)]
rm(df.patients, df.population)
</pre>
<p>Eventually, we can check whether the differences in the level of distress between both samples are still significant.</p>
<pre>
pacman::p_load(tableone)
table4 &lt;- CreateTableOne(vars = c('Age', 'Sex', 'Distress'), 
                         data = df.match, 
                         factorVars = 'Sex', 
                         strata = 'Sample')
table4 &lt;- print(table4, 
                printToggle = FALSE, 
                noSpaces = TRUE)
kable(table4[,1:3],  
      align = 'c', 
      caption = 'Table 4: Comparison of matched samples')
</pre>
<table>
<thead>
<tr>
<th align="left"></th>
<th align="center">Patients</th>
<th align="center">Population</th>
<th align="center">p</th>
</tr>
</thead>
<tbody>
<tr>
<td align="left">n</td>
<td align="center">250</td>
<td align="center">250</td>
<td align="center"></td>
</tr>
<tr>
<td align="left">Age (mean (sd))</td>
<td align="center">53.71 (13.88)</td>
<td align="center">53.65 (13.86)</td>
<td align="center">0.961</td>
</tr>
<tr>
<td align="left">Sex = Female (%)</td>
<td align="center">77 (30.8)</td>
<td align="center">77 (30.8)</td>
<td align="center">1.000</td>
</tr>
<tr>
<td align="left">Distress (mean (sd))</td>
<td align="center">22.86 (11.38)</td>
<td align="center">24.13 (11.88)</td>
<td align="center">0.222</td>
</tr>
</tbody>
</table>
<p>With a p-value of 0.222, Student’s t-test does not indicate significant differences anymore. Thus, PSM helped to avoid an alpha mistake.</p>
<p>PS 1: The packages used in this blog post can be loaded/installed using the following code:</p>
<pre>
pacman::p_load(knitr, wakefield, MatchIt, tableone, captioner)
</pre>
<p>PS 2: Thanks very much to my colleague <a target="_blank" rel="nofollow" href="http://medpsy.uniklinikum-leipzig.de/medpsych.site,postext,mitarbeiter,a_id,5503.html">Katharina Kuba</a> for telling me about the <code>MatchIt</code> package.</p>
<ul></ul>
<p><strong>
<p>Related Post</p>
<p>
<ol>
<li><a title="Permanent Link toR for Publication by Page Piccinini: Lesson 1 – R Basics" target="_blank" rel="bookmark" href="http://datascienceplus.com/r-for-publication-by-page-piccinini-lesson-1-r-basics/">R for Publication by Page Piccinini: Lesson 1 – R Basics</a></li>
<li><a title="Permanent Link toR for Publication by Page Piccinini: Lesson 0 – Introduction and Set-up" target="_blank" rel="bookmark" href="http://datascienceplus.com/r-for-publication-by-page-piccinini-lesson-0-introduction-and-set-up/">R for Publication by Page Piccinini: Lesson 0 – Introduction and Set-up</a></li>
<li><a title="Permanent Link toHow to export Regression results from R to MS Word" target="_blank" rel="bookmark" href="http://datascienceplus.com/how-to-export-regression-results-from-r-to-ms-word/">How to export Regression results from R to MS Word</a></li>
<li><a title="Permanent Link toLearn R by Intensive Practice" target="_blank" rel="bookmark" href="http://datascienceplus.com/learn-r-by-intensive-practice/">Learn R by Intensive Practice</a></li>
<li><a title="Permanent Link toLearn R from the Ground Up" target="_blank" rel="bookmark" href="http://datascienceplus.com/learn-r-from-the-ground-up/">Learn R from the Ground Up</a></li>
</ol></strong>
<div>
<a target="_blank" href="http://feeds.feedburner.com/~ff/datascienceplus?a=ykPpvmhuhVM:0axfn0gGh0w:D7DqB2pKExk"><img src="http://feeds.feedburner.com/~ff/datascienceplus?i=ykPpvmhuhVM:0axfn0gGh0w:D7DqB2pKExk" border="0"></a> <a target="_blank" href="http://feeds.feedburner.com/~ff/datascienceplus?a=ykPpvmhuhVM:0axfn0gGh0w:qj6IDK7rITs"><img src="http://feeds.feedburner.com/~ff/datascienceplus?d=qj6IDK7rITs" border="0"></a> <a target="_blank" href="http://feeds.feedburner.com/~ff/datascienceplus?a=ykPpvmhuhVM:0axfn0gGh0w:yIl2AUoC8zA"><img src="http://feeds.feedburner.com/~ff/datascienceplus?d=yIl2AUoC8zA" border="0"></a>
</div>
<p><div><div><div><div data-width="55" data-show-faces="false" data-send="true" data-layout="box_count" data-height="62" data-href="http://www.r-bloggers.com/how-to-use-r-for-matching-samples-propensity-score/"></div></div><div><a data-url="http://www.r-bloggers.com/how-to-use-r-for-matching-samples-propensity-score/" data-counturl="http://www.r-bloggers.com/how-to-use-r-for-matching-samples-propensity-score/" data-count="vertical" data-text="How to use R for matching samples (propensity score)" class="twitter-share-button" data-via="rbloggers" href="https://twitter.com/share"></a></div><div></div></div><div></div></div>
<div>
<div>To <strong>leave a comment</strong> for the author, please follow the link and comment on their blog: <strong><a href="http://datascienceplus.com/how-to-use-r-for-matching-samples-propensity-score/"> DataScience+</a></strong>.</div>
<hr>
<a rel="nofollow" href="http://www.r-bloggers.com/">R-bloggers.com</a> offers <strong><a rel="nofollow" href="http://feedburner.google.com/fb/a/mailverify?uri=RBloggers">daily e-mail updates</a></strong> about <a title="The R Project for Statistical Computing" rel="nofollow" href="http://www.r-project.org/">R</a> news and <a title="R tutorials" rel="nofollow" href="http://www.r-bloggers.com/?s=tutorial">tutorials</a> on topics such as: <a title="Data science" rel="nofollow" href="http://www.r-bloggers.com/?s=data%20science">Data science</a>, <a title="Big Data" rel="nofollow" href="http://www.r-bloggers.com/?s=Big%20Data">Big Data, <a title="R jobs" rel="nofollow" href="http://www.r-users.com/">R jobs</a>, visualization (<a title="ggplot and ggplot2 tutorials" rel="nofollow" href="http://www.r-bloggers.com/?s=ggplot2">ggplot2</a>, <a title="Boxplots using lattice and ggplot2 tutorials" rel="nofollow" href="http://www.r-bloggers.com/?s=boxplot">Boxplots</a>, <a title="Maps and gis" rel="nofollow" href="http://www.r-bloggers.com/?s=map">maps</a>, <a title="Animation in R" rel="nofollow" href="http://www.r-bloggers.com/?s=animation">animation</a>), programming (<a title="RStudio IDE for R" rel="nofollow" href="http://www.r-bloggers.com/?s=RStudio">RStudio</a>, <a title="Sweave and literate programming" rel="nofollow" href="http://www.r-bloggers.com/?s=sweave">Sweave</a>, <a title="LaTeX in R" rel="nofollow" href="http://www.r-bloggers.com/?s=LaTeX">LaTeX</a>, <a title="SQL and databases" rel="nofollow" href="http://www.r-bloggers.com/?s=SQL">SQL</a>, <a title="Eclipse IDE for R" rel="nofollow" href="http://www.r-bloggers.com/?s=eclipse">Eclipse</a>, <a title="git and github, Version Control System" rel="nofollow" href="http://www.r-bloggers.com/?s=git">git</a>, <a title="Large data in R using Hadoop" rel="nofollow" href="http://www.r-bloggers.com/?s=hadoop">hadoop</a>, <a title="Web Scraping of google, facebook, yahoo, twitter and more using R" rel="nofollow" href="http://www.r-bloggers.com/?s=Web+Scraping">Web Scraping</a>) statistics (<a title="Regressions and ANOVA analysis tutorials" rel="nofollow" href="http://www.r-bloggers.com/?s=regression">regression</a>, <a title="principal component analysis tutorial" rel="nofollow" href="http://www.r-bloggers.com/?s=PCA">PCA</a>, <a title="Time series" rel="nofollow" href="http://www.r-bloggers.com/?s=time+series">time series</a>, <a title="finance trading" rel="nofollow" href="http://www.r-bloggers.com/?s=trading">trading</a>) and more...
<img height="1" alt="" width="1" src="http://feeds.feedburner.com/~r/RBloggers/~4/3No3Vq-cMdM"></a></div>