Skip to content

Instantly share code, notes, and snippets.

@mwaskom

mwaskom/lmer_foo.ipynb

Created September 27, 2012 05:02
Show Gist options
  • Save mwaskom/3792250 to your computer and use it in GitHub Desktop.
Save mwaskom/3792250 to your computer and use it in GitHub Desktop.
Download ZIP
Playing with lmer in ipynb
Display the source blob
Display the rendered blob
Raw
lmer_foo.ipynb
{
"metadata": {
"name": "lmer_foo"
},
"nbformat": 3,
"nbformat_minor": 0,
"worksheets": [
{
"cells": [
{
"cell_type": "code",
"collapsed": true,
"input": [
"%load_ext rmagic\n",
"%%R library(lme4)"
],
"language": "python",
"metadata": {},
"outputs": [
{
"output_type": "display_data",
"text": [
"Loading required package: Matrix\n",
"Loading required package: lattice\n",
"\n",
"Attaching package: 'Matrix'\n",
"\n",
"The following object(s) are masked from 'package:base':\n",
"\n",
" det\n",
"\n",
"\n",
"Attaching package: 'lme4'\n",
"\n",
"The following object(s) are masked from 'package:stats':\n",
"\n",
" AIC, BIC\n",
"\n",
"Warning messages:\n",
"1: package 'lme4' was built under R version 2.13.2 \n",
"2: package 'Matrix' was built under R version 2.13.2 \n"
]
}
],
"prompt_number": 2
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"n_subs = 30\n",
"intercepts = 10 + 2 * randn(n_subs)\n",
"slopes = 1.5 + randn(n_subs)\n",
"dv = concatenate([intercepts, intercepts + slopes])\n",
"iv = repeat([0, 1], n_subs)\n",
"subj = tile([\"s%02d\" % i for i in range(n_subs)], 2)"
],
"language": "python",
"metadata": {},
"outputs": [],
"prompt_number": 40
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"from scipy.stats import zscore\n",
"dv_z = zeros_like(dv)\n",
"for s in range(n_subs):\n",
" idx = subj == \"s%02d\" % s\n",
" dv_z[idx] = zscore(dv[idx])"
],
"language": "python",
"metadata": {},
"outputs": [],
"prompt_number": 42
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"%Rpush dv dv_z iv subj"
],
"language": "python",
"metadata": {},
"outputs": [],
"prompt_number": 48
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"%%R \n",
"m = lm(dv ~ iv)\n",
"print(summary(m))\n",
"print(logLik(m))"
],
"language": "python",
"metadata": {},
"outputs": [
{
"output_type": "display_data",
"text": [
"\n",
"Call:\n",
"lm(formula = dv ~ iv)\n",
"\n",
"Residuals:\n",
" Min 1Q Median 3Q Max \n",
"-4.050 -1.608 -0.028 1.377 5.292 \n",
"\n",
"Coefficients:\n",
" Estimate Std. Error t value Pr(>|t|) \n",
"(Intercept) 9.9152 0.3891 25.486 <2e-16 ***\n",
"iv 1.2769 0.5502 2.321 0.0238 * \n",
"---\n",
"Signif. codes: 0 \u2018***\u2019 0.001 \u2018**\u2019 0.01 \u2018*\u2019 0.05 \u2018.\u2019 0.1 \u2018 \u2019 1 \n",
"\n",
"Residual standard error: 2.131 on 58 degrees of freedom\n",
"Multiple R-squared: 0.08498,\tAdjusted R-squared: 0.0692 \n",
"F-statistic: 5.386 on 1 and 58 DF, p-value: 0.02383 \n",
"\n",
"'log Lik.' -129.5124 (df=3)\n"
]
}
],
"prompt_number": 59
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"%%R \n",
"m = lmer(dv ~ iv + (1 + iv | subj))\n",
"print(m)\n",
"print(logLik(m))"
],
"language": "python",
"metadata": {},
"outputs": [
{
"output_type": "display_data",
"text": [
"Linear mixed model fit by REML \n",
"Formula: dv ~ iv + (1 + iv | subj) \n",
" AIC BIC logLik deviance REMLdev\n",
" 232.4 244.9 -110.2 218.9 220.4\n",
"Random effects:\n",
" Groups Name Variance Std.Dev. Corr \n",
" subj (Intercept) 3.59278 1.89546 \n",
" iv 0.17659 0.42022 0.362 \n",
" Residual 0.57116 0.75575 \n",
"Number of obs: 60, groups: subj, 30\n",
"\n",
"Fixed effects:\n",
" Estimate Std. Error t value\n",
"(Intercept) 9.9152 0.3726 26.61\n",
"iv 1.2769 0.2097 6.09\n",
"\n",
"Correlation of Fixed Effects:\n",
" (Intr)\n",
"iv -0.121\n",
"'log Lik.' -110.1846 (df=6)\n"
]
}
],
"prompt_number": 60
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"%%R\n",
"m = lm(dv_z ~ iv)\n",
"print(summary(m))\n",
"print(logLik(m))"
],
"language": "python",
"metadata": {},
"outputs": [
{
"output_type": "display_data",
"text": [
"\n",
"Call:\n",
"lm(formula = dv_z ~ iv)\n",
"\n",
"Residuals:\n",
" Min 1Q Median 3Q Max \n",
" -1.6 -0.4 0.0 0.4 1.6 \n",
"\n",
"Coefficients:\n",
" Estimate Std. Error t value Pr(>|t|) \n",
"(Intercept) -0.6000 0.1486 -4.039 0.00016 ***\n",
"iv 1.2000 0.2101 5.712 4.05e-07 ***\n",
"---\n",
"Signif. codes: 0 \u2018***\u2019 0.001 \u2018**\u2019 0.01 \u2018*\u2019 0.05 \u2018.\u2019 0.1 \u2018 \u2019 1 \n",
"\n",
"Residual standard error: 0.8137 on 58 degrees of freedom\n",
"Multiple R-squared: 0.36,\tAdjusted R-squared: 0.349 \n",
"F-statistic: 32.62 on 1 and 58 DF, p-value: 4.049e-07 \n",
"\n",
"'log Lik.' -71.7477 (df=3)\n"
]
}
],
"prompt_number": 61
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"%%R\n",
"m = lmer(dv_z ~ iv + (1 | subj) + (iv -1 | subj), REML=F)\n",
"print(m)\n",
"print(logLik(m))"
],
"language": "python",
"metadata": {},
"outputs": [
{
"output_type": "display_data",
"text": [
"Linear mixed model fit by maximum likelihood \n",
"Formula: dv_z ~ iv + (1 | subj) + (iv - 1 | subj) \n",
" AIC BIC logLik deviance REMLdev\n",
" 153.5 164 -71.75 143.5 147.5\n",
"Random effects:\n",
" Groups Name Variance Std.Dev.\n",
" subj (Intercept) 0.00 0.0 \n",
" subj iv 0.00 0.0 \n",
" Residual 0.64 0.8 \n",
"Number of obs: 60, groups: subj, 30\n",
"\n",
"Fixed effects:\n",
" Estimate Std. Error t value\n",
"(Intercept) -0.6000 0.1461 -4.108\n",
"iv 1.2000 0.2066 5.809\n",
"\n",
"Correlation of Fixed Effects:\n",
" (Intr)\n",
"iv -0.707\n",
"'log Lik.' -71.7477 (df=5)\n"
]
}
],
"prompt_number": 66
},
{
"cell_type": "code",
"collapsed": false,
"input": [],
"language": "python",
"metadata": {},
"outputs": []
}
],
"metadata": {}
}
]
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment