[R-sig-ME] gamm4 model formulation clarification

[Reinhold Kliegl]

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
>