This file demonstrates fitting a nonlinear discounting multilevel model to data when there are one or two variables that are assumed to predict the level of discounting. These examples do not test for interactions.

The Basics - Loading data, summarizing, plotting

library(nlme)
library(lattice)

Be sure that you either place the command and data files in your working directory or CHANGE your working directory (see R menus). To check the current working directory, use “getwd()”.

myd<-read.csv("Young E1 DelayDiscounting.dat")
summary(myd)
##     Subject        Delay        Indifference.point Sex      VGExperience  
##  01a    :  5   Min.   :   7.0   Min.   :  0        f:113   Min.   : 0.00  
##  01b    :  5   1st Qu.:  30.0   1st Qu.:402        m: 63   1st Qu.: 6.00  
##  01c    :  5   Median :  90.0   Median :812                Median :13.00  
##  02a    :  5   Mean   : 411.9   Mean   :679                Mean   :11.59  
##  02b    :  5   3rd Qu.: 365.0   3rd Qu.:989                3rd Qu.:17.00  
##  02d    :  5   Max.   :1825.0   Max.   :999                Max.   :21.00  
##  (Other):146                                               NA's   :10     
##     Cluster     
##  Min.   :1.000  
##  1st Qu.:1.000  
##  Median :1.000  
##  Mean   :1.561  
##  3rd Qu.:2.000  
##  Max.   :3.000  
##  NA's   :5

Note, if a categorical variable is designated with numbers, R will assume it’s continuous. If so, do something like this: myd$condition<-as.factor(myd$condition) before continuing. We’ll need to do this for Cluster in a later example.

For this data set, I have some very noisy data from a titration procedure. Variables included two between-subject predictors: sex (categorical) and video game experience (continuous). The data also includes missing values for VG experience and some missing delays.

I’ll start with an analysis involving none of the predictors to help me find good starting values.

mynls<-nls(Indifference.point~1000/(1+exp(logk)*Delay), data=myd, start=c(logk=-2))
summary(mynls) # suggests that -6 is a good starting value for logk
## 
## Formula: Indifference.point ~ 1000/(1 + exp(logk) * Delay)
## 
## Parameters:
##      Estimate Std. Error t value Pr(>|t|)    
## logk  -6.2735     0.1864  -33.65   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 357.1 on 175 degrees of freedom
## 
## Number of iterations to convergence: 14 
## Achieved convergence tolerance: 7.572e-06

If you want to explore various starting values of k, you may experience some unhelpful “solutions.”

Fit the multilevel hyperbolic model.

myr<-nlme(Indifference.point~1000/(1+exp(logk)*Delay), fixed=logk~1, random=logk~1|Subject, data=myd, start=c(logk=-6))
summary(myr)
## Nonlinear mixed-effects model fit by maximum likelihood
##   Model: Indifference.point ~ 1000/(1 + exp(logk) * Delay) 
##  Data: myd 
##       AIC      BIC    logLik
##   2359.89 2369.402 -1176.945
## 
## Random effects:
##  Formula: logk ~ 1 | Subject
##             logk Residual
## StdDev: 2.732835  139.972
## 
## Fixed effects: logk ~ 1 
##          Value Std.Error  DF   t-value p-value
## logk -6.159484 0.4664613 139 -13.20471       0
## 
## Standardized Within-Group Residuals:
##          Min           Q1          Med           Q3          Max 
## -7.019870117 -0.475215699 -0.009234761  0.259459548  2.835331276 
## 
## Number of Observations: 176
## Number of Groups: 37
hist(coef(myr)$logk)

The following highlights the messiness of the original data and the quality of the fits.

xyplot(fitted(myr)+Indifference.point~Delay|Subject, data=myd, type="a")