[R-sig-ME] gamm4 model formulation clarification

[Mike Lawrence]

Hi folks,

I have some data that form a 2x2x15 design, where the 15 level
variable is a discrete sampling of a ratio variable with clear
non-linearities (see bottom for a dput() of the means). I came across
gamm4 tonight and it looks like it will help tackle this data, but I'm
not sure how to tell it to let the smooth of the 15 level variable
vary as a function of BOTH of the other predictor variables. As a
hack, I created a dummy 4 level variable that represented the
combination of the 2x2 level variables, but I'm not positive that this
was the right thing to do. Any feedback would be greatly appreciated.
Here's how I had things set up:

> str(ab)
'data.frame':	49668 obs. of  5 variables:
 $ id       : Factor w/ 20 levels "5","6","7","8",..: 1 1 1 1 1 1 1 1 1 1 ...
 $ design   : Factor w/ 2 levels "CD","NCD": 1 1 1 1 1 1 1 1 1 1 ...
 $ ddB      : Factor w/ 2 levels "0ddB","+ddB": 1 1 1 1 1 1 1 1 1 1 ...
 $ soa      : num  300 300 300 300 300 300 300 300 300 300 ...
 $ rt       : num  441 373 440 290 221 ...
 $ g        : Factor w/ 4 levels "+ddB CD","+ddB NCD",..: 3 3 3 3 3 3
3 3 3 3 ...

#id is the random effect (human participants)
#design and ddB are 2-level fixed effects
#soa is a 15 level fixed effect
#rt is the data to predict (there are dozens of observations for each
combination of the random and fixed effects)

#fit using dummy variable "g" to get different soa smooths per
combination of ddB and design
fit = gamm4(
	data = ab
	, formula = rrt ~ ddB*design+s(soa,by=g)
	, random = ~ (1|id)
	, family = 'gaussian'
)


#here's dput() ouput of the 2x2x15 means, revealing that soa is
clearly non-linear
> dput(means)
structure(list(ddB = structure(c(1L, 2L, 1L, 2L, 1L, 2L, 1L,
2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L,
2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L,
2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L,
2L, 1L, 2L, 1L, 2L), .Label = c("0ddB", "+ddB"), class = "factor"),
    design = structure(c(1L, 1L, 2L, 2L, 1L, 1L, 2L, 2L, 1L,
    1L, 2L, 2L, 1L, 1L, 2L, 2L, 1L, 1L, 2L, 2L, 1L, 1L, 2L, 2L,
    1L, 1L, 2L, 2L, 1L, 1L, 2L, 2L, 1L, 1L, 2L, 2L, 1L, 1L, 2L,
    2L, 1L, 1L, 2L, 2L, 1L, 1L, 2L, 2L, 1L, 1L, 2L, 2L, 1L, 1L,
    2L, 2L, 1L, 1L, 2L, 2L), .Label = c("CD", "NCD"), class = "factor"),
    soa = c(-175, -175, -175, -175, -125, -125, -125, -125, -75,
    -75, -75, -75, -25, -25, -25, -25, 25, 25, 25, 25, 75, 75,
    75, 75, 125, 125, 125, 125, 175, 175, 175, 175, 250, 250,
    250, 250, 350, 350, 350, 350, 450, 450, 450, 450, 550, 550,
    550, 550, 700, 700, 700, 700, 900, 900, 900, 900, 1100, 1100,
    1100, 1100), rt = structure(c(444.93273739708, 441.513208123373,
    471.546335977687, 472.283526755609, 444.056771928409, 442.461563636396,
    472.166623352385, 474.462330421383, 445.513188139913, 444.966221088285,
    475.669898472348, 461.315736508479, 442.982086750377, 435.355068015326,
    458.89264807265, 455.638816561315, 434.340063293089, 415.709247937572,
    450.359216604288, 438.114012378908, 420.982354512772, 404.811090521404,
    442.721222728774, 430.421981079446, 406.373845346035, 393.411137886923,
    439.687373941982, 427.611295261772, 398.185439780545, 384.007719475387,
    429.451304040456, 431.531830923294, 390.486982286177, 380.880066145206,
    431.287499227013, 435.926925139256, 382.937376664574, 382.020452210116,
    439.606939778423, 441.751714989440, 384.546631449636, 385.477927260379,
    440.7286749068, 442.790235121405, 387.589278670798, 386.68589658312,
    440.384534575679, 440.376618739428, 392.649929780933, 392.496542853778,
    442.673509587221, 446.148838198613, 401.042352162140, 394.475542684099,
    441.072702415870, 440.157897723193, 406.059905100442, 398.292572833949,
    443.899793077701, 445.715779415161), .Dim = c(60L, 1L), .Dimnames = list(
        NULL, NULL))), .Names = c("ddB", "design", "soa", "rt"
), row.names = c(NA, 60L), class = "data.frame")

#visualise the means
library(ggplot2)
ggplot(
	data = means
	, mapping = aes(
		x = soa
		, y = rt
		, colour = design
		, linetype = ddB
	)
)+
geom_line()+
geom_point(shape=21,fill='white')

-- 
Mike Lawrence
Graduate Student
Department of Psychology
Dalhousie University

Looking to arrange a meeting? Check my public calendar:
http://tr.im/mikes_public_calendar

~ Certainty is folly... I think. ~