Modeling non-linear change

Modeling non-linear change over time

• Choosing a functional form
  • Function must be adequate to data
  • Dynamic consistency
  • Prediction
• Higher-order polynomials
• Orthogonal polynomials

Function must be adequate to data

• Easy plotting residuals vs. fitted: ggplot(fortify(m), aes(.fitted, .resid)) + geom_point()

Two kinds of non-linearity

1. Non-linear in variables (Time): \(Y_{ij} = \beta_{0i} + \beta_{1i} \cdot Time_{j} + \beta_{2i} \cdot Time^2_{j} + \epsilon_{ij}\)
2. Non-linear in parameters (b): \(Y = \frac{p-b}{1+exp(4 \cdot \frac{s}{p-b} \cdot (c-t))} + b\)

• Dynamic Consistency: The model of the average is equal to the average of the individual models
  • Recall: random effects must have a mean of 0
  • Average of individual deviations is 0
  • So, the average of the individual models will have 0 deviation from the model of the average
  • Not true of some tempting functional forms

Example of lack of dynamic consistency

Logistic power peak function (Scheepers, Keller, & Lapata, 2008) fit to semantic competition data (Mirman & Magnuson, 2009).

• Standard statistical inference assumes that the group average represents central tendency of individuals
• For dynamically inconsistent models this is not true
• Can’t make inferences about central tendencies without dynamic consistency

Prediction: Two kinds

Using higher-order polynomials

• Can model any curve shape
• Dynamically consistent
• Bad at capturing asymptotic behavior
  • Don’t extrapolate
  • Try to avoid long flat sections

• How to choose polynomial order?
• Curve shape
• Statistical: include only terms that statistically improve model fit
• Theoretical: include only terms that are predicted to matter

Natural vs. Orthogonal polynomials

• Natural Polynomials: Correlated time terms
• Orthogonal Polynomials: Uncorrelated time terms
  • Need to specify range and order

Interpreting orthogonal polynomial terms

Intercept (\(\beta_0\)): Overall average

Interpreting orthogonal polynomial terms

Intercept (\(\beta_0\)): Overall average

Linear (\(\beta_1\)): Overall slope

Quadratic (\(\beta_2\)): Centered rise and fall rate

Cubic, Quartic, … (\(\beta_3, \beta_4, ...\)): Inflection steepness

Example

Effect of transitional probability on word learning

summary(WordLearnEx)

##     Subject       TP          Block         Accuracy    
##  244    : 10   Low :280   Min.   : 1.0   Min.   :0.000  
##  253    : 10   High:280   1st Qu.: 3.0   1st Qu.:0.667  
##  302    : 10              Median : 5.5   Median :0.833  
##  303    : 10              Mean   : 5.5   Mean   :0.805  
##  305    : 10              3rd Qu.: 8.0   3rd Qu.:1.000  
##  306    : 10              Max.   :10.0   Max.   :1.000  
##  (Other):500

ggplot(WordLearnEx, aes(Block, Accuracy, color=TP)) + 
  stat_summary(fun.data=mean_se, geom="pointrange") + 
  stat_summary(fun.y=mean, geom="line")

Example: Prepare for GCA

Create 2nd-order orthogonal polynomial

WordLearn.gca <- code_poly(WordLearnEx, predictor="Block", poly.order=2, orthogonal=TRUE)

summary(WordLearn.gca)

##     Subject       TP          Block         Accuracy      Block.Index  
##  244    : 10   Low :280   Min.   : 1.0   Min.   :0.000   Min.   : 1.0  
##  253    : 10   High:280   1st Qu.: 3.0   1st Qu.:0.667   1st Qu.: 3.0  
##  302    : 10              Median : 5.5   Median :0.833   Median : 5.5  
##  303    : 10              Mean   : 5.5   Mean   :0.805   Mean   : 5.5  
##  305    : 10              3rd Qu.: 8.0   3rd Qu.:1.000   3rd Qu.: 8.0  
##  306    : 10              Max.   :10.0   Max.   :1.000   Max.   :10.0  
##  (Other):500                                                           
##      poly1            poly2       
##  Min.   :-0.495   Min.   :-0.348  
##  1st Qu.:-0.275   1st Qu.:-0.261  
##  Median : 0.000   Median :-0.087  
##  Mean   : 0.000   Mean   : 0.000  
##  3rd Qu.: 0.275   3rd Qu.: 0.174  
##  Max.   : 0.495   Max.   : 0.522  
## 

Example: Fit the models

library(lme4)
#fit base model
m.base <- lmer(Accuracy ~ (poly1+poly2) + (poly1+poly2 | Subject), 
               data=WordLearn.gca, REML=F)
#add effect of TP on intercept 
m.0 <- lmer(Accuracy ~ (poly1+poly2) + TP + (poly1+poly2 | Subject), 
            data=WordLearn.gca, REML=F)
#add effect on slope
m.1 <- lmer(Accuracy ~ (poly1+poly2) + TP + TP:poly1 + (poly1+poly2 | Subject), 
            data=WordLearn.gca, REML=F)
#add effect on quadratic
m.2 <- lmer(Accuracy ~ (poly1+poly2)*TP + (poly1+poly2 | Subject), 
            data=WordLearn.gca, REML=F)

Example: Compare the model fits

anova(m.base, m.0, m.1, m.2)

## Data: WordLearn.gca
## Models:
## m.base: Accuracy ~ (poly1 + poly2) + (poly1 + poly2 | Subject)
## m.0: Accuracy ~ (poly1 + poly2) + TP + (poly1 + poly2 | Subject)
## m.1: Accuracy ~ (poly1 + poly2) + TP + TP:poly1 + (poly1 + poly2 | 
## m.1:     Subject)
## m.2: Accuracy ~ (poly1 + poly2) * TP + (poly1 + poly2 | Subject)
##        Df  AIC  BIC logLik deviance Chisq Chi Df Pr(>Chisq)  
## m.base 10 -331 -288    175     -351                          
## m.0    11 -330 -283    176     -352  1.55      1      0.213  
## m.1    12 -329 -277    176     -353  0.36      1      0.550  
## m.2    13 -333 -276    179     -359  5.95      1      0.015 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Example: Parameter estimates with p-values

get_pvalues(m.2)

##               Estimate Std..Error  t.value  p.normal p.normal.star
## (Intercept)   0.778525    0.02173 35.83141 0.000e+00           ***
## poly1         0.286315    0.03779  7.57682 3.531e-14           ***
## poly2        -0.050849    0.03319 -1.53217 1.255e-01              
## TPHigh        0.052961    0.03073  1.72357 8.478e-02             .
## poly1:TPHigh  0.001075    0.05344  0.02012 9.839e-01              
## poly2:TPHigh -0.116455    0.04693 -2.48121 1.309e-02             *

Example: Plot the model fit

ggplot(WordLearn.gca, aes(Block, Accuracy, color=TP)) + 
  stat_summary(fun.data=mean_se, geom="pointrange") + 
  stat_summary(aes(y=fitted(m.2)), fun.y=mean, geom="line") +
  theme_bw() + expand_limits(y=c(0.5, 1)) + scale_x_continuous(breaks = 1:10)

Questions?

Random effects

• Keep it maximal: A full or maximal random effect structure is when all of the factors that could hypothetically vary across individual observational units are allowed to do so.
• Incomplete random effects can inflate false alarms
• Full random effects can produce convergence problems (we’ll come back to this issue in a little while)

Within-subject effects

Example: Target fixation in spoken word-to-picure matching (VWP)

ggplot(TargetFix, aes(Time, meanFix, color=Condition, fill=Condition)) +
  stat_summary(fun.y=mean, geom="line") +
  stat_summary(fun.data=mean_se, geom="ribbon", color=NA, alpha=0.3) +
  theme_bw() + expand_limits(y=c(0,1)) + legend_positioning(c(1,1)) +
  labs(y="Fixation Proportion", x="Time since word onset (ms)")

Within-subject effects: fit full GCA model

m.full <- lmer(meanFix ~ (poly1+poly2+poly3)*Condition + #fixed effects
                 (poly1+poly2+poly3 | Subject) + #random effects of Subject
                 (poly1+poly2+poly3 | Subject:Condition), #random effects of Subj by Cond
               data=TargetFix.gca, REML=F)

## singular fit

get_pvalues(m.full)

##                      Estimate Std..Error   t.value  p.normal p.normal.star
## (Intercept)         0.4773228    0.01385 34.457777 0.000e+00           ***
## poly1               0.6385604    0.05994 10.654180 0.000e+00           ***
## poly2              -0.1095979    0.03849 -2.847573 4.405e-03            **
## poly3              -0.0932612    0.02330 -4.002200 6.276e-05           ***
## ConditionLow       -0.0581122    0.01879 -3.093227 1.980e-03            **
## poly1:ConditionLow  0.0003188    0.06579  0.004846 9.961e-01              
## poly2:ConditionLow  0.1635455    0.05393  3.032545 2.425e-03            **
## poly3:ConditionLow -0.0020869    0.02704 -0.077168 9.385e-01

Within-subject effects: random effects

# str(ranef(m.full))
VarCorr(m.full)

##  Groups            Name        Std.Dev. Corr             
##  Subject:Condition (Intercept) 0.0405                    
##                    poly1       0.1404   -0.43            
##                    poly2       0.1124   -0.33  0.72      
##                    poly3       0.0417    0.13 -0.49 -0.43
##  Subject           (Intercept) 0.0124                    
##                    poly1       0.1195    0.91            
##                    poly2       0.0165   -0.42 -0.76      
##                    poly3       0.0421   -0.85 -0.99  0.83
##  Residual                      0.0438

What is being estimated?

• Random variance and covariance
• Unit-level random effects (but constrained to have mean = 0)

This is why df for parameter estimates are poorly defined in MLR

Alternative random effects structure

m.left <- lmer(meanFix ~ (poly1+poly2+poly3)*Condition + #fixed effects
                ((poly1+poly2+poly3)*Condition | Subject), #random effects
              data=TargetFix.gca, REML=F)

## Warning in commonArgs(par, fn, control, environment()): maxfun < 10 *
## length(par)^2 is not recommended.

## singular fit

get_pvalues(m.left)

##                      Estimate Std..Error   t.value  p.normal p.normal.star
## (Intercept)         0.4773228    0.01679 28.434031 0.000e+00           ***
## poly1               0.6385604    0.05146 12.409015 0.000e+00           ***
## poly2              -0.1095979    0.03717 -2.948228 3.196e-03            **
## poly3              -0.0932612    0.02062 -4.523841 6.073e-06           ***
## ConditionLow       -0.0581122    0.02110 -2.753958 5.888e-03            **
## poly1:ConditionLow  0.0003188    0.07477  0.004264 9.966e-01              
## poly2:ConditionLow  0.1635455    0.06274  2.606891 9.137e-03            **
## poly3:ConditionLow -0.0020869    0.03334 -0.062598 9.501e-01

Alternative random effects structure requires more parameters

Convergence problems

Warning message: 
In (function (fn, par, lower = rep.int(-Inf, n), upper = rep.int(Inf, : 
    failure to converge in 10000 evaluations

Need to simplify random effects

• Remove random effects of higher-order terms

  Outcome ~ (poly1+poly2+poly3)*Condition + (poly1+poly2+poly3 | Subject)
  Outcome ~ (poly1+poly2+poly3)*Condition + (poly1+poly2 | Subject)

• Remove correlation between random effects

  Outcome ~ (poly1+poly2+poly3)*Condition + (1 | Subject) + 
    (0+poly1 | Subject) + (0+poly2 | Subject) + (0+poly3 | Subject)

Participants as fixed vs. random effects

In general, participants should be treated as random effects.

This captures the typical assumption of random sampling from some population to which we wish to generalize.

Pooling

http://tjmahr.github.io/plotting-partial-pooling-in-mixed-effects-models/

Participants as fixed vs. random effects

In general, participants should be treated as random effects.

This captures the typical assumption of random sampling from some population to which we wish to generalize.

Shrinkage

http://tjmahr.github.io/plotting-partial-pooling-in-mixed-effects-models/

Participants as fixed vs. random effects

In general, participants should be treated as random effects.

This captures the typical assumption of random sampling from some population to which we wish to generalize.

Exceptions are possible: e.g., neurological/neuropsychological case studies where the goal is to characterize the pattern of performance for each participant, not generalize to a population.

Lunch break: noon-1pm

• Try GCA on your own data
• Ask questions