[R-sig-ME] gamm4 model formulation clarification

[Reinhold Kliegl]

Sorry, the following line should be:
constrasts(ab$g) <- cmat

On Sat, Jul 31, 2010 at 11:30 PM, Reinhold Kliegl
<reinhold.kliegl at gmail.com> wrote:
> Perhaps assign contrasts to g. For example,
> # ... ... contrast estimates
> cmat <- matrix(c(  -1/2, -1/2, +1/2, +1/2,   # Main effect 1
>                              -1/2, +1/2, -1/2, +1/2,   # Main effect 2
>                              +1/2, -1/2, -1/2, +1/2),  4,  3,    # Interaction
>                               dimnames=list(c("A1", "A2", "A3", "A4"),
>                                c(".34-12", ".24-13", ".14-23")))
>
> constrasts(g) <- cmat
>  fit = gamm4(
>      data = ab
>      , formula = rrt ~ g+s(soa,by=g)
>     , random = ~ (1|id)
>     , family = 'gaussian'
>   )
>
> Reinhold Kliegl
>
> On Sat, Jul 31, 2010 at 6:10 AM, Mike Lawrence <Mike.Lawrence at dal.ca> wrote:
>> 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. ~
>>
>> _______________________________________________
>> R-sig-mixed-models at r-project.org mailing list
>> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
>>
>