Week 8, Lec 3: Issues in Mixed Models Mixed models issues raised in - - PowerPoint PPT Presentation

Week 8, Lec 3: Issues in Mixed Models

• Mixed models issues raised in class.
• HW-5 with ‘kv0.csv’
• A Generalised LMM (GLMM) analysis and some

complications, e.g., corr(int, slope) = ±1.

• More on p-values.
• An example with nested factors
• How to implement and test/compare models in

R? Use of ML vs REML deviance.

1

Robert: ¡AIC ¡is ¡a ¡principled ¡index!

• Suppose ¡the ¡true ¡likelihood ¡of ¡the ¡data, ¡y, ¡is ¡g(y), ¡and ¡
• ur ¡model ¡for ¡y ¡contains ¡k ¡parameters, ¡θ, ¡and ¡yields ¡

the ¡(modeled) ¡likelihood ¡of ¡the ¡data, ¡y, ¡is ¡f(y, ¡θ). ¡ ¡We ¡ would ¡want ¡our ¡model ¡and ¡our ¡ML ¡es=mate, ¡θ*, ¡of ¡θ ¡ to ¡be ¡such ¡that ¡the ¡‘distance’ ¡between ¡g(y) ¡and ¡f(y, ¡ θ*) ¡to ¡be ¡a ¡minimum. ¡ ¡Because ¡we ¡do ¡not ¡know ¡g(y), ¡ we ¡can’t ¡calculate ¡this ¡‘distance’, ¡but ¡we ¡can ¡try ¡to ¡ approximate ¡it ¡in ¡some ¡way. ¡

• Akaike ¡introduced ¡a ¡‘distance’ ¡index ¡based ¡on ¡

informa=on ¡theory, ¡the ¡Kullback-­‑Leibler ¡(KL) ¡index. ¡ ¡

2 ¡

2I(θ) = E log g(y) f (y,θ) ⎛ ⎝ ⎜ ⎞ ⎠ ⎟ .

Robert: ¡AIC ¡is ¡a ¡principled ¡index!

• Akaike ¡showed ¡that, ¡apart ¡from ¡a ¡constant ¡term, ¡(i) ¡the ¡

deviance ¡of ¡a ¡model ¡(previously ¡defined) ¡is ¡a ¡biased ¡ es=mate ¡of ¡2I(θ), ¡and ¡(ii) ¡the ¡bias ¡is ¡approximately ¡ equal ¡to ¡twice ¡the ¡dimensionality, ¡k, ¡of ¡θ. ¡ ¡Hence, ¡an ¡ approx ¡unbiased ¡es=mate ¡of ¡the ¡K-­‑L ¡distance ¡is: ¡

¡AIC ¡= ¡Deviance ¡+ ¡2k ¡= ¡Deviance ¡+ ¡2(p ¡+ ¡1), ¡

• where ¡p ¡is ¡the ¡number ¡of ¡explanatory ¡variables ¡in ¡a ¡

regression, ¡and ¡the ¡‘1’ ¡refers ¡to ¡the ¡intercept. ¡

• AIC ¡is ¡well-­‑behaved ¡when ¡n ¡>> ¡p. ¡If, ¡say, ¡p ¡= ¡n/2 ¡or ¡

larger, ¡AIC ¡could ¡decrease ¡as ¡p ¡increases ¡– ¡giving ¡an ¡ unfair ¡advantage ¡to ¡complex ¡models! ¡

3 ¡

2I(θ) = E log g(y) f (y,θ) ⎛ ⎝ ⎜ ⎞ ⎠ ⎟ .

Use ¡of ¡AIC ¡

• AIC ¡can ¡be ¡used ¡to ¡assess ¡the ¡rela=ve ¡fit ¡of ¡non-­‑nested ¡

models ¡– ¡the ¡smaller ¡AIC ¡is, ¡the ¡beSer. ¡ ¡(RE)ML ¡ deviance ¡is ¡used ¡to ¡test ¡nested ¡models; ¡the ¡difference ¡ in ¡deviance ¡between ¡models ¡has ¡the ¡chi-­‑square ¡ distrn, ¡and ¡this ¡is ¡the ¡basis ¡of ¡the ¡test. ¡

• Suppose ¡model, ¡M*, ¡has ¡the ¡minimum ¡AIC. ¡ ¡Is ¡another ¡

model, ¡M, ¡‘as ¡good ¡as’ ¡M*? ¡ ¡Let ¡d ¡= ¡AICM ¡– ¡AICmin. ¡ ¡A ¡ rule-­‑of-­‑thumb ¡is: ¡ ¡

– ‘Yes’, ¡if ¡0 ¡< ¡d ¡< ¡2; ¡ ¡ – ‘Probably ¡not’, ¡if ¡4 ¡< ¡d ¡< ¡7; ¡and ¡ ¡ – ‘No’, ¡if ¡10 ¡< ¡d. ¡

• Source: ¡hSp://myweb.uiowa.edu/cavaaugh/ms_lec_2_ho.pdf ¡

Ian: ¡Is ¡var(Int) ¡changed ¡by ¡rescaling ¡the ¡IV?

• Consider ¡then ¡intercept-­‑only ¡model ¡with ¡fixed ¡effects ¡

factor, ¡T, ¡and ¡random ¡effect, ¡S, ¡ ¡

Xijk = µ + ai + b*Tj + eijk.

• The ¡random ¡effects ¡variance, ¡var(Int) ¡= ¡σa

2, ¡denotes ¡a ¡

substan=ve ¡quan=ty, ¡namely, ¡varia=on ¡across ¡Ss ¡in ¡their ¡ ‘average’ ¡response, ¡that ¡cannot ¡depend ¡on ¡which ¡level ¡of ¡ T ¡is ¡chosen ¡as ¡the ¡origin. ¡(However, ¡rescaling ¡can ¡affect ¡ the ¡correlaKon ¡between ¡2 ¡random ¡effects.) ¡

d1$timec = d1$time - 2; d1$weightc = scale(d1$weight, scale = F) res3a = lmer(weightc ~ timec + (1 | id), d1)

• The ¡results ¡are ¡the ¡same ¡as ¡with ¡the ¡uncentered ¡model. ¡

5 ¡

Yuan ¡Chang ¡& ¡Monica: ¡Can ¡adding ¡a ¡‘random ¡ slope’ ¡for ¡X ¡remove ¡the ¡fixed ¡effect ¡of ¡X?

• Consider ¡the ¡‘intercept-­‑only’ ¡and ¡‘intercept+slope’ ¡

models: ¡

res3 = lmer(weight ~ time + (1 | id), d1) res5 = lmer(weight ~ time + (1 + time | id), d1)

• Assuming ¡a ¡random ¡slope ¡for ¡>me ¡is ¡the ¡same ¡as ¡

assuming ¡a ¡S*>me ¡interac=on. ¡ ¡If ¡there ¡is ¡a ¡significant ¡ interac=on, ¡the ¡fixed ¡effect ¡of ¡>me ¡may ¡or ¡may ¡not ¡be ¡

• significant. ¡

6 ¡

Zeynep: ¡Isn’t ¡the ¡Compound ¡Symmetry ¡(CS) ¡assn. ¡ violated ¡in ¡longitudinal ¡designs?

• The ¡CS ¡assn ¡implies ¡that ¡the ¡correla=on ¡across ¡Ss ¡

between ¡the ¡scores ¡at ¡T ¡= ¡1, ¡2, ¡… ¡sa=sfy: ¡cor(x1i, ¡x2i) ¡= ¡ cor(x1i, ¡x3i) ¡= ¡cor(x3i, ¡x2i) ¡= ¡…. ¡

• In ¡a ¡longitudinal ¡design ¡with ¡T ¡= ¡>me, ¡it ¡is ¡likely ¡that ¡

successive ¡observa=ons ¡will ¡be ¡more ¡highly ¡correlated ¡ than ¡observa=ons ¡separated ¡by ¡many ¡=me-­‑points. ¡ ¡ Ojen ¡a ¡reasonable ¡alterna=ve ¡is ¡the ¡assump=on ¡that: ¡

¡cor(xt,i, ¡xt+k,i) ¡= ¡ρk, ¡ρ ¡< ¡1, ¡

• and ¡depends ¡only ¡on ¡the ¡‘distance’, ¡k, ¡between ¡
• bserva=ons. ¡ ¡This ¡model ¡is ¡called ¡an ¡autoregressive ¡

model ¡of ¡order ¡1, ¡AR(1). ¡ ¡It ¡applies ¡to ¡mixed ¡models ¡and ¡ to ¡univariate ¡Kme-­‑series. ¡ ¡(See ¡Russ ¡Poldrack’s ¡lecture ¡

• n ¡11/19/14.) ¡

7 ¡

8

CS and other assumptions about the variance- covariance structure in Mixed Models

• In a repeated measures design or in a random effects

design, let cij be the covariance between two

• bservations, Xi and Xj, from the same subject.
• If there are only within-subject factors, then a

correlation among repeated measures (which is a violation of the “I” in the LINE model) does not change the definition of the F ratios for testing the effects of the within-subjects factors.

• If there is a between-subjects factor, we need to use

the Mixed Models procedures.

9

• With this notation, cii is the variance of Xi. We now show

a few of the possible covariance matrices, C = (cij), that may apply to your study. For simplicity, we consider the case, i, j = 1, 2, 3.

• Compound symmetry: This is appropriate when the

errors in Xi and Xj are correlated, and the correlation is the same for all pairs of different responses. This condition is usually met when we randomize the order in which the subject is exposed to the different within-subject factor levels.

2 1

1 1 , | | 1. 1 C ρ ρ σ ρ ρ ρ ρ ρ ⎛ ⎞ ⎜ ⎟ = < ⎜ ⎟ ⎜ ⎟ ⎝ ⎠

10

• Autoregressive structure, AR(1): This is appropriate

when the subject is exposed to the different within-subject factor levels in a fixed order. In that case, the errors in Xi and Xj are more highly correlated if the two responses are adjacent than if they are farther apart.

2 2 2 2

1 1 , | | 1. 1 C ρ ρ σ ρ ρ ρ ρ ρ ⎛ ⎞ ⎜ ⎟ = < ⎜ ⎟ ⎜ ⎟ ⎝ ⎠

11

• Unstructured covariance matrix: This most general

model can be tested against special cases, e.g., compound symmetry, by using the goodness-of-fit index, (-2)*log likelihood, that is given in lme() and lmer()

• utput. SPSS, but not lmer(), makes it easy to specify

covariance structures other than CS!

2 1 12 13 2 3 12 2 23 2 13 23 3

. C σ σ σ σ σ σ σ σ σ ⎛ ⎞ ⎜ ⎟ = ⎜ ⎟ ⎜ ⎟ ⎝ ⎠

12

HW-5 using ‘kv0.csv’

• The fictitious data set, ‘kv0.csv’, is based on
• Dr. Katerina Velanova’s FYP data that she

shared with us when she was a Psych 252 student.

• I found this acknowledgement in an old

handout: Data and analyses courtesy of Katerina Velanova, to whom Psy 252 owes a debt of gratitude [11/8/96].

• Relevant R script: ‘skv1.r’

13

• Ten subjects were run in a divided attention

condition (attn = 1), and 10 different subjects were run in a focused attention condition (attn = 2). Each subject was then tested on a word task (anagrams) that had a unique solution (numsol = 1), two solutions (numsol = 2), or multiple solutions (numsol = 3). The dependent measure was a memory score (higher numbers reflect better performance).

• This design involves both between- and within-

factors and is probably the most useful design for psychologists.

14

Attention SubjID Numsol=1 Numsol=2 Numsol=3 1.00 2.00 4.00 7.00 2.00 3.00 4.00 5.00 3.00 3.00 5.00 6.00 4.00 5.00 7.00 5.00 Attn = 1 5.00 4.00 5.00 8.00 6.00 5.00 5.00 6.00 7.00 5.00 4.50 6.00 8.00 5.00 7.00 8.00 9.00 2.00 3.00 7.00 10.00 6.00 5.00 6.00 1.00 6.00 5.00 6.00 2.00 8.00 9.00 8.00 3.00 6.00 5.00 9.00 4.00 8.00 8.00 7.00 Attn = 2 5.00 8.00 8.00 7.00 6.00 6.00 8.00 7.00 7.00 7.00 7.00 6.00 8.00 7.00 8.00 6.00 9.00 5.00 6.00 6.00 10.00 6.00 6.00 5.00

15

16

#Scripts to analyse kv0.csv; in ‘skv1.r’ source('makerm.r') library(lme4) library(nlme) d0 = read.csv('kv0.csv', header = TRUE) con1 = cbind(c(-2,1,1), c(0,-1,1)) #Convert d0 into long form with make.rm() #Now we know how to use melt() {ggplot2}! Now we know how to use melt() {ggplot2}!

dl1 = make.rm(constant = 1:2, repeated = 3:5, data = d0) colnames(dl1) = c('suid','atten','score','nsol')

17

dl1$nsol = factor(dl1$nsol, labels=c('one','two','many')) dl1$atten = factor(dl1$atten, labels=c('div','foc')) dl1$suid = factor(dl1$suid) contrasts(dl1$atten) = c(-1,1) contrasts(dl1$nsol) = cbind(lin=c(-1,0,1), quad=c(-1,2,-1))

18

rs.aov = aov(score ~ atten*nsol + Error(suid), dl1) print(summary(rs.aov))

Error: suid Df Sum Sq Mean Sq F value Pr(>F) atten 1 42.504 42.504 17.598 0.0005442 *** Residuals 18 43.475 2.415 Error: Within Df Sum Sq Mean Sq F value Pr(>F) nsol 2 14.408 7.2042 6.5909 0.003639 ** atten:nsol 2 15.408 7.7042 7.0483 0.002612 ** Residuals 36 39.350 1.0931

19

rs.lme1 = lmer(score ~ atten*nsol + (1|suid), dl1) print(summary(rs.lme1))

Linear mixed model fit by REML Formula: score ~ atten * nsol + (1 | suid) Data: dl1 AIC BIC logLik deviance REMLdev 213.5 230.2 -98.73 185.1 197.5 Random effects: Groups Name Variance Std.Dev. suid (Intercept) 0.44074 0.66388 Residual 1.09306 1.04549 Number of obs: 60, groups: suid, 20 Fixed effects: Estimate Std. Error t value (Intercept) 5.958333 0.200628 29.698 atten1 0.841667 0.200628 4.195 nsollin 0.600000 0.165307 3.630 nsolquad 0.008333 0.095440 0.087 atten1:nsollin -0.600000 0.165307 -3.630 atten1:nsolquad 0.091667 0.095440 0.960

20

rs.lme2a = lmer(score ~ atten * nsol + (1 + nsol|suid), dl1) print(summary(rs.lme2a))

Linear mixed model fit by REML Formula: score ~ atten * nsol + (1 + nsol | suid) Data: dl1 AIC BIC logLik deviance REMLdev 215.5 242.8 -94.77 176.4 189.5 Random effects: Groups Name Variance Std.Dev. Corr suid (Intercept) 0.663354 0.81447 nsollin 0.459614 0.67795 -0.401 nsolquad 0.069408 0.26345 0.494 -0.804 Residual 0.425219 0.65209 Number of obs: 60, groups: suid, 20

21

Fixed effects: Estimate Std. Error t value (Intercept) 5.958333 0.200636 29.697 atten1 0.841667 0.200636 4.195 nsollin 0.600000 0.183334 3.273 nsolquad 0.008333 0.083751 0.100 atten1:nsollin -0.600000 0.183334 -3.273 atten1:nsolquad 0.091667 0.083751 1.095

Describe the data!

22 ¡

Judd, ¡C.M., ¡Weslall, ¡J., ¡& ¡Kenny, ¡D.A. ¡(2012). ¡Trea=ng ¡s=muli ¡as ¡a ¡random ¡factor ¡ in ¡social ¡psychology: ¡A ¡new ¡and ¡comprehensive ¡solu=on ¡to ¡a ¡pervasive ¡but ¡largely ¡ ignored ¡problem. ¡Journal ¡of ¡Personality ¡and ¡Social ¡Psychology ¡ ¡ Baayen, ¡R.H., ¡Davidson, ¡D.J., ¡& ¡Bates, ¡D.M. ¡(2007). ¡Mixed-­‑effects ¡modeling ¡with ¡ crossed ¡random ¡effects ¡for ¡subjects ¡and ¡items. ¡ ¡Journal ¡of ¡Memory ¡and ¡

• Language. ¡

¡ Baayen, ¡R.H. ¡(2008) ¡Analyzing ¡linguis>c ¡data: ¡A ¡prac>cal ¡introduc>on ¡to ¡sta>s>cs. ¡ [Chapter ¡7 ¡on ¡mixed ¡models]. ¡PDF ¡available ¡online ¡at ¡: ¡hSp://www.ualberta.ca/ ~baayen/publica=ons/baayenCUPstats.pdf ¡ ¡ Barr, ¡D.J., ¡Levy, ¡R., ¡Scheepers, ¡C., ¡& ¡Tily, ¡H.J. ¡(2013). ¡Random ¡effects ¡structure ¡for ¡ confirmatory ¡hypothesis ¡tes=ng: ¡Keep ¡it ¡maximal. ¡Journal ¡of ¡Memory ¡and ¡ Language, ¡68, ¡255-­‑278. ¡ ¡ Grimm, ¡K.J., ¡Ram, ¡N., ¡& ¡Hamagami, ¡F. ¡(2011). ¡ ¡Nonlinear ¡growth ¡curves ¡in ¡ developmental ¡research. ¡ ¡Child ¡Development, ¡82, ¡1357-­‑1371. ¡ ¡ ¡

Further ¡Readings ¡

GLMM: ¡Binary ¡logisKc ¡regression ¡

• In ¡the ¡‘study’ ¡phase ¡of ¡a ¡memory ¡expt, ¡S ¡learns ¡many ¡
• words. ¡ ¡In ¡the ¡‘test’ ¡phase, ¡S ¡is ¡presented ¡with ¡a ¡word ¡
• n ¡each ¡trial, ¡and ¡S’s ¡task ¡is ¡to ¡respond ¡‘old’, ¡if ¡the ¡

presented ¡word ¡was ¡previously ¡studied; ¡and ¡‘new’, ¡

• therwise. ¡ ¡The ¡2 ¡categories ¡of ¡presented ¡sKmuli ¡are ¡

OLD ¡(i.e., ¡previously ¡studied) ¡and ¡NEW ¡(i.e., ¡new). ¡

• In ¡the ¡Theory ¡of ¡Signal ¡Detec=on ¡(TSD), ¡the ¡data ¡are ¡
• jen ¡summarised ¡by ¡2 ¡condi=onal ¡probabili=es: ¡pH ¡= ¡

P(‘old’ ¡| ¡OLD), ¡the ¡‘hit ¡rate’; ¡and ¡pFA ¡= ¡P(‘old’ ¡| ¡NEW), ¡ the ¡‘false ¡alarm ¡rate’. ¡ ¡A ¡S’s ¡ability ¡or ¡capacity ¡or ¡ sensi>vity ¡is ¡usually ¡defined ¡as ¡a ¡‘generalised ¡ difference’ ¡between ¡the ¡2 ¡rates, ¡e.g., ¡d0 ¡= ¡pH ¡– ¡pFA, ¡d’ ¡ = ¡zH ¡– ¡zFA, ¡or ¡d* ¡= ¡logit(pH) ¡– ¡logit(pFA). ¡ ¡The ¡last, ¡d*, ¡is ¡ easily ¡obtained ¡from ¡fixng ¡a ¡logisKc ¡regression. ¡

23 ¡

Binary ¡logisKc ¡regression ¡

• Suppose ¡s=muli ¡vary ¡within ¡each ¡category ¡(e.g., ¡

words ¡or ¡faces ¡differ ¡in ¡how ¡memorable ¡they ¡are, ¡if ¡ they ¡are ¡OLD, ¡and ¡in ¡how ¡familiar ¡they ¡seem, ¡if ¡they ¡ are ¡NEW). ¡ ¡

• Let ¡us ¡code ¡S’s ¡response ¡as ¡‘1’, ¡if ¡‘old’, ¡and ¡‘0’, ¡if ¡

‘new’. ¡ ¡Binary ¡logisKc ¡regression ¡allows ¡us ¡ conveniently ¡to ¡explore ¡effects ¡of ¡various ¡predictors ¡

• f ¡S’s ¡response, ¡both ¡fixed-­‑effects ¡and ¡random-­‑

effects ¡predictors ¡that ¡vary ¡across ¡trials ¡(this ¡moreso ¡ than ¡TSD, ¡in ¡which ¡we ¡pool ¡data ¡across ¡trials ¡to ¡get ¡ pH ¡and ¡pFA). ¡

¡

24 ¡

Crossed random effects

• N subjects memorise m faces each under k
• conditions. Then ‘subno’ and ‘face’ are random

effects that are crossed factorially. For each subno*face, we have k repeated measures. E.g., (suid=1, faceid=f1) appears k times in Row

• 1. Ditto for (suid=i, faceid=fj)
• Suid f1

f1.c1 f2.c1 … fm.c1 f1 f1.c2 …fm.c2 … f1 f1.ck … fm.ck 1 2.1 3.5 … 2 4.0 1.3 … 3 … … N …

25

The ¡effects ¡of ¡an=cipatory ¡stress ¡on ¡ associa=ve ¡memory ¡retrieval ¡

Stephanie ¡Gagnon, ¡Anthony ¡Wagner ¡ November ¡2013 ¡

Steph ¡Gagnon’s ¡FYP, ¡2013 ¡

• Study ¡phase: ¡20 ¡Ss ¡study ¡64 ¡faces, ¡each ¡paired ¡with ¡an ¡

image ¡of ¡a ¡nature ¡scene ¡(i.e., ¡SCENE); ¡and ¡another ¡64 ¡ faces, ¡each ¡paired ¡with ¡an ¡image ¡of ¡an ¡object ¡(i.e., ¡ OBJECT). ¡Thus, ¡there ¡are ¡128 ¡OLD ¡faces. ¡

• Test ¡phase: ¡Ss ¡see ¡all ¡128 ¡OLD ¡faces ¡and ¡64 ¡NEW ¡faces. ¡

They ¡indicate ¡if ¡each ¡face ¡was ¡paired ¡in ¡the ¡study ¡phase ¡ with ¡a ¡scene ¡(‘scene’), ¡with ¡an ¡object ¡(‘object’), ¡with ¡ “either ¡a ¡scene ¡or ¡an ¡object ¡but ¡I ¡can’t ¡decide ¡ which” ¡(‘old’), ¡or ¡with ¡neither ¡a ¡scene ¡nor ¡an ¡object ¡ (‘new’). ¡ ¡Tes=ng ¡is ¡done ¡under ¡threat ¡of ¡shock ¡(‘threat’) ¡or ¡ not ¡(‘safe’); ¡this ¡is ¡a ¡within-­‑S ¡factor. ¡

• Simulated ¡data ¡are ¡in ¡the ¡file, ¡‘faces2.csv’. ¡To ¡simplify, ¡collapse ¡

‘scene’, ¡‘object’ ¡and ¡‘old’ ¡into ¡‘old’. ¡ ¡Thus ¡the ¡response ¡is ¡binary, ¡ ‘old’ ¡vs ¡‘new’. ¡

27 ¡

Methods: ¡Behavior ¡& ¡Physiology ¡

Study Phase Test Phase

Results: ¡Behavioral ¡

• Retrieval ¡of ¡object ¡associates ¡was ¡impaired ¡under ¡stress ¡(p<0.01), ¡while ¡

memory ¡for ¡place ¡associates ¡was ¡unaffected ¡(p>0.1; ¡interac=on ¡p=0.01). ¡

GLMM ¡

• Let ¡pn ¡= ¡P(‘old’) ¡on ¡the ¡n’th ¡trial ¡with ¡s=mulus, ¡S>mn ¡

(S>mn ¡= ¡0, ¡on ¡NEW ¡trials, ¡and ¡S>mn ¡= ¡1 ¡on ¡OLD ¡trials); ¡ and ¡let ¡logit(pn) ¡= ¡log(pn/(1 ¡– ¡pn)). ¡

• Then ¡the ¡binary ¡logis=c ¡model ¡is ¡

¡logit(pn) ¡= ¡b0 ¡+ ¡b1*S>mn ¡ ¡(*) ¡

• On ¡NEW ¡trials, ¡P(‘old’) ¡= ¡pFA, ¡by ¡defini=on, ¡and ¡on ¡OLD ¡

trials, ¡P(‘old’) ¡= ¡pH. ¡ ¡So, ¡on ¡puxng ¡S>mn ¡= ¡0 ¡ ¡and ¡S>mn ¡= ¡ 1, ¡respec=vely, ¡in ¡Eq. ¡(*), ¡we ¡get: ¡

¡logit(pFA) ¡= ¡b0 ¡, ¡and ¡

¡logit(pH) ¡= ¡b0 ¡+ ¡b1, ¡implying ¡that ¡ ¡logit(pH) ¡– ¡logit(pFA) ¡= ¡b1 ¡. ¡ This ¡is ¡very ¡convenient! ¡ ¡The ¡coefficient, ¡b1, ¡of ¡S>mn ¡in ¡ the ¡logis=c ¡regression ¡is ¡the ¡sensiKvity ¡parameter ¡ ¡

30 ¡

GLMM ¡

• To ¡generalise ¡by ¡including ¡Subject ¡& ¡Face, ¡let ¡pijn ¡= ¡

P(‘old’) ¡for ¡Subject ¡i ¡and ¡Face ¡j, ¡on ¡the ¡n’th ¡trial ¡with ¡ s=mulus, ¡S>mn; ¡and ¡let ¡logit(pijn) ¡= ¡log(pijn/(1 ¡– ¡pijn)). ¡

• Then ¡a ¡binary ¡logis=c ¡model ¡is, ¡e.g., ¡

¡logit(pijn) ¡= ¡b00 ¡+ ¡b0i ¡+ ¡aj ¡+ ¡(b10 ¡+ ¡b1i)*S>mn ¡(**) ¡

• In ¡this ¡model, ¡aj ¡is ¡a ¡random ¡variable ¡with ¡sd, ¡σf, ¡that ¡

indexes ¡varia=on ¡across ¡faces ¡in ¡ ‘memorability’ ¡(assumed ¡the ¡same ¡for ¡the ¡3 ¡types ¡of ¡ faces). ¡ ¡Similarly, ¡b0i ¡and ¡b1i ¡are ¡(possibly ¡correlated) ¡ random ¡variables ¡indexing ¡varia=on ¡across ¡Subjects. ¡

31 ¡

Analysis ¡of ¡‘faces2.csv’ ¡

32 ¡

suid ¡faceid ¡typesh ¡respshort ¡ ¡cond ¡ ¡ ¡ 1 ¡1 ¡SCENE ¡scene ¡ ¡safe ¡ 1 ¡2 ¡SCENE ¡scene ¡ ¡safe ¡ 1 ¡3 ¡SCENE ¡scene ¡ ¡… ¡ 1 ¡4 ¡SCENE ¡old ¡ ¡ 1 ¡5 ¡SCENE ¡object ¡ ¡ 1 ¡6 ¡SCENE ¡old ¡ ¡ 1 ¡7 ¡SCENE ¡scene ¡ ¡ 1 ¡8 ¡SCENE ¡scene ¡ ¡ 1 ¡9 ¡SCENE ¡new ¡ … ¡… ¡ ¡

Recode ‘typesh’ as OLD vs NEW, and ‘respshort’ as ‘old’ vs ‘new’.

33 ¡

d0 = read.csv('faces2.csv') d0$respsh2 = ifelse(d0$respshort == 'new', 0, 1) d0$typesh2 = ifelse(d0$typesh==‘NEW’,0,1) d0$suid = as.factor(d0$suid) d0$faceid = as.factor(d0$faceid) rs0 = with(d0, table(typesh, respshort))[stimtype1, responses1] print(rs0) print(round(prop.table(rs0, 1), 3))

‘sgagnon2c.r’

¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡respshort [cond = safe] typesh scene object old new SCENE 0.612 0.123 0.134 0.130 OBJECT 0.130 0.655 0.134 0.080 ¡ ¡ ¡ ¡ ¡ ¡NEW 0.133 0.127 0.098 0.642

¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡respshort [cond = threat]

typesh scene object old new SCENE 0.635 0.113 0.144 0.108 OBJECT 0.114 0.581 0.127 0.178 NEW 0.116 0.148 0.093 0.643

res1 = glmer(respsh2 ~ typesh2 + cond + (1 + typesh2 + cond | suid) + (1 | faceid), family = binomial, data = d0) Note: lmer(…, family = binomial, …) is

• deprecated. Use glmer() instead.

res4 = glmer(respsh2 ~ typesh2c * cond + (1 + typesh2c + cond | suid) + (1 | faceid), family = binomial, data = d0) res10 = glmer(respsh2 ~ typesh2c * cond + (1 | suid) + (0 + typesh2c | suid) + (1 | faceid), family = binomial, data = d0) print(anova(res10, res4))

Formula: respsh2 ~ typesh2 + cond + (1 + typesh2 + cond | suid) + (1 | faceid) AIC BIC logLik deviance 7145.058 7214.522 -3562.529 7125.058 Random effects: Groups Name Variance Std.Dev. Corr faceid (Intercept) 0.119250 0.34533 suid (Intercept) 0.023735 0.15406 typesh2 0.129813 0.36030 -1.00 condthreat 0.001014 0.03185 1.00 -1.00 Number of obs: 7680, groups: faceid, 192; suid, 20 Fixed effects: Estimate Std. Error z value Pr(>|z|) (Intercept) -0.51683 0.07502 -6.889 5.6e-12 *** typesh2 2.62445 0.11352 23.119 < 2e-16 *** condthreat -0.17845 0.06014 -2.967 0.003 **

Formula: respsh2 ~ typesh2 * cond + (1 + typesh2 + cond | suid) + (1 | faceid) AIC BIC logLik deviance 7138.205 7214.615 -3558.103 7116.205 Random effects: Groups Name Variance Std.Dev. Corr faceid (Intercept) 0.119692 0.34596 suid (Intercept) 0.021868 0.14788 typesh2 0.131206 0.36222 -1.00 condthreat 0.001787 0.04227 1.00 -1.00 Number of obs: 7680, groups: faceid, 192; suid, 20 Fixed effects: Estimate Std. Error z value Pr(>|z|) (Intercept) -0.602227 0.080076 -7.521 5.45e-14 *** typesh2 2.811920 0.130723 21.511 < 2e-16 *** condthreat -0.005388 0.083591 -0.064 0.94860 typesh2:condthreat -0.357988 0.119914 -2.985 0.00283 **

Formula: respsh2 ~ typesh2c * cond + (1 + typesh2c | suid) + (1 | faceid) AIC BIC logLik deviance 7132.697 7188.268 -3558.349 7116.697 Random effects: Groups Name Variance Std.Dev. Corr faceid (Intercept) 0.1197 0.3460 suid (Intercept) 0.0000 0.0000 typesh2c 0.1283 0.3582 NaN Number of obs: 7680, groups: faceid, 192; suid, 20 Fixed effects: Estimate Std. Error z value Pr(>|z|) (Intercept) 0.80040 0.05127 15.611 < 2e-16 *** typesh2c 2.80802 0.13016 21.573 < 2e-16 *** condthreat -0.18037 0.05996 -3.008 0.00263 ** typesh2c:condthreat -0.35373 0.11992 -2.950 0.00318 **

Removing ¡cond ¡as ¡a ¡random ¡slope ¡does ¡not ¡‘solve’ ¡the ¡corr ¡= ¡NaN ¡ ¡ problem, ¡probably ¡because ¡var(int) ¡= ¡0! ¡

Formula: respsh2 ~ typesh2c * cond + (1 | suid) + (0 + typesh2c | suid) + (1 | faceid) AIC BIC logLik deviance 7130.697 7179.322 -3558.349 7116.697 Random effects: Groups Name Variance Std.Dev. faceid (Intercept) 1.197e-01 3.460e-01 suid typesh2c 1.283e-01 3.582e-01 suid.1 (Intercept) 9.839e-10 3.137e-05 Number of obs: 7680, groups: faceid, 192; suid, 20 Fixed effects: Estimate Std. Error z value Pr(>|z|) (Intercept) 0.80041 0.05127 15.611 < 2e-16 *** typesh2c 2.80804 0.13016 21.574 < 2e-16 *** condthreat -0.18037 0.05996 -3.008 0.00263 ** typesh2c:condthreat -0.35372 0.11992 -2.950 0.00318 ** typesh2c is centered typesh2: typesh2c = (-½, ½) for (NEW, OLD). ‘Main’ effect of cond is affected.

Correlated ¡vs ¡Uncorrelated ¡random ¡effects ¡

• rs1 = glmer(resp ~ type*cond + (1 + type | suid) + (1

| faceid), data = d)

• rs2 = glmer(resp ~ type*cond + (1 | suid) + (0 + type

| suid) + (1 | faceid), data = d)

• lmer ¡and ¡glmer ¡assume ¡that ¡separate ¡random ¡effects, ¡i.e., ¡the ¡

terms ¡contained ¡in ¡round ¡brackets, ¡such ¡as, ¡(1 | suid) ¡, ¡are ¡

• independent. ¡ ¡In ¡rs2 ¡above, ¡the ¡‘int’ ¡in ¡(1 | suid) ¡and ¡the ¡

‘slope’ ¡in ¡(0 + type | suid) ¡are ¡assumed ¡to ¡be ¡

• uncorrelated. ¡ ¡Thus, ¡in ¡calcula=ng ¡the ¡likelihood ¡of ¡the ¡data, ¡the ¡

assump=on ¡of ¡independence ¡is ¡made, ¡and ¡the ¡likelihood ¡func=on ¡ contains ¡no ¡parameter ¡for ¡a ¡correla=on ¡between ¡the ¡‘int’ ¡and ¡ ‘slope’. ¡ ¡In ¡contrast, ¡in ¡rs1, ¡the ¡'int' ¡and ¡the ¡'slope' ¡are ¡assumed ¡to ¡ be ¡possibly ¡correlated ¡across ¡Subjects, ¡and ¡the ¡likelihood ¡func=on ¡ does ¡contain ¡a ¡parameter ¡for ¡the ¡correla=on ¡between ¡the ¡2 ¡

• effects. ¡
• Indeed, ¡anova(rs2, rs1) allows ¡us ¡to ¡test ¡whether ¡we ¡

should ¡include ¡a ¡correla=on ¡between ¡‘int’ ¡and ¡‘slope’ ¡in ¡our ¡ model ¡of ¡the ¡random ¡effects. ¡ ¡ ¡

• A ¡‘model’ ¡that ¡explains ¡why ¡corr(int, ¡slope) ¡= ¡-­‑1, ¡

whether ¡dummy ¡or ¡effect ¡coding ¡is ¡used. ¡

• Assume ¡all ¡Ss ¡have ¡the ¡same ¡P(‘old’ ¡| ¡OLD). ¡

logit(p) ¡ s=mulus ¡ Dummy ¡code: ¡ ¡ ¡ ¡ ¡0 ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡1 ¡ Effect ¡code: ¡ ¡ ¡ ¡ ¡ ¡-­‑0.5 ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡0 ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡0.5 ¡ ¡ ¡ ¡ ¡ ¡ ¡ Note ¡that ¡centering, ¡as ¡in ¡ Effect ¡coding, ¡does ¡not ¡solve ¡ the ¡corr(int, ¡slope) ¡= ¡-­‑1 ¡ ¡

• problem. ¡

What ¡value ¡is ¡added ¡by ¡considering ¡random ¡ effects? ¡

• Typically ¡our ¡main ¡interest ¡is ¡in ¡the ¡fixed ¡effects, ¡e.g., ¡

whether ¡the ¡effect ¡of ¡threat ¡is ¡the ¡same ¡for ¡object-­‑ tagged ¡and ¡scene-­‑tagged ¡faces? ¡ ¡It ¡follows ¡that ¡the ¡ main ¡value ¡of ¡LMM’s ¡is ¡in ¡accurately ¡assessing ¡the ¡ fixed ¡effects, ¡given ¡the ¡correla=ons ¡in ¡the ¡responses ¡ that ¡result ¡from ¡repeatedly ¡measuring ¡Subjects ¡and ¡

• Faces. ¡ ¡We ¡can ¡be ¡confident ¡that ¡our ¡conclusions ¡

about ¡the ¡fixed ¡effects ¡generalise ¡to ¡the ¡popula=on ¡of ¡ Subjects ¡and ¡to ¡the ¡popula=on ¡of ¡Faces. ¡

• There ¡are ¡secondary ¡benefits, ¡as ¡well, ¡from ¡

construc=ng ¡the ¡‘best-­‑fixng’ ¡random ¡effects ¡model. ¡ ¡ E.g., ¡in ¡what ¡ways ¡do ¡Subjects ¡vary? ¡

In ¡what ¡ways ¡do ¡Subjects ¡vary? ¡ ¡

• In ¡this ¡study, ¡Subjects ¡can ¡vary ¡in ¡their ¡

‘intercept’ ¡(e.g., ¡their ¡false-­‑alarm ¡rates), ¡their ¡ sensi=vity ¡to ¡the ¡OLD ¡vs ¡NEW ¡categories ¡of ¡faces ¡(the ¡ ‘slope’, ¡d’), ¡or ¡their ¡sensi=vity ¡to ¡the ¡safe ¡vs ¡threat ¡ condi=ons ¡(another ¡‘slope’ ¡parameter). ¡ ¡Do ¡Subjects ¡ vary ¡appreciably ¡in ¡all ¡3 ¡parameters? ¡ ¡Is ¡the ¡varia=on ¡ across ¡Subjects ¡in ¡‘intercept’ ¡correlated ¡with ¡that ¡in ¡ ‘slope’? ¡

• These ¡ques=ons ¡can ¡mo=vate ¡future ¡studies ¡of ¡(i) ¡

the ¡(e.g., ¡neural, ¡physiological) ¡antecedents ¡or ¡ correlates ¡of ¡individual ¡differences ¡in ¡performance; ¡ and ¡(ii) ¡the ¡physical ¡or ¡social ¡aspects ¡of ¡a ¡face ¡that ¡ correlate ¡with ¡the ¡‘memorability’ ¡of ¡a ¡face. ¡

P-­‑values ¡

• For ¡con=nuity ¡with ¡an ¡earlier ¡lecture, ¡here ¡are ¡a ¡few ¡

comments ¡on ¡p-­‑values. ¡This is a relatively unimportant topic, but one that may arise in the journal review process.

• The ¡following ¡exchange ¡(Feb-­‑Aug, ¡2012) ¡is ¡

informa=ve. ¡ ¡Source: ¡ http://stats.stackexchange.com/questions/22988/ significant-effect-in-lme4-mixed-model

• ECII: ¡I ¡use ¡lme4 ¡in ¡R ¡to ¡fit ¡the ¡mixed ¡model ¡

¡lmer(value~status+(1|experiment))) ¡

• where ¡value ¡is ¡con=nuous, ¡status ¡and ¡experiment ¡

are ¡factors. ¡How ¡can ¡I ¡know ¡that ¡the ¡effect ¡of ¡ status ¡is ¡significant? ¡R ¡reports ¡only ¡t-­‑values ¡and ¡ not ¡p-­‑values. ¡ ¡

• Ben ¡Bolker: ¡I ¡wd ¡suggest ¡lme(), ¡because ¡you ¡

don't ¡need ¡any ¡of ¡the ¡stuff ¡that ¡lmer() ¡offers ¡ (higher ¡speed, ¡handling ¡of ¡crossed ¡random ¡effects, ¡ GLMMs ¡...). ¡lme() ¡should ¡give ¡you ¡exactly ¡the ¡ same ¡coefficient ¡and ¡variance ¡es=mates ¡but ¡will ¡ also ¡compute ¡df ¡and ¡p-­‑values ¡for ¡you ¡(which ¡do ¡ make ¡sense ¡in ¡a ¡"classical" ¡design ¡such ¡as ¡you ¡ appear ¡to ¡have). ¡ ¡

• You ¡may ¡also ¡want ¡to ¡consider ¡the ¡random ¡term ¡

~status|experiment ¡(allowing ¡for ¡variaKon ¡of ¡ status ¡effects ¡across ¡blocks, ¡or ¡equivalently ¡ including ¡a ¡status-­‑by-­‑experiment ¡interacKon). ¡ Posters ¡above ¡are ¡also ¡correct ¡that ¡your ¡t ¡sta=s=cs ¡ are ¡so ¡large ¡that ¡your ¡p-­‑value ¡will ¡definitely ¡be ¡ <0.05, ¡but ¡I ¡can ¡imagine ¡you ¡would ¡like ¡"real" ¡p-­‑

• values. ¡

library(nlme) m1 <- lme(value ~ status, random = ~ 1 | experiment, data=mydata) anova(m1) m2 <- lme(value ~ status, random = ~ status | experiment, data=mydata)

• John: ¡I ¡don't ¡know ¡about ¡this ¡answer. ¡lmer ¡could ¡

just ¡as ¡easily ¡report ¡the ¡same ¡kinds ¡of ¡p-­‑values ¡but ¡ doesn't ¡for ¡valid ¡reasons. ¡I ¡guess ¡it's ¡the ¡comment ¡ that ¡there ¡are ¡any ¡"real" ¡p-­‑values ¡here ¡that ¡bugs ¡

• me. ¡You ¡could ¡argue ¡that ¡you ¡can ¡find ¡one ¡possible ¡

cutoff, ¡and ¡that ¡any ¡reasonable ¡cutoff ¡is ¡passed. ¡ But ¡you ¡can't ¡argue ¡there's ¡a ¡real ¡p-­‑value. ¡

• [In ¡another ¡post] ¡… ¡you ¡therefore ¡shouldn't ¡report ¡

an ¡exact ¡p-­‑value, ¡only ¡that ¡your ¡cutoff ¡is ¡passed. ¡… ¡ You ¡should ¡really ¡consider ¡the ¡[model ¡comparison ¡ approach] ¡because ¡the ¡whole ¡idea ¡of ¡just ¡s=cking ¡ with ¡a ¡pass ¡point ¡for ¡p-­‑values ¡as ¡the ¡final ¡and ¡most ¡ important ¡informa=on ¡to ¡extract ¡from ¡your ¡data ¡is ¡ generally ¡misguided. ¡

• For testing fixed effects with t, we need the df of

the error MS in the denominator of t. Although

• ne can put bounds on this df, in general, the

exact value of the df is not known. Thus exact p- values cannot be calculated.

• This is what is meant by no ‘real’ p-values.

Bolker’s point is that, for balanced designs in which one can assume Compound Symmetry,

• ne can calculate the right error MS and its df.

(We in Psy 252 have these EMS Tables in slides from previous years, but they are now deprecated!) In these cases, we do get ‘real’ p- values.

• Ben ¡Bolker: ¡ ¡For ¡a ¡classical ¡design ¡(balanced, ¡

nested, ¡etc.) ¡I ¡think ¡I ¡can ¡indeed ¡argue ¡that ¡ there's ¡a ¡real ¡p-­‑value, ¡i.e. ¡a ¡probability ¡of ¡ gexng ¡an ¡es=mate ¡of ¡beta ¡of ¡an ¡observed ¡ magnitude ¡or ¡greater ¡if ¡the ¡null ¡hypothesis ¡ (beta=0) ¡were ¡true ¡... ¡lme4 ¡doesn't ¡provide ¡ these ¡denominator ¡df, ¡I ¡believe, ¡because ¡it's ¡ harder ¡to ¡detect, ¡in ¡general ¡from ¡an ¡lme4 ¡ model ¡structure, ¡when ¡the ¡model ¡specified ¡is ¡

• ne ¡where ¡some ¡heuris=c ¡for ¡compu=ng ¡a ¡

classical ¡denominator ¡df ¡would ¡work ¡... ¡

50

• The current thinking in mixed-effects modeling is

that “H0: σb

2 = 0” is not meaningful (How can a

variance ever be 0?!). Thus, F tests of this null are not meaningful.

• Model comparison: Is the likelihood of the data

is sig greater if the parameter, σb

2, were included

in the model than if it were not. That is, one shd test these 2 nested models to see if the more complex model that includes the extra parameter, σb

2, is warranted. For this, use the

• utput for the 2 models, noting that higher

logLik is better:,

• Or lower -logLik is better!

AIC BIC

logLik 104.3624 107.7689 -49.18119

51

AIC BIC logLik 104.3624 107.7689 -49.18119

• Differences in -2*logLik are distributed as chi-

square with df = difference in # of parameters.

• Often the interesting models can be arrayed from

least to most complex (i.e., can be nested). If so, using -2*logLik to compare models is satisfying because one has a formal test.

• But interesting models are sometimes not nested.

In that case, one can use AIC or BIC to compare models (lower is better, noting -3 < -2!). AIC and BIC penalise models for (i) poor fit (indexed by ‘deviance’ = logLik), and (ii) complexity (indexed by the # of parameters).

52

Comparing Mixed Models: Benoît’s guide

You ¡always ¡want ¡deviance ¡/ ¡AIC ¡to ¡be ¡low. ¡Differences ¡in ¡deviance ¡are ¡tested ¡as ¡a ¡ chi-­‑square ¡with ¡the ¡difference ¡in ¡# ¡of ¡parameters ¡as ¡the ¡df’s. ¡

Are the 2 models nested? Use AIC

No

Same fixed effect structure? Use REML/ML Deviance Use ML Deviance

No Yes Yes

No formal test Chi- square

Week 9, Lec 2: Issues in Mixed Models

• More examples of mixed models to prepare us

for the issues, some novel, that we’ll confront in

• ur own research. Feel free to consult the lme4

Help docs, listserv’s. etc.

• A Generalised LMM (GLMM) analysis and some

complications, e.g., corr(int, slope) = ±1.

• More on p-values.
• How to implement and test/compare models in

R? Use of ML vs REML deviance.

• An example with nested factors

53

On ¡Jan ¡2, ¡2013, ¡at ¡4:47 ¡AM, ¡Carlos ¡Gias ¡wrote: ¡ Hi, ¡ ¡I ¡am ¡new ¡to ¡mixed ¡models ¡and ¡not ¡sure ¡how ¡to ¡design ¡a ¡model ¡for ¡my ¡data ¡to ¡ use ¡with ¡the ¡lmer ¡func=on. ¡ ¡ I ¡am ¡trying ¡to ¡study ¡the ¡effect ¡of ¡a ¡drug ¡(treatment) ¡in ¡subjects ¡along ¡=me. ¡For ¡ every ¡subject ¡we ¡measure ¡ac=vity ¡from ¡a ¡(variable) ¡number ¡of ¡neurons ¡at ¡ different ¡=me ¡points. ¡Therefore, ¡the ¡neuron ¡measurement ¡is ¡nested ¡within ¡

• subject. ¡I ¡would ¡also ¡like ¡to ¡use ¡the ¡baseline ¡ac=vity ¡measurement ¡as ¡a ¡covariate ¡

for ¡each ¡neuron. ¡This ¡is ¡a ¡small ¡sample ¡of ¡the ¡dataset: ¡[ET’s ¡Note: ¡The ¡natural ¡ (only?) ¡code ¡for ¡neurons ¡wd ¡be ¡1, ¡2, ¡… ¡for ¡each ¡neuron.] ¡

subject treatment neuron baseline time activity

1 1 1 3.06 1 7.02 1 1 1 3.06 2 6 1 1 1 3.06 3 9 1 1 2 3.00 1 5 1 1 2 3.00 2 6 1 1 2 3.00 3 9 2 2 3 4.77 1 3 2 2 3 4.77 2 2 2 2 3 4.77 3 1 3 2 3 2.14 1 2.03 3 2 3 2.14 2 2 3 2 3 2.14 3 1.5

I ¡was ¡wondering ¡if ¡the ¡following ¡would ¡be ¡an ¡appropriate ¡model. ¡ ¡ ac=vity ¡~ ¡baseline ¡+ ¡treatment*=me ¡+ ¡(1|subject:neuron) ¡+ ¡(1|=me) ¡ ¡ ¡ I ¡am ¡also ¡wondering ¡how ¡I ¡could ¡deal ¡with ¡the ¡problem ¡of ¡missing ¡data. ¡Is ¡ there ¡a ¡func=on/package ¡that ¡could ¡be ¡helpful ¡in ¡this ¡design? ¡ ¡ I ¡hope ¡you ¡can ¡help. ¡ ¡ ¡ Best ¡regards, ¡ ¡ Carlos ¡Gias ¡

¡

carlos, ¡ ¡

• 1. ¡ ¡are ¡the ¡=me ¡intervals ¡between ¡=me=1, ¡=me=2 ¡and ¡=me=3 ¡

known ¡to ¡you, ¡and ¡are ¡they ¡the ¡same ¡across ¡subjects? ¡ ¡that ¡is, ¡do ¡ you ¡expect ¡a ¡=me ¡effect ¡of ¡treatment? ¡ ¡it ¡appears ¡that ¡you ¡do ¡ expect ¡an ¡effect ¡-­‑ ¡is ¡it ¡linear ¡or ¡quadra=c? ¡ ¡e.g., ¡maybe ¡'baseline' ¡is ¡ '=me ¡= ¡0', ¡'=me=1' ¡is ¡'=me ¡= ¡10 ¡minutes', ¡'=me=2' ¡is ¡'=me ¡= ¡20 ¡ minutes', ¡etc. ¡ ¡if ¡this ¡is ¡true, ¡i ¡wd ¡not ¡model ¡'=me' ¡as ¡a ¡random ¡ effect, ¡only ¡as ¡a ¡fixed ¡effect, ¡and ¡i ¡wd ¡include ¡the ¡baseline ¡ac=vity ¡ as ¡part ¡of ¡the ¡=me-­‑series ¡for ¡each ¡neuron. ¡ ¡it's ¡unlikely ¡that ¡'=me' ¡ is ¡best ¡modeled ¡as ¡a ¡random ¡effect, ¡but ¡maybe ¡i'm ¡missing ¡

• something. ¡

¡

• 2. ¡ ¡for ¡the ¡random ¡effects, ¡i ¡might ¡model ¡'neuron' ¡as ¡nested ¡

within ¡'subject', ¡altho ¡'subject:neuron' ¡might ¡be ¡jus=fiable ¡and ¡

• interpretable. ¡

i ¡wd ¡try ¡the ¡following ¡nested ¡models: ¡ ¡ Model1: ¡ ¡ac=vity ¡~ ¡treatment ¡* ¡=me ¡+ ¡(1 ¡| ¡subject ¡/ ¡neuron) ¡ Model1a: ¡ac=vity ¡~ ¡treatment ¡* ¡=me ¡+ ¡(1 ¡+ ¡=me ¡| ¡subject ¡/ ¡neuron) ¡ ¡ Model2: ¡ac=vity ¡~ ¡treatment ¡* ¡poly(=me, ¡2) ¡+ ¡(1 ¡| ¡subject ¡/ ¡neuron) ¡ Model2a: ¡ac=vity ¡~ ¡treatment ¡* ¡poly(=me, ¡2) ¡+ ¡(1 ¡+ ¡=me ¡| ¡subject ¡/ ¡ neuron) ¡ Model2b: ¡ac=vity ¡~ ¡treatment ¡* ¡poly(=me, ¡2) ¡+ ¡(1 ¡+ ¡poly(=me, ¡2) ¡| ¡ subject ¡/ ¡neuron) ¡ ¡ good ¡luck! ¡ ¡ ¡ ewart ¡

¡

ewart, ¡thanks ¡for ¡your ¡quick ¡reply. ¡It ¡has ¡been ¡extremely ¡useful. ¡ ¡ To ¡answer ¡your ¡ques=ons ¡

• 1. ¡ ¡are ¡the ¡Kme ¡intervals ¡between ¡Kme=1, ¡Kme=2 ¡and ¡Kme=3 ¡

known ¡to ¡you, ¡and ¡are ¡they ¡the ¡same ¡across ¡subjects? ¡that ¡is, ¡do ¡ you ¡expect ¡a ¡Kme ¡effect ¡of ¡treatment? ¡ ¡

• ­‑ ¡The ¡intervals ¡are ¡the ¡same ¡across ¡subjects ¡and ¡they ¡are ¡known ¡(in ¡

fact ¡your ¡guess ¡was ¡right!): ¡1-­‑> ¡10mins, ¡2-­‑>20 ¡mins... ¡etc ¡(baseline ¡-­‑> ¡ 0mins). ¡ ¡

• ­‑ ¡Yes, ¡I ¡am ¡expec=ng ¡to ¡see ¡a ¡=me ¡effect ¡of ¡treatment. ¡

¡ is ¡it ¡linear ¡or ¡quadraKc? ¡ ¡

• ­‑ ¡This ¡is ¡not ¡known ¡in ¡advance ¡so ¡I ¡would ¡start ¡assuming ¡the ¡simplest ¡
• case. ¡ ¡

¡

¡

¡if ¡this ¡is ¡true, ¡i ¡wd ¡not ¡model ¡'Kme' ¡as ¡a ¡random ¡effect, ¡only ¡as ¡a ¡ fixed ¡effect, ¡ ¡

• ­‑ ¡Thanks ¡for ¡the ¡sugges=on. ¡I ¡think ¡you ¡are ¡right, ¡it ¡would ¡be ¡more ¡

appropriate ¡to ¡model ¡it ¡as ¡a ¡fixed ¡effect. ¡ ¡ ¡ ¡i ¡wd ¡include ¡the ¡baseline ¡acKvity ¡as ¡part ¡of ¡the ¡Kme-­‑series ¡for ¡each ¡

• neuron. ¡ ¡
• ­‑ ¡Could ¡this ¡improve ¡the ¡analysis ¡in ¡some ¡way? ¡In ¡any ¡case, ¡due ¡to ¡

adequate ¡randomiza=on, ¡I ¡am ¡not ¡expec=ng ¡a ¡strong ¡effect ¡of ¡using ¡ the ¡baseline ¡measurement ¡in ¡the ¡model. ¡ ¡ ¡

• 2. ¡ ¡for ¡the ¡random ¡effects, ¡i ¡might ¡model ¡'neuron' ¡as ¡nested ¡within ¡

'subject', ¡altho ¡'subject:neuron' ¡might ¡be ¡jusKfiable ¡and ¡

• interpretable. ¡
• ­‑ ¡I ¡actually ¡meant ¡to ¡use ¡the ¡nested ¡symbol ¡(/) ¡so ¡thanks ¡for ¡poin=ng ¡

this ¡out! ¡ ¡

I ¡hope ¡I ¡am ¡not ¡abusing ¡your ¡=me ¡but ¡I ¡would ¡be ¡grateful ¡if ¡you ¡had ¡any ¡sugges=ons ¡

• n ¡missing ¡values. ¡My ¡ini=al ¡thought ¡was ¡to ¡use ¡the ¡program ¡Amelia ¡II ¡(from ¡Gary ¡

King) ¡but ¡it ¡doesn't ¡handle ¡a ¡nested ¡structure. ¡Are ¡you ¡aware ¡of ¡any ¡packages ¡that ¡ could ¡help? ¡ Regards, ¡Carlos ¡ ¡ Date: ¡ ¡January ¡2, ¡2013 ¡9:34:00 ¡AM ¡PST ¡ carlos, ¡one ¡of ¡the ¡virtues ¡of ¡lmer() ¡is ¡that ¡it ¡accepts ¡different ¡numbers ¡of ¡

• bserva=ons ¡from ¡different ¡subjects. ¡ ¡when ¡some ¡subjects ¡have ¡less ¡data ¡than ¡
• thers, ¡we ¡say ¡the ¡design ¡is ¡'unbalanced'; ¡it ¡isn't ¡that ¡data ¡are ¡'missing' ¡from ¡an ¡
• therwise ¡'balanced' ¡design. ¡ ¡so ¡the ¡issue ¡doesn't ¡arise. ¡

¡ maybe ¡you ¡cd ¡try ¡using ¡the ¡baseline ¡data ¡in ¡a ¡couple ¡of ¡ways, ¡and ¡see ¡what ¡

• happens. ¡ ¡one ¡way ¡is ¡to ¡treat ¡it ¡as ¡the ¡=me=0 ¡observa=on ¡in ¡a ¡=me-­‑series. ¡ ¡another ¡

wd ¡be ¡to ¡use ¡it ¡as ¡a ¡covariate. ¡ ¡i ¡just ¡saw ¡onkelinx's ¡reply ¡sugges=ng ¡offset(baseline); ¡ he's ¡an ¡expert, ¡so ¡you ¡shd ¡try ¡this ¡approach ¡as ¡well. ¡ ¡good ¡luck! ¡ ewart ¡ [ET’s ¡note: ¡offset(X) ¡is ¡used ¡as ¡a ¡term ¡in ¡a ¡linear ¡predictor ¡to ¡s=pulate ¡that ¡X ¡is ¡a ¡ predictor ¡in ¡the ¡regression ¡with ¡a ¡known ¡coeff ¡of ¡1, ¡rather ¡than ¡an ¡es=mated ¡coeff.] ¡

61

We ¡end ¡with ¡a ¡slide ¡from ¡Week ¡1, ¡Lec ¡1, ¡sugges=ng ¡how ¡ ¡

• ur ¡interests ¡ought ¡to ¡affect ¡our ¡research: ¡

Workhorse ¡func=ons, ¡depending ¡on ¡your ¡design: ¡ ¡

rs1 = lm(score ~ time + group, d) rs2 = glm(score0 ~ time + group, family = ‘binomial’, d) rs3 = lmer(score ~ time + group + (1 | suid), d) rs4 = glmer(score0 ~ time + group + (1 | suid), family = binomial, d)

Week 9, Lec 3: Issues in Mixed Models

• When to, and why, use a mixed model? (Helpful

prototype: indep groups vs related groups t-test)

• Plausible (mixed) models for analysing the data?

The rationale for these models?

• Implement and test/compare models in R, using

ML vs REML deviance.

• Is a fixed effect sig? Use rule-of-thumb: ‘If |t| >

2, reject H0’. What about p-values?

• Interpret random effects, i.e., vars & covars.
• Theory of Maximum Likelihood (ML) estimation

62

63

Attention SubjID Numsol=1 Numsol=2 Numsol=3 1.00 2.00 4.00 7.00 2.00 3.00 4.00 5.00 3.00 3.00 5.00 6.00 4.00 5.00 7.00 5.00 Attn = 1 5.00 4.00 5.00 8.00 6.00 5.00 5.00 6.00 7.00 5.00 4.50 6.00 8.00 5.00 7.00 8.00 9.00 2.00 3.00 7.00 10.00 6.00 5.00 6.00 1.00 6.00 5.00 6.00 2.00 8.00 9.00 8.00 3.00 6.00 5.00 9.00 4.00 8.00 8.00 7.00 Attn = 2 5.00 8.00 8.00 7.00 6.00 6.00 8.00 7.00 7.00 7.00 7.00 6.00 8.00 7.00 8.00 6.00 9.00 5.00 6.00 6.00 10.00 6.00 6.00 5.00

Reprise ¡of ¡ ‘kv.csv’ ¡

Fixed vs Random effects

• The fixed effects are the attn and nsol main

effects, and the attn*nsol interaction.

• ‘suid’ is the random effects factor. What

variables are assumed to differ across Ss?

• For levels 1 and 2 of attn, let a1 and a2 (a1 + a2

= 0) be the fixed effects of attn; let b1 and b2 be the effects of nsol (assumed quantitative) within levels 1 and 2, respectively, of attn.

• Subject j at level 1 of attn has an ‘intercept’ of

a1 + a1j , and we assume that a1j is a random effect that is Normally distributed with mean 0 and variance σ0

• 2. Ditto for a2j.

64

65

‘Intercept only’ and ‘intercept + slope’ models for the effect of nsol on score: 3 randomly chosen Ss

Score ¡ nsol ¡ nsol ¡

Plot each S’s data to see which model is more appropriate.

• We might also assume that Ss vary in the

effect of nsol on score. For Subject j at level 1

• f attn, b1 + b1j is the effect of nsol, where b1j is

a random effect (‘slope’) that is N(0, σ1

2). Ditto

for b2j.

• If Subjects vary on both a1j and b1j, then it is

likely that, across Ss, there will be a covariance, σ01, between a1j and b1j. Thus, the random effects parameters of the full ‘intercept + slope’ model are σ0

2, σ1 2, σ01; as well as the

error variance, σe

2.

• The 2 models we consider are the simple,

‘intercept only’ model, and the less simple, ‘intercept + slope’ model.

66

67

These plots suggest that there is appreciable variation in both intercept and slope. In addition, one can see the fixed effects.