SocialTV.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Social facilitation of laughter and smiles in preschool children\n",
    "## Data analysis \n",
    "\n",
    "#### Caspar Addyman, Charlotte Folquist, Lenka Levakova, Sarah Rees\n",
    "#### Birkbeck University of London\n",
    "\n",
    "#### Abstract\n",
    "\n",
    "Abstract\n",
    "======\n",
    "Surprisingly little research has investigated the social dimension of laughter in preschool children.  This experiment studied children’s responses to amusing video clips in the presence or absence of peers. The sample consisted of 9 boys and 11 girls aged 31-49 months (mean 39.8, SD 4.2) who watched three cartoons under three different conditions: individually, in pairs, or in groups of 6 or 8. The social viewing conditions showed significantly higher numbers of laughs and smiles than the individual viewing condition. On average children laughed eight times more in company as on their own, and smiled almost three times as much. No differences were found between pairs and groups, and no association was found between subjective funniness ratings and group size. This suggests that the presence of even a single social partner can change behaviour in response to humorous material. It supports the idea that laughter and smiles are primarily flexible social signals rather than reflexive responses to humour.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# necessary packages\n",
    "install.packages(c('readr','reshape2','ez','effects','ggplot2','ggsignif', 'Rmisc', 'svglite','pwr','sjmisc'), repos='http://cran.us.r-project.org')\n",
    "library(tidyverse)\n",
    "library(readxl) \n",
    "library(ggplot2)\n",
    "library(ggsignif)\n",
    "library(reshape2)\n",
    "library(ez)\n",
    "library(effects)\n",
    "library(Rmisc)\n",
    "library(svglite)\n",
    "library(pwr)\n",
    "library(sjPlot)\n",
    "library(sjmisc)\n",
    "library()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Power Calculations\n",
    "\n",
    "Some useful information found here\n",
    "http://r-video-tutorial.blogspot.co.uk/2017/07/power-analysis-and-sample-size.html"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#first write a small helper function to calculate partial eta squared..\n",
    "eta2_p <- function(dfeffect,dferror,FVAL) {\n",
    "  #If we have a reported effect\n",
    "  # F(dfeffect,dferror) = FVAL\n",
    "  # then partial eta square is calculated as\n",
    "  \n",
    "  eta = (FVAL * dfeffect) / (FVAL * dfeffect + dferror) \n",
    "  # see Lakens(2013) doi:10.3389/fpsyg.2013.00863\n",
    "  return(eta)\n",
    "}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We calculate the power based on the effects found in the original Chapman (1973) study of social laughter in 7 & 8 year olds.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#What was the effect size in Chapman (1973)?\n",
    "#For laughter F( 2,74) = 27.86, p < .005]\n",
    "EtaSQ = eta2_p(2,74,27.86) # returns 0.4295405\n",
    "\n",
    "#For smiling [ F( 2,74 ) = 49.43\n",
    "eta2_p(2,74,49.43) # returns 0.5719079\n",
    "\n",
    "##To calculate the group size we convert etasq into f squared aka f2\n",
    "Chapmanf2 = EtaSQ / (1-EtaSQ) \n",
    "\n",
    "##the pwr package lets us calculate the group size \n",
    "pwr.f2.test(u=2, f2=Chapmanf2, sig.level=0.05,power = .80)\n",
    "\n",
    "#this gives v = 13.2 which implies a mimumum group size\n",
    "# of 17 (so that degrees of freedom is greater than 13.2)\n",
    "#alternatively we have chosen a group size of 20 then our predicted power is\n",
    "pwr.f2.test(u=2, v=17, f2=Chapmanf2, sig.level=0.05)\n",
    "\n",
    "#If we were to use a more conservative general large effect size\n",
    "# suggested by Cohen (1988) as provided by the cohen.ES function\n",
    "# then we would need a larger sample\n",
    "cohenLarge = cohen.ES(test=\"f2\",size=\"large\")$effect.size\n",
    "pwr.f2.test(u=2, f2=cohenLarge, sig.level=0.05,power = .80)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Load in the data\n",
    "First load in the raw data from the spreadsheet. \n",
    "Note, to make this work on your machine you will need to set the working directory appropriatelty with `setwd`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#load raw data\n",
    "setwd('C:\\\\Users\\\\cas\\\\Dropbox\\\\projects\\\\Baby Laughter\\\\SocialTV')\n",
    "alldata<- read_excel('SocialTV.xlsx', 'RatingData')\n",
    "\n",
    "#specify which variables are factors\n",
    "alldata$SEX <- factor(alldata$SEX)\n",
    "alldata$ORDER <- factor(alldata$ORDER)\n",
    "#have a quick look at the data \n",
    "head(alldata)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Basic descriptive stats & correlations"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "##############################\n",
    "# descriptive stats & correlations\n",
    "##############################\n",
    "\n",
    "#Descriptive stats\n",
    "mean(alldata$`Laughs Groups`)\n",
    "sd(alldata$`Laughs Groups`)\n",
    "mean(alldata$`Laughs Pairs`)\n",
    "sd(alldata$`Laughs Pairs`)\n",
    "mean(alldata$`Laughs Indiv`)\n",
    "sd(alldata$`Laughs Indiv`)\n",
    "\n",
    "# we use sjmisc to output descriptive stats table\n",
    "descr(alldata[,c(3,6:11)],out='txt')\n",
    "\n",
    "# use sjPlot to output a nicely formatted correlation table.\n",
    "sjt.corr(alldata[,c(3,6:11)],p.numeric = TRUE, corr.method = \"pearson\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Some helper functions to display the graphs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "##Social TV paper Addyman et al 2017\n",
    "##Some helper functions to display graphs.\n",
    "\n",
    "#' Takes the wide (1 row per person) version of the data and makes a long \n",
    "#' subset with a single DV (1 col for grouping variable, 1 col for DV)\n",
    "#' In output grouping variable is \"groupVar\", DV is \"value\"\n",
    "#' \n",
    "#' @param allData original data\n",
    "#' @param col1 1st DV column, etc\n",
    "groupLongData <- function(allData,col1,col2,col3){\n",
    "  #subset and rearrange laugh data for anova and plotting\n",
    "  #convert to long format \n",
    "  data.long <- melt(allData,\n",
    "                    id.vars = c(\"ID\",\"SEX\",\"ORDER\"),\n",
    "                    measure.vars = c(col1, col2,col3),\n",
    "                    variable.name =\"groupVar\")\n",
    "  #indicate the factor variable\n",
    "  data.long$groupVar <- factor(data.long$groupVar)\n",
    "  return(data.long)\n",
    "}\n",
    "\n",
    "#' Takes three columns from dataframe and compares them in one way Anova\n",
    "#' while also plotting a box plot with overlayed individual data points.\n",
    "#' \n",
    "#' @param data.Long Long format version of the data.\n",
    "#' @param varName Name of the grouping Variable.\n",
    "#' @examples\n",
    "#' plotGroupData(longLaughs, \"Smiles\")\n",
    "plotGroupData <- function (data.long, varName) {\n",
    "  #plot groups as box plot and individual values\n",
    "  ylabel<-paste(\"Number of \", varName,sep = \"\")\n",
    "  p<-ggplot(data=data.long,aes(groupVar,value)) +\n",
    "    geom_boxplot() +\n",
    "    geom_point(aes(groupVar,value), size=3,position = position_jitter(width = .1) ) +\n",
    "    geom_signif(comparisons = list(c(2,3),c(1,3),c(1,2)), step_increase = .12, map_signif_level = TRUE, test=\"t.test\", test.args = list(paired = TRUE)) +\n",
    "    ylab(ylabel) +\n",
    "    xlab(NULL) +\n",
    "    theme_bw(base_size =18)\n",
    "  return(p)\n",
    "}\n",
    "\n",
    "#' plots a standard interaction plot for one-way anova but taking care to\n",
    "#' use standard errors suitable for repeated errors (like SPSS)\n",
    "#' \n",
    "#' @param data.Long Long format version of the data.\n",
    "#' @param varName Name of the grouping Variable.\n",
    "#' @examples\n",
    "#' plotOneWayANOVA(groupLaughter, \"Smiles\") \n",
    "plotOneWayANOVA <- function (longData, varName) {\n",
    "  #calculate repeated measures standard errors\n",
    "  laughsummary<-summarySEwithin(longData, measurevar = \"value\",betweenvars = NULL,withinvars =  \"groupVar\")\n",
    "  # Standard error of the mean\n",
    "  ylabel<-paste(\"Number of \", varName,sep = \"\")\n",
    "  p<-ggplot(laughsummary, aes(x=groupVar, y=value, group=1)) + \n",
    "    geom_errorbar(aes(ymin=value-se, ymax=value+se), width=.1) +\n",
    "    geom_signif(comparisons = list(c(2,3),c(1,3),c(1,2)), step_increase = .12, map_signif_level = TRUE, test=\"t.test\", test.args = list(paired = TRUE)) +\n",
    "    geom_line() +\n",
    "    geom_point(size=3) +\n",
    "    ylab(ylabel) +\n",
    "    xlab(NULL) +\n",
    "    theme_bw(base_size =18)\n",
    "  return(p)\n",
    "}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Number of laughs per group\n",
    "\n",
    "We need to reshape the data before we plot and analyse it. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "#To visualise the data we need to reorder it. \n",
    "longLaughs<-groupLongData(alldata, \"Laughs Groups\", \"Laughs Pairs\", \"Laughs Indiv\")\n",
    "\n",
    "#plot data in two different ways \n",
    "p1<-plotGroupData(longLaughs,\"Laughs\")\n",
    "p2<-plotOneWayANOVA(longLaughs,\"Laughs\")\n",
    "show(p1)\n",
    "show(p2)\n",
    "##show plots side by side\n",
    "#grid.arrange(p1, p2, ncol=2)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### One Way ANOVA\n",
    "Now we run the stats. First we run a 3X3 ANOVA with presentation condition as the factor at three levels (Groups, Pairs, Indiv) and order at three levels (A,B,C) and include sex as a factor\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#initial analysis including sex as variable\n",
    "ezANOVA(longLaughs,\n",
    "        dv = value,\n",
    "        wid = ID,\n",
    "        between = .(ORDER, SEX),\n",
    "        within = groupVar )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Next run the analysis from the paper without sex as a factor"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# 3x3 mixed ANOVA from paper \n",
    "ezANOVA(longLaughs,\n",
    "        dv = value,\n",
    "        wid = ID,\n",
    "        between = .(ORDER),\n",
    "        within = groupVar )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Planned comparisons"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#run pairwise t-tests\n",
    "t.test(alldata$`Laughs Groups`,alldata$`Laughs Pairs`, paired = TRUE)\n",
    "t.test(alldata$`Laughs Groups`,alldata$`Laughs Indiv`, paired = TRUE)\n",
    "t.test(alldata$`Laughs Pairs`, alldata$`Laughs Indiv`, paired = TRUE)\n",
    "\n",
    "#but we'd like bonferoni corrected p-values\n",
    "pairwise.t.test(longLaughs$value,longLaughs$groupVar , \n",
    "                p.adjust.method =  \"bonferroni\", \n",
    "                paired = TRUE)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Number of Smiles\n",
    "\n",
    "Next we run the same analyses for the number of smiles per condition."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "longSmiles<-groupLongData(alldata, \"Smiles Groups\", \"Smiles Pairs\", \"Smiles Indiv\")\n",
    "plotGroupData(longSmiles,\"Smiles\")\n",
    "plotOneWayANOVA(longSmiles,\"Smiles\")\n",
    "ezANOVA(longSmiles,\n",
    "        dv = value,\n",
    "        wid = ID,\n",
    "        within = groupVar )\n",
    "t.test(alldata$`Smiles Groups`,alldata$`Smiles Pairs`, paired = TRUE)\n",
    "t.test(alldata$`Smiles Groups`,alldata$`Smiles Indiv`, paired = TRUE)\n",
    "t.test(alldata$`Smiles Pairs`, alldata$`Smiles Indiv`, paired = TRUE)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Combined number of laughs and smiles\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "longBoth<-groupLongData(allData, \"SmilesLaughs Groups\", \"SmilesLaughs Pairs\", \"SmilesLaughs Indiv\")\n",
    "plotGroupData(longBoth,\"Laughs and Smiles\")\n",
    "plotOneWayANOVA(longBoth,\"Laughs and Smiles\")\n",
    "ezANOVA(longBoth,\n",
    "        dv = value,\n",
    "        wid = ID,\n",
    "        within = groupVar )\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Ratings of funniness\n",
    "\n",
    "After each video each child was asked to rate the video as \"Very Funny\", \"Quite Funny\" or \"Not Funny\" with a verbal request and using the following visual rating scale:\n",
    "<img src=\"funninessRating.png\" height=\"450\" width=\"260\">\n",
    "\n",
    "We were interested in whether their ratings differed according to the conditions in which they viewed the video or how much they laughed."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#But do they find it funny?\n",
    "#Let's look at the ratings?\n",
    "\n",
    "longRatings<-groupLongData(alldata,\"Funniness Groups\", \"Funniness Pairs\", \"Funniness Indiv\")\n",
    "#let R know this is categorical data\n",
    "longRatings$value<-factor(longRatings$value)\n",
    "#lets collect responses into table (pick the appropriate columns)\n",
    "tableRatings<-table(longRatings[,4:5])\n",
    "#and a data frame for plotting\n",
    "dfsummaryRatings<-as.data.frame(tableRatings)\n",
    "\n",
    "ggplot(dfsummaryRatings,aes(groupVar,Freq,fill=value)) +\n",
    "  geom_bar(stat=\"identity\") + \n",
    "  ylab(\"Count\") +\n",
    "  xlab(NULL) +\n",
    "  theme_bw(base_size =14)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#is it significant?\n",
    "chisq.test(tableRatings)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "longRatings<-rename(longRatings,c(\"value\"=\"Rating\"))\n",
    "#Finally\n",
    "#are all videos equally funny?\n",
    "videoRatings<-groupLongData(alldata,\"Rating Video1\", \"Rating Video2\", \"Rating Video3\")\n",
    "#let R know this is categorical data\n",
    "videoRatings$value<-factor(videoRatings$value)\n",
    "#lets collect responses into table (pick the appropriate columns)\n",
    "tableVideoRatings<-table(videoRatings[,4:5])\n",
    "#and a data frame for plotting\n",
    "dfsummaryRatings<-as.data.frame(tableVideoRatings)\n",
    "ggplot(dfsummaryRatings,aes(groupVar,Freq,fill=value)) +\n",
    "  geom_bar(stat=\"identity\") + \n",
    "  ylab(\"Count\") +\n",
    "  xlab(NULL) +\n",
    "  theme_bw(base_size =14)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#is it significant?\n",
    "chisq.test(tableVideoRatings)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A reviewer pointed out that chi squared test is overly conservative because our rating scale is nominal and we have a repeated measure with children giving a rating in 3 conditions. Therefore, after some research it seemed like a logistic regression was appropriate. We use the repolr package."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#if we are being strictly accurate then repolr package is best way to test ordinal rating data with repeated measures!\n",
    "longRatings$valueNum<-as.numeric(longRatings$value)\n",
    "longRatings$groupVarNum<-as.numeric(longRatings$groupVar)\n",
    "longRatings$ORDERNum<-as.numeric(longRatings$ORDER)\n",
    "funi<-repolr(valueNum~groupVarNum*ORDERNum,subjects = \"ID\", times=c(1,2,3),data=longRatings,categories =3)\n",
    "summary.repolr(funi)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Alternatively, we could to see if funniness is a predictor of amount of laughter in an one-way ANOVA."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "#another check look at amount of laughter by funniess rating\n",
    "laughsvsfunny<-cbind(longRatings,longLaughs$value)\n",
    "laughsvsfunny[\"groupVar\"]=NULL\n",
    "laughsvsfunny<-rename(laughsvsfunny,c(\"Rating\"= \"groupVar\", \"longLaughs$value\"=\"value\"))\n",
    "\n",
    "\n",
    "ezANOVA(laughsvsfunny, dv=value,wid =ID, between=groupVar)\n",
    "\n",
    "p<-plotOneWayANOVA(laughsvsfunny,varName=\"Laughs\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "R",
   "language": "R",
   "name": "ir"
  },
  "language_info": {
   "codemirror_mode": "r",
   "file_extension": ".r",
   "mimetype": "text/x-r-source",
   "name": "R",
   "pygments_lexer": "r",
   "version": "3.4.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}