Visualize pre-post comparison of intervention #rstats

My sjPlot-package was just updated on CRAN, introducing a new function called sjp.emm.int to plot estimated marginal means (least-squares means) of linear models with interaction terms. Or: plotting adjusted means of an ANCOVA.

The idea to this function came up when we wanted to analyze the effect of an intervention (an educational programme on knowledge about mental disorders and associated stigma) between two groups: a “treatmeant” group (city) where a campaign on mental disorders was conducted and another city without this campaign. People from both cities were asked about their attitudes and knowledge about specific mental disorders at t0 before the campaign started in the one city. Some month later (t1), again people from both cities were asked the same questions. The intention was to see a) whether there were differences in knowledge and pro-social attidutes of people towards mental disorders and b) if the compaign successfully reduces stigma and increases knowledge.

To analyse these questions, we used an ANCOVA with knowledge and stigma score as dependent variables, “city” and “time” (t0 versus t1) as predictors and adjusted for covariates like age, sex, education etc. The estimated marginal means (or least-squares means) show you the differences of the dependent variable.

Here’s an example plot, quickly done with the sjp.emm.int function:
sjpemmint

Since the data is not publicly available, I’ve set an an documentation with reproducable examples (though those example do not fit very well…).

The latest development snapshot of my package is available on GitHub.

BTW: You may have noticed that this function is quite similar to the sjp.lm.int function for visually interpreting interaction terms in linear models…

Simply creating various scatter plots with ggplot #rstats

Inspired by these two postings, I thought about including a function in my package for simply creating scatter plots.

In my package, there’s a function called sjp.scatter for creating scatter plots. To reproduce these examples, first load the package and then attach the sample data set:

data(efc)

The simplest function call is by just providing two variables, one for the x- and one for the y-axis:

sjp.scatter(efc$c160age, efc$e17age)

which plots following graph:
sct_01

If you have continuous variables with a larger scale, you shouldn’t have problems with overplotting or overlaying dots. However, this problem usually occurs, if you have variables with just a few categories (factor levels). The function automatically estimates the amount of overlaying dots and then automatically jitters them, like in following example, which also includes a marginal rug-plot:

sjp.scatter(efc$e16sex,efc$neg_c_7, efc$c172code, showRug=TRUE)

sct_02

The same plot, when auto-jittering is turned off, would look like this:

sjp.scatter(efc$e16sex,efc$neg_c_7, efc$c172code,
            showRug=TRUE, autojitter=FALSE)

sct_03

You can also add a grouping variable. The scatter plot is then “divided” into as many groups as indicated by the grouping variable. In the next example, two variables (elder’s and carer’s age) are grouped by different dependency levels of the elderly. Additionally, a fitted line for each group is plotted:

sjp.scatter(efc$c160age,efc$e17age, efc$e42dep, title="Scatter Plot",
            legendTitle=sji.getVariableLabels(efc)['e42dep'],
            legendLabels=sji.getValueLabels(efc)[['e42dep']],
            axisTitle.x=sji.getVariableLabels(efc)['c160age'],
            axisTitle.y=sji.getVariableLabels(efc)['e17age'],
            showGroupFitLine=TRUE)

sct_04

If the groups are difficult to distinguish in a single plot area, the graph can be faceted by groups. This is shown in the last example, where the same scatter plot as above is plotted with facets for each group:

sjp.scatter(efc$c160age,efc$e17age, efc$e42dep, title="Scatter Plot",
            legendTitle=sji.getVariableLabels(efc)['e42dep'],
            legendLabels=sji.getValueLabels(efc)[['e42dep']],
            axisTitle.x=sji.getVariableLabels(efc)['c160age'],
            axisTitle.y=sji.getVariableLabels(efc)['e17age'],
            showGroupFitLine=TRUE, useFacetGrid=TRUE, showSE=TRUE)

sct_05

Find a complete overview of the various function options in the package-help or at inside-r.

Comparing multiple (g)lm in one graph #rstats

It’s been a while since a user of my plotting-functions asked whether it would be possible to compare multiple (generalized) linear models in one graph (see comment). While it is already possible to compare multiple models as table output, I now managed to build a function that plots several (g)lm-objects in a single ggplot-graph.

The following examples are take from my sjPlot package which is available on CRAN. Once you’ve installed the package, you can run one of the examples provided in the function’s documentation:

# prepare dummy variables for binary logistic regression
y1 <- ifelse(swiss$Fertility<median(swiss$Fertility), 0, 1)
y2 <- ifelse(swiss$Infant.Mortality<median(swiss$Infant.Mortality), 0, 1)
y3 <- ifelse(swiss$Agriculture<median(swiss$Agriculture), 0, 1)

# Now fit the models. Note that all models share the same predictors
# and only differ in their dependent variable (y1, y2 and y3)
fitOR1 <- glm(y1 ~ swiss$Education+swiss$Examination+swiss$Catholic,
              family=binomial(link="logit"))
fitOR2 <- glm(y2 ~ swiss$Education+swiss$Examination+swiss$Catholic,
              family=binomial(link="logit"))
fitOR3 <- glm(y3 ~ swiss$Education+swiss$Examination+swiss$Catholic,
              family=binomial(link="logit"))

# plot multiple models
sjp.glmm(fitOR1, fitOR2, fitOR3)

multiodds1

Thanks to the help of a stackoverflow user, I now know that the order of aes-parameters matters in case you have dodged positioning of geoms on a discrete scale. An example: I use following code in my function ggplot(finalodds, aes(y=OR, x=xpos, colour=grp, alpha=pa)) to apply different colours to each model and setting an alpha-level for geoms depending on the p-level. If the alpha-aes would appear before the colour-aes, the order of lines representing a model may be different for different x-values (see stackoverflow example).

Another more appealing example (not reproducable, since it relies on data from a current research project):
multiodds2

And finally an example where p-levels are represented by different shapes and non-significant odds have a lower alpha-level:
multiodds3

sjPlot 0.9 (data visualization package) now on CRAN #rstats

Since version 0.8, my package for data visualization using ggplot has been released on the Comprehensive R Archive Network (CRAN), which means you can simply install the package with install.packages("sjPlot").

Last week, version 0.9 was released. Binaries are already available for OS X and Windows, and source code for Linux. Further updates will no longer be announced on this blog (except for new functions which may be described in dedicated blog postings), so please use the update function in order make sure you are using the latest package version.

sjPlot – data visualization for statistics (in social science) #rstats

I’d like to announce the release of version 0.7 of my R package for data visualization and give a small overview of this package (download and installation instructions can be found on the package page, detailed examples on RPubs).

What does this package do?
In short, the functions in this package mostly do two things:

  1. compute basic or advanced statistical analyses
  2. either plot the results as ggplot-diagram or print them as html-table

However, meanwhile the amount of functions has increased, hence you’ll also find some utility functions beside the plotting functions.

How does this package help me?
Basically, this package either helps those users…

  • who have difficulties using and/or understanding all possibilities that ggplot offers to create plots, simply by providing intuitive function parameters, which allow for manipulating the appearance of plots; or
  • who don’t want to set up complex ggplot-object each time from the scratch; or
  • like quick inspections of (basic) statistics via (html-)tables that are shown in the GUI viewer pane or default browser; or
  • want easily create beautiful table outputs that can be imported in office applications.

Furthermore, for advanced users, each functions returns either the prepared ggplot-object (in case of sjp-plotting functions) or the HTML-tables (in case of sjt-table-output functions), which than can be manipulated even further (for instance, for ggplot-objects, you can specify certain parameters that cannot be modified via the sjPlot package or html-tables could be integrated into knitr-documents).

What are all these functions about?
There’s a certain naming convention for the functions:

  • sjc – collection of functions useful for carrying out cluster analyses
  • sji – collection of functions for data import and manipulation
  • sjp – collection plotting functions, the “core” of this package
  • sjt – collection of function that create (HTML) table outputs (instead of ggplot-graphics
  • sju – collection of statistical utility functions

Use cases?

  • You can plot results of Anova, correlations, histograms, box plots, bar plots, (generalized) linear models, likert scales, PCA, proportional tables as bar chart etc.
  • You can create plots to analyse model assumptions (lm, glm), predictor interactions, multiple contigency tables etc.
  • You can create table outputs instead of graphs for most plotting functions
  • With the import and utility functions, you can, for instance, extract beta coefficients of linear models, convert numeric scales into grouped factors, perform statistical tests, import SPSS data sets (and retrieve variable and value labels from the importet data), convert factors to numeric variables (and vice versa)…

Final remarks
At the bottom of my package page you’ll find some examples of selected functions that have been published on this blog before I created the package. Furthermore, the package includes a sample dataset from one of my research projects. Once the package is installed, you can test each function by running the examples. All news and recent changes can be found in the NEWS section of the package help (type ?sjPlot to access the help file inside R).

I tried to write a very comprehensive documentation for each function and their parameters, hopefully this will help you using my package…

Any comments, suggestions etc. are very welcome!

sjPlotting functions now as package available #rstats

This weekend I had some time to deal with package building in R. After some struggling, I now managed to setup RStudio, Roxygen and MikTex properly so I can compile my collection of R-scripts into a package that even succeeds the package check.

Downloads (package and manual) as well as package description are available at the package information page!

Since the packages successfully passed the package check and a manual could also be created, I’ll probably submit my package to the CRAN. Currently, I’m only able to compile the source and the Windows binaries of the package, because at home I use RStudio on my Mac with OS X 10.9 Mavericks. It seems that there’s an issue with the GNU Tar on Mavericks, which is needed to compile the OS X binaries… I’m not sure whether it’s enough to just submit the source the the CRAN.

Anyway, please check out my package and let me know if you encounter any problems or if you have suggestions on improving the documentation etc.

Open questions

  • How do I write an “ü” in the R-documentation (needed for my family name in author information)? The documentation is inside the R-files, the RD-files are created using Roxygen.
  • How do I include datasets inside an R-package? I would like to include an SPSS-dataset (.sav-File), so I can make the examples of my sji.XYZ functions running… (currently they’re outcommented so the package will compile and pass its check properly)
  • How to include a change log inside R-packages?

Visual interpretation of interaction terms in linear models with ggplot #rstats

I haven’t used interaction terms in (generalized) linear model quite often yet. However, recently I have had some situations where I tried to compute regression models with interaction terms and was wondering how to interprete the results. Just looking at the estimates won’t help much in such cases.

One approach used by some people is to compute the regressions with subgroups for each category of one interaction term. Let’s say predictor A has a 0/1 coding and predictor B is a continuous scale from 1 to 10, you fit a model for all cases with A=0 (hence excluding A from the model, no interaction of A and B), and for all cases with A=1 and compare the estimates of predictor B in each fitted model. This may give you an impression under which condition (i.e. in which subgroup) A has a stronger effect on B (higher interaction), but of course you don’t have the correct estimate values compared to a fitted model that includes both the interaction terms A and B.

Another approach is to calculate the results of y by hand, using the formula:
y = b0 + b1*predictorA + b2*predictorB + b3*predictorA*predictorB
This is quite complex and time-comsuming, especially if both predictors have several categories. However, this approach gives you a correct impression of the interaction between A and B. I investigated further on this topic and found this nice blogpost on interpreting interactions in regression (and a follow up), which explains very well how to calculate and interprete interaction terms.

Based on this knowledge, I thought of an automatization of calculating and visualizing interaction terms in linear models using R and ggplot.

Read on …

Plotting Likert-Scales (net stacked distributions) with ggplot #rstats

Update Thanks to Forrest for finding and fixing a bug. Scripts have been updated!

Update 2 Scripts have been updated because item ordering was still buggy. Hope everything is fixed now. Very helpful in this context was the new debug feature of RStudio, that also keeps track of all variables and their content and allows step-by-step execution of your code.

First of all, credits for this script must go to Ethan Brown, whose ideas for creating Likert scales like plots with ggplot built the core of the sjp.likert function in my package.

All I did was some visual tweaking like having positive percentage values on both sides of the x-axis, adding value labels and so on… You can pass a lot of different parameters to modify the graphical output. Please refer to my blog postings on R to get some impressions of how to tweak the plot (and/or look into the script header, which includes a description of all parameters).

Now to some examples. First, install the package from CRAN and load it with library(sjPlot). Then run following code:

likert_2 <- data.frame(as.factor(sample(1:2, 500, replace=T, prob=c(0.3,0.7))),
                       as.factor(sample(1:2, 500, replace=T, prob=c(0.6,0.4))),
                       as.factor(sample(1:2, 500, replace=T, prob=c(0.25,0.75))),
                       as.factor(sample(1:2, 500, replace=T, prob=c(0.9,0.1))),
                       as.factor(sample(1:2, 500, replace=T, prob=c(0.35,0.65))))
levels_2 <- list(c("Disagree", "Agree"))
items <- list(c("Q1", "Q2", "Q3", "Q4", "Q5"))
sjp.likert(likert_2, legendLabels=levels_2, axisLabels.x=items, orderBy="neg")

2-items Likert scale, ordered by "negative" categories.

2-items Likert scale, ordered by “negative” categories.


What you see above is a scale with two dimensions, ordered from highest “negative” category to lowest. If you leave out the orderBy parameter, the plot uses the normal item order:

likert_4 <- data.frame(as.factor(sample(1:4, 500, replace=T, prob=c(0.2,0.3,0.1,0.4))),
                       as.factor(sample(1:4, 500, replace=T, prob=c(0.5,0.25,0.15,0.1))),
                       as.factor(sample(1:4, 500, replace=T, prob=c(0.25,0.1,0.4,0.25))),
                       as.factor(sample(1:4, 500, replace=T, prob=c(0.1,0.4,0.4,0.1))),
                       as.factor(sample(1:4, 500, replace=T, prob=c(0.35,0.25,0.15,0.25))))
levels_4 <- list(c("Strongly disagree", "Disagree", "Agree", "Strongly Agree"))
items <- list(c("Q1", "Q2", "Q3", "Q4", "Q5"))
sjp.likert(likert_4, legendLabels=levels_4, axisLabels.x=items)

4-category-Likert-scale, ordered by items.

4-category-Likert-scale, ordered by items.


And finally, a plot with a different color set and items ordered from highest positive answer to lowest.

likert_6 <- data.frame(as.factor(sample(1:6, 500, replace=T, prob=c(0.2,0.1,0.1,0.3,0.2,0.1))),
                       as.factor(sample(1:6, 500, replace=T, prob=c(0.15,0.15,0.3,0.1,0.1,0.2))),
                       as.factor(sample(1:6, 500, replace=T, prob=c(0.2,0.25,0.05,0.2,0.2,0.2))),
                       as.factor(sample(1:6, 500, replace=T, prob=c(0.2,0.1,0.1,0.4,0.1,0.1))),
                       as.factor(sample(1:6, 500, replace=T, prob=c(0.1,0.4,0.1,0.3,0.05,0.15))))
levels_6 <- list(c("Very strongly disagree", "Strongly disagree", "Disagree", "Agree", "Strongly Agree", "Very strongly agree"))
items <- list(c("Q1", "Q2", "Q3", "Q4", "Q5"))
sjp.likert(likert_6, legendLabels=levels_6, barColor="brown", axisLabels.x=items, orderBy="pos")
6-category-Likert-scale with different color set and ordered by "positive" categories.

6-category-Likert-scale with different color set and ordered by “positive” categories.

If you need to plot stacked frequencies that have no “negative” and “positive”, but only one direction, you can also use the sjp.stackfrq function. Given that you use the likert-data frames from the above examples, you can run following code to plot stacked frequencies for scales that range from “low” to “high” and not from “negative” to “positive”.

levels_42 <- list(c("Independent", "Slightly dependent", "Dependent", "Severely dependent"))
levels_62 <- list(c("Independent", "Slightly dependent", "Dependent", "Very dependent", "Severely dependent", "Very severely dependent"))
sjp.stackfrq(likert_4, legendLabels=levels_42, axisLabels.x=items)
sjp.stackfrq(likert_6, legendLabels=levels_62, axisLabels.x=items)

This produces following two plots:

Stacked frequencies of 4-category-items.

Stacked frequencies of 4-category-items.


Stacked frequencies of 6-category-items.

Stacked frequencies of 6-category-items.

That’s it!

Plotting principal component analysis with ggplot #rstats

This script was almost written on parallel to the sjPlotCorr script because it uses a very similar ggplot-base. However, there’s also a very nice posting over at Martin’s Bio Blog which show alternative approaches on plotting PCAs.

Anyway, if you download the sjPlotPCA.R script, you can easily plot a PCA with varimax rotation like this:

likert_4 <- data.frame(sample(1:4, 500, replace=T, prob=c(0.2,0.3,0.1,0.4)),
                       sample(1:4, 500, replace=T, prob=c(0.5,0.25,0.15,0.1)),
                       sample(1:4, 500, replace=T, prob=c(0.4,0.15,0.25,0.2)),
                       sample(1:4, 500, replace=T, prob=c(0.25,0.1,0.4,0.25)),
                       sample(1:4, 500, replace=T, prob=c(0.1,0.4,0.4,0.1)),
                       sample(1:4, 500, replace=T,),
                       sample(1:4, 500, replace=T, prob=c(0.35,0.25,0.15,0.25)))
colnames(likert_4) <- c("V1", "V2", "V3", "V4", "V5", "V6", "V7")
source("../lib/sjPlotPCA.R")
sjp.pca(likert_4)

So, all you have to do is creating a data frame where each column represents one variable / case and pass this data frame to the function. This will result in something like this:

PCA of 7 variables resulting in 3 extracted factors. Cronbach's Alpha value of each "factor scale" printed at bottom.

PCA of 7 variables resulting in 3 extracted factors (varimax rotation). Cronbach’s Alpha value of each “factor scale” printed at bottom.

The script automatically calculates the Cronbach’s Alpha value for each “factor scale”, assuming that the variables with the highest factor loading belongs to this factor. The amount of factors is calculated according to the Kaiser criterion. You can also create a plot of this calcuation by setting the parameter plotEigenvalues=TRUE.

The next small example shows two plots and uses a computed PCA as paramater:

pca <- prcomp(na.omit(likert_4), retx=TRUE, center=TRUE, scale.=TRUE)
sjp.pca(pca, plotEigenvalues=TRUE, type="circle")

Eigenvalue plot determining amount of factors (Kaiser criterion)

Eigenvalue plot determining amount of factors (Kaiser criterion)


Same PCA plot as above, with PCA object instead of data frame as parameter.

Same PCA plot as above, with PCA object instead of data frame as parameter.


Note that when using a PCA object as parameter and no data frame, the Cronbach’s Alpha value cannot be calculated.

That’s it! The source is available on my download page.

Examples for sjPlotting functions, including correlations and proportional tables with ggplot #rstats

Sometimes people ask me how the examples of my plotting functions I show here can be reproduced without having a SPSS data set (or at least, without having the data set I use because it’s not public yet). So I started to write some examples that run “out of the box” and which I want to present you here. Furthermore, two new plotting functions are introduced: plotting correlations and plotting proportional tables on a percentage scale.

As always, you can find the latest version of my R scripts on my download page.

Following plotting functions will be described in this posting:

  • Plotting proportional tables: sjPlotPropTable.R
  • Plotting correlations: sjPlotCorr.R
  • Plotting frequencies: sjPlotFrequencies.R
  • Plotting grouped frequencies: sjPlotGroupFrequencies.R
  • Plotting linear model: sjPlotLinreg.R
  • Plotting generalized linear models: sjPlotOdds.R

Please note that I have changed function and parameter names in order to have consistent, logical names across all functions!

At the end of this posting you will find some explanation on the different parameters that allow you to fit the plotting results to your needs…

Continue reading this post…