Skip to content

Latest commit

 

History

History
293 lines (164 loc) · 9.66 KB

README.md

File metadata and controls

293 lines (164 loc) · 9.66 KB

ROBUSTA

Author: Eitan Hemed

robusta is a statistics package in Python3 providing an interface to many common statistical analyses, performed using through R and RPY2.

PLEASE NOTE robusta is under active development and is supplied as-is with no guarantees.

Installation

Install with pip using pip install robusta-stats, see also Installation.

Documentation

See here.

Usage

For the most recent, thorough tutorial in the different features of robusta, head on to Google Colab.

Some of the features are shown below.

Importing the library and loading data

This could take ~15 seconds as many R libraries are imported under the hood. If you begin with an empty R environment the first you import robusta should take at least a couple of minutes, as R dependencies will be installed.

import robusta as rst

First, define a helper function used to pretty-print output of dataframes when converting the notebook to .md (credit).

from tabulate import tabulate
import IPython.display as d

def md_print_df(df):
    md = tabulate(df, headers='keys', tablefmt='pipe')
    md = md.replace('|    |','| %s |' % (df.index.name if df.index.name else ''))
    return d.Markdown(md)

First off, we need data. Using robusta we can import R built-in and some imported datasets. You can get a full list of the datasets, similarly to calling to data() with no input arguments in R.

md_print_df(rst.get_available_datasets().tail())
Package Item Description
284 ARTool Higgins1990Table5 Split-plot Experiment Examining Effect of Moisture and Fertilizer on Dry Matter in Peat Pots
285 ARTool Higgins1990Table1.art Aligned Rank Transformed Version of Higgins1990Table1
286 ARTool Higgins1990Table1 Synthetic 3x3 Factorial Randomized Experiment
287 ARTool ElkinABC Synthetic 2x2x2 Within-Subjects Experiment
288 ARTool ElkinAB Synthetic 2x2 Within-Subjects Experiment

We can import a dataset using rst.load_dataset

iris = rst.load_dataset('iris')
md_print_df(iris.head())
dataset_rownames Sepal.Length Sepal.Width Petal.Length Petal.Width Species
0 1 5.1 3.5 1.4 0.2 setosa
1 2 4.9 3 1.4 0.2 setosa
2 3 4.7 3.2 1.3 0.2 setosa
3 4 4.6 3.1 1.5 0.2 setosa
4 5 5 3.6 1.4 0.2 setosa

Running statistical analyses

Analyses are performed through using designated model objects that also store the . The model objects are returned through calls to the function API. In this example we create a model (m) object by calling t2samples. m will be used to fit the statistical model, returning the results object.

Here is a paired-samples t-test using the Students' sleep dataset previously loaded:

# Create the model
m = rst.groupwise.T2Samples(
    data=rst.load_dataset('sleep'), independent='group', 
    dependent='extra', subject='ID', paired=True, tail='less')

# Dataframe format of the results
md_print_df(m.report_table())
t df p-value Cohen-d Low Cohen-d Cohen-d High
1 -4.06213 9 0.00141645 -2.11801 -1.28456 -0.414622
# Textual report of the results - copy and paste into your results section!
m.report_text()
't(9) = -4.06, p = 0.001'

We can reset the models in order to update the model parameters and re-fit it. In this example, we run the same model an an independent samples t-test:

m.reset(paired=False, assume_equal_variance=True, refit=True)
md_print_df(m.report_table())
t df p-value Cohen-d Low Cohen-d Cohen-d High
1 -1.86081 18 0.0395934 -1.73882 -0.832181 0.0954595

Bayesian t-tests

bayes_t2samples and bayes_t1sample allow you to calculate Bayes factors or sample from the posterior distribution:

m = rst.groupwise.BayesT2Samples(
        data=rst.load_dataset('mtcars'), subject='dataset_rownames',
        dependent='mpg', independent='am', prior_scale=0.5,
        paired=False)

md_print_df(m.report_table())
model bf error
0 Alt., r=0.5 71.3861 7.97835e-07
# Test different null intervals and prior values:
m.reset(prior_scale=0.1, null_interval=[0, 0.5], refit=True)
print(f'{m.report_text()}\n\n')
md_print_df(m.report_table())
Alt., r=0.1 [BF1:0 = 18.64, Error = 0.001%]
model bf error
0 Alt., r=0.1 18.6411 2.33663e-05

Analysis of variance

use Anova to run between, within or mixed-design ANOVA, we load the anxiety dataset for the next demonstrations.

For non-parametric ANOVAs see KruskalWallisTest, FriedmanTest and AlignedRanksTest

# Load the dataset and modify it from a 'wide' to 'long' format dataframe
anxiety = rst.load_dataset('anxiety').set_index(['id', 'group']
                                           ).filter(regex='^t[1-3]$').stack().reset_index().rename(
    columns={0: 'score',
             'level_2': 'time'})
md_print_df(anxiety.head())
id group time score
0 1 grp1 t1 14.1
1 1 grp1 t2 14.4
2 1 grp1 t3 14.1
3 2 grp1 t1 14.5
4 2 grp1 t2 14.6
m = rst.groupwise.Anova(
        data=anxiety, subject='id',
        dependent='score', between='group', within='time')
md_print_df(m.report_table())
R[write to console]: Contrasts set to contr.sum for the following variables: group
Term p-value Partial Eta-Squared F df1 df2
1 group 0.019 0.172 4.35 2 42
2 time 0.001 0.904 394.91 1.79 75.24
3 group:time 0.001 0.84 110.19 3.58 75.24

Similarly, we run the model usign only the between subject term (group). As the model was already generated we can simpyl drop the within-subject term:

m.reset(within=None, refit=True)
md_print_df(m.report_table())
R[write to console]: Contrasts set to contr.sum for the following variables: group
Term p-value Partial Eta-Squared F df1 df2
1 group 0.019 0.172 4.35 2 42

R and many other statistical packages (e.g., statsmodels support a formula interface to fit statistical models. Here it is shown that a model can also be specified by the formula kwargs rather than specifying dependent, between etc. The formula indicates that the score column is regressed by the time variable, with observations nested within the id column.

m.reset(formula='score~time|id', refit=True)
md_print_df(m.report_table())
Term p-value Partial Eta-Squared F df1 df2
1 time 0.001 0.601 66.23 1.15 50.55

We can also run a similar, bayesian ANOVA using BayesAnova comparing the specified terms to the null model:

m = rst.groupwise.BayesAnova(data=anxiety, within='time',
                             dependent='score', subject='id')
md_print_df(m.report_table())
model bf error
0 time 496.129 7.82496e-05

Work in progress and planned features

robusta includes several other features that are either under development or planned for the future.

Currently under work

  • Regressions and correlations modules

Planned

  • Sequential analysis plots (inspired by JASP)

How to contribute

All help is welcome, but currently there are no specific guidelines. Please contact Eitan Hemed