Overview

This script works through calculation of fear sensitivity scores using a Hyperbolic Tangent Growth Model. The objective is to obtain a measure of between-person differences in children’s fear sensitivity using data about children’s expression of fearful behavior during a laboratory visit wherein they completed six novelty tasks (e.g., novel stranger type tasks). We use the method described in
Buss, K., Davis, E., Ram, N., & Coccia, M. (in press). Dysregulated fear, social inhibition and respiratory sinus arrythmia: A replication and extension. Child Development.

The original implementation was done in SAS. Here, we work out the details of the R-based implementation, and extend the calculations to a new sample.

Outline

This script covers …

A. Describing the repeated measures data
B. Fitting the hyperbolic tangent growth model
C. Calculating fear sensitivity scores

Preliminaries

Loading libraries used in this script.

library(reshape2)
library(car)
library(psych)
library(ggplot2)
library(nlstools)
library(nlme)
library(lme4)

Reading in the data

#Loading the data
#read in the .csv file 
orig <- read.csv(file="lpadata_052917.csv",header=TRUE)

#replacing missing values (coded as -99) with NA
orig[orig == -99] <- NA

#fixing names
names(orig)[1] <- c("id")

#Looking at the top few rows of the wide data.
head(orig,6)
##   id     clown       ps       sa       sw    robot   spider
## 1  1  2.255639 20.39216  0.00000 26.66667 23.18681 14.47368
## 2  4 17.209302 37.31884  0.00000 31.33858 34.59459 48.72340
## 3  5 94.434783 74.78261 57.24138 39.04000 40.34783 21.90476
## 4  6  4.758065 43.75940       NA 43.69748 65.56701 12.26415
## 5  7  4.814815 69.84000 63.02326       NA 55.23810 31.57895
## 6  8 19.726962 20.30000  0.00000 29.28571 64.84848 27.40741

Reshaping the data

Reshape from wide to long for plotting etc.

#melting to long 
data_long <- melt(orig,
                 id.vars="id", 
                 variable.name="task",
                 na.rm=FALSE,
                 value.name = "fear")

#sorting for easy viewing
# order by id and task
data_long <- data_long[order(data_long$id,data_long$task), ]

Creating threat level variables (numeric variables)

#threat level ordering based on Kristin's theoretical ordering
data_long$threatlevel_theory <- recode(data_long$task, "'clown'=1; 'ps'=2; 'sw'=3; 'sa'=4; 'robot'=5; 'spider'=6", as.factor.result=FALSE)

#setting level of factor variable 
levels(data_long$task) <- c("clown","ps","sw","sa","robot","spider")

#threat level ordering based on empirical ordering (from original analysis)
data_long$threatlevel_empirical <- recode(data_long$task, "'clown'=3; 'ps'=4; 'sw'=1; 'sa'=2;  'robot'=5; 'spider'=6", as.factor.result=FALSE)

#creating new empirical order task variable  
data_long$taskempirical <- factor(data_long$task, 
                                  levels = c("sw","sa","clown","ps","robot","spider"))

#looking at the top few rows of the long data
head(data_long,18)
##      id   task      fear threatlevel_theory threatlevel_empirical
## 1     1  clown  2.255639                  1                     3
## 236   1     ps 20.392157                  2                     4
## 471   1     sw  0.000000                  4                     1
## 706   1     sa 26.666667                  3                     2
## 941   1  robot 23.186813                  5                     5
## 1176  1 spider 14.473684                  6                     6
## 2     4  clown 17.209302                  1                     3
## 237   4     ps 37.318841                  2                     4
## 472   4     sw  0.000000                  4                     1
## 707   4     sa 31.338583                  3                     2
## 942   4  robot 34.594595                  5                     5
## 1177  4 spider 48.723404                  6                     6
## 3     5  clown 94.434783                  1                     3
## 238   5     ps 74.782609                  2                     4
## 473   5     sw 57.241379                  4                     1
## 708   5     sa 39.040000                  3                     2
## 943   5  robot 40.347826                  5                     5
## 1178  5 spider 21.904762                  6                     6
##      taskempirical
## 1            clown
## 236             ps
## 471             sw
## 706             sa
## 941          robot
## 1176        spider
## 2            clown
## 237             ps
## 472             sw
## 707             sa
## 942          robot
## 1177        spider
## 3            clown
## 238             ps
## 473             sw
## 708             sa
## 943          robot
## 1178        spider
#reshaping from long to wide
data_wide <- dcast(data_long, id ~ task, value.var="fear")

A. Describing the data

Descriptives

Basic descriptives of the 6-task data.

#chekcing data frame formats
str(data_long)
## 'data.frame':    1410 obs. of  6 variables:
##  $ id                   : int  1 1 1 1 1 1 4 4 4 4 ...
##  $ task                 : Factor w/ 6 levels "clown","ps","sw",..: 1 2 3 4 5 6 1 2 3 4 ...
##  $ fear                 : num  2.26 20.39 0 26.67 23.19 ...
##  $ threatlevel_theory   : num  1 2 4 3 5 6 1 2 4 3 ...
##  $ threatlevel_empirical: num  3 4 1 2 5 6 3 4 1 2 ...
##  $ taskempirical        : Factor w/ 6 levels "sw","sa","clown",..: 3 4 1 2 5 6 3 4 1 2 ...
str(data_wide)
## 'data.frame':    235 obs. of  7 variables:
##  $ id    : int  1 4 5 6 7 8 9 10 11 12 ...
##  $ clown : num  2.26 17.21 94.43 4.76 4.81 ...
##  $ ps    : num  20.4 37.3 74.8 43.8 69.8 ...
##  $ sw    : num  0 0 57.2 NA 63 ...
##  $ sa    : num  26.7 31.3 39 43.7 NA ...
##  $ robot : num  23.2 34.6 40.3 65.6 55.2 ...
##  $ spider: num  14.5 48.7 21.9 12.3 31.6 ...
#means, sd, etc.
describe(data_wide)
##        vars   n    mean      sd  median trimmed     mad  min     max
## id        1 235 3330.85 3101.59 6035.00 3339.56 1031.89 1.00 6820.00
## clown     2 235   25.90   23.29   22.12   23.07   27.50 0.22   98.86
## ps        3 234   37.32   21.71   37.36   36.45   23.89 0.09   99.79
## sw        4 229   27.67   21.06   26.12   26.00   23.72 0.00   97.47
## sa        5 230   25.21   17.53   26.43   24.06   16.43 0.00   80.68
## robot     6 233   51.27   25.92   51.80   51.96   29.43 0.00   99.68
## spider    7 233   51.05   22.39   50.75   51.52   26.03 0.00   99.69
##          range  skew kurtosis     se
## id     6819.00 -0.10    -1.99 202.33
## clown    98.64  0.84    -0.11   1.52
## ps       99.70  0.40    -0.35   1.42
## sw       97.47  0.59    -0.23   1.39
## sa       80.68  0.53     0.33   1.16
## robot    99.68 -0.18    -0.85   1.70
## spider   99.69 -0.14    -0.74   1.47
#boxplot by task
ggplot(data=data_long[which(data_long$id > 0), ], aes(x=task, y=fear)) + 
  geom_boxplot(notch = TRUE) +
  stat_summary(fun.y="mean", geom="point", shape=23, size=3, fill="white") +
  labs(x = "Task", y = "Fear")
## Warning: Removed 16 rows containing non-finite values (stat_boxplot).
## Warning: Removed 16 rows containing non-finite values (stat_summary).

#boxplot by task empirical
ggplot(data=data_long[which(data_long$id > 0), ], aes(x=taskempirical, y=fear)) + 
  geom_boxplot(notch = TRUE) +
  stat_summary(fun.y="mean", geom="point", shape=23, size=3, fill="white") +
  labs(x = "Task", y = "Fear")
## Warning: Removed 16 rows containing non-finite values (stat_boxplot).

## Warning: Removed 16 rows containing non-finite values (stat_summary).

#correlations
# Correlations 
round(cor(data_wide[,c("clown","ps","sa","sw","robot","spider")], 
          use="complete.obs",method="spearman"),2)
##        clown   ps   sa   sw robot spider
## clown   1.00 0.29 0.34 0.30  0.28   0.37
## ps      0.29 1.00 0.25 0.15  0.12   0.07
## sa      0.34 0.25 1.00 0.19  0.16   0.08
## sw      0.30 0.15 0.19 1.00  0.17   0.11
## robot   0.28 0.12 0.16 0.17  1.00   0.52
## spider  0.37 0.07 0.08 0.11  0.52   1.00
#pairs plot from the psych library
pairs.panels(data_wide[,c("clown","ps","sa","sw","robot","spider"),])