rePCA_demo.ipynb

{
  "cells": [
    {
      "cell_type": "markdown",
      "source": [
        "# rePCA() Demo\n### Reinhold Kliegl\n### 2020-02-26\n## Background\n\nIn the R-based `lme4` package, `rePCA()` results are not invariant with linear transformation of predictor\nif based on (relative) covariance matrix. In the Julia-based `MixedModels` package, by default they are based\non the correlation, not the covariance matrix. The covariance-based version is available with \n`MixedModels.PCA(model, corr=false)`. \n\n## Packages"
      ],
      "metadata": {}
    },
    {
      "outputs": [],
      "cell_type": "code",
      "source": [
        "using DataFrames, DataFramesMeta, MixedModels, RCall\nusing StatsBase, StatsModels, BenchmarkTools"
      ],
      "metadata": {},
      "execution_count": null
    },
    {
      "cell_type": "markdown",
      "source": [
        "## Sleepstudy"
      ],
      "metadata": {}
    },
    {
      "outputs": [],
      "cell_type": "code",
      "source": [
        "MixedModels.datasets()\n\nsleepstudy =  MixedModels.dataset(:sleepstudy)\n\ncolnames = [\"Subj\", \"days\", \"reaction\"]\nrename!(sleepstudy, Symbol.(colnames))\n\nsleepstudy = @linq sleepstudy |>\n             transform(days2 = :days .*2,\n                       days3 = :days ./2);"
      ],
      "metadata": {},
      "execution_count": null
    },
    {
      "cell_type": "markdown",
      "source": [
        "## LMMs\n\nLMM `fm2` and `fm3` use a linear transformation of `days`"
      ],
      "metadata": {}
    },
    {
      "outputs": [],
      "cell_type": "code",
      "source": [
        "fm1 = fit(LinearMixedModel, @formula(reaction ~ 1 + days  + (1 + days | Subj)), sleepstudy)\nfm2 = fit(LinearMixedModel, @formula(reaction ~ 1 + days2 + (1 + days2 | Subj)), sleepstudy)\nfm3 = fit(LinearMixedModel, @formula(reaction ~ 1 + days3 + (1 + days3 | Subj)), sleepstudy)"
      ],
      "metadata": {},
      "execution_count": null
    },
    {
      "cell_type": "markdown",
      "source": [
        "## Inspect rePCA property \n\nWhen we use `corr=true` -- normalized cumulative proportion of variance is the same"
      ],
      "metadata": {}
    },
    {
      "outputs": [],
      "cell_type": "code",
      "source": [
        "fm1.rePCA"
      ],
      "metadata": {},
      "execution_count": null
    },
    {
      "outputs": [],
      "cell_type": "code",
      "source": [
        "fm2.rePCA"
      ],
      "metadata": {},
      "execution_count": null
    },
    {
      "outputs": [],
      "cell_type": "code",
      "source": [
        "fm3.rePCA"
      ],
      "metadata": {},
      "execution_count": null
    },
    {
      "cell_type": "markdown",
      "source": [
        "Equivalent info for `MixedModels.PCA(fm, corr=true)`"
      ],
      "metadata": {}
    },
    {
      "outputs": [],
      "cell_type": "code",
      "source": [
        "MixedModels.PCA(fm1, corr=true)"
      ],
      "metadata": {},
      "execution_count": null
    },
    {
      "outputs": [],
      "cell_type": "code",
      "source": [
        "MixedModels.PCA(fm2, corr=true)"
      ],
      "metadata": {},
      "execution_count": null
    },
    {
      "outputs": [],
      "cell_type": "code",
      "source": [
        "MixedModels.PCA(fm3, corr=true)"
      ],
      "metadata": {},
      "execution_count": null
    },
    {
      "cell_type": "markdown",
      "source": [
        "... but not for `MixedModels.PCA(fm, corr=false)`"
      ],
      "metadata": {}
    },
    {
      "outputs": [],
      "cell_type": "code",
      "source": [
        "MixedModels.PCA(fm1, corr=false)"
      ],
      "metadata": {},
      "execution_count": null
    },
    {
      "outputs": [],
      "cell_type": "code",
      "source": [
        "MixedModels.PCA(fm2, corr=false)"
      ],
      "metadata": {},
      "execution_count": null
    },
    {
      "outputs": [],
      "cell_type": "code",
      "source": [
        "MixedModels.PCA(fm3, corr=false)"
      ],
      "metadata": {},
      "execution_count": null
    },
    {
      "cell_type": "markdown",
      "source": [
        "The normalized cumulative proportions of variance reported for `corr=false` \nagree with results from `lme4::rePCA()` as shown below. \n\n## MixedModel.PCA() default options \n\nThe `show()` provides additional information. Here are the options.\n\n```\nBase.show(pca::PCA;\n          ndigitsmat=2, ndigitsvec=2, ndigitscum=4,\n          covcor=true, loadings=true, variances=false, stddevs=false\n```\n\nFor example, to obtain only inequivalent cumulative proportions of variance."
      ],
      "metadata": {}
    },
    {
      "outputs": [],
      "cell_type": "code",
      "source": [
        "fm1_pca_f = MixedModels.PCA(fm1, corr=false);\nshow(fm1_pca_f.Subj, covcor=false, loadings=false, variances=false, stddevs=false)  \n# Note: Must add grouping identifier!\n\nfm2_pca_f = MixedModels.PCA(fm2, corr=false);\nshow(fm2_pca_f.Subj, covcor=false, loadings=false, variances=false, stddevs=false)\n\nfm3_pca_f = MixedModels.PCA(fm3, corr=false);\nshow(fm3_pca_f.Subj, covcor=false, loadings=false, variances=false, stddevs=false)"
      ],
      "metadata": {},
      "execution_count": null
    },
    {
      "cell_type": "markdown",
      "source": [
        "## R-based lme4::rePCA() results"
      ],
      "metadata": {}
    },
    {
      "outputs": [],
      "cell_type": "code",
      "source": [
        "R\"\"\"\nlibrary(lme4)\n\nDays2 <- sleepstudy$Days*2\nDays3 <- sleepstudy$Days/2\n\nfm1 <- lmer(Reaction ~ 1 + Days + (1 + Days | Subject), sleepstudy, REML=FALSE)\nfm2 <- lmer(Reaction ~ Days2 + (1 + Days2 | Subject), sleepstudy, REML=FALSE)\nfm3 <- lmer(Reaction ~ Days3 + (1 + Days3 | Subject), sleepstudy, REML=FALSE)\n\nprint(summary(rePCA(fm1)))\nprint(summary(rePCA(fm2)))  \nprint(summary(rePCA(fm3)))\n\"\"\";"
      ],
      "metadata": {},
      "execution_count": null
    }
  ],
  "nbformat_minor": 2,
  "metadata": {
    "language_info": {
      "file_extension": ".jl",
      "mimetype": "application/julia",
      "name": "julia",
      "version": "1.3.1"
    },
    "kernelspec": {
      "name": "julia-1.3",
      "display_name": "Julia 1.3.1",
      "language": "julia"
    }
  },
  "nbformat": 4
}