www.phidot.org

I have been looking at the DM's for an analysis I'm working on and I'm confused as to how an interaction (:) model vs. a fully-crossed model (*) should look like (I don't know if this makes a difference but "ev" is a time-varying individual covariate).

Fully crossed model

Code: Select all

> ev.res[[4]]$design.matrix # S(~winter * ev)f(~1)
                      S:(Intercept) S:winter S:ev   S:winter:ev f:(Intercept)
S gfaA0 c1 c1 a0 t1   "1"           "1"      "ev1"  "ev1"       "0"         
S gfaA0 c1 c1 a1 t2   "1"           "0"      "ev2"  "0"         "0"         
S gfaA0 c1 c1 a2 t3   "1"           "0"      "ev3"  "0"         "0"         
S gfaA0 c1 c1 a3 t4   "1"           "0"      "ev4"  "0"         "0"         
S gfaA0 c1 c1 a4 t5   "1"           "1"      "ev5"  "ev5"       "0"         
S gfaA0 c1 c1 a5 t6   "1"           "0"      "ev6"  "0"         "0"         
S gfaA0 c1 c1 a6 t7   "1"           "0"      "ev7"  "0"         "0"         
S gfaA0 c1 c1 a7 t8   "1"           "0"      "ev8"  "0"         "0"         
S gfaA0 c1 c1 a8 t9   "1"           "1"      "ev9"  "ev9"       "0"         
S gfaA0 c1 c1 a9 t10  "1"           "0"      "ev10" "0"         "0"         
S gfaA0 c1 c1 a10 t11 "1"           "0"      "ev11" "0"         "0"         
S gfaA0 c1 c1 a11 t12 "1"           "0"      "ev12" "0"         "0"         
S gfaA0 c1 c1 a12 t13 "1"           "1"      "ev13" "ev13"      "0"         
S gfaA0 c1 c1 a13 t14 "1"           "0"      "ev14" "0"         "0"         
S gfaA0 c1 c1 a14 t15 "1"           "0"      "ev15" "0"         "0"         
S gfaA0 c1 c1 a15 t16 "1"           "0"      "ev16" "0"         "0"         
S gfaA0 c1 c1 a16 t17 "1"           "1"      "ev17" "ev17"      "0"         
S gfaA0 c1 c1 a17 t18 "1"           "0"      "ev18" "0"         "0"         
S gfaA0 c1 c1 a18 t19 "1"           "0"      "ev19" "0"         "0"         
S gfaA0 c1 c1 a19 t20 "1"           "0"      "ev20" "0"         "0"         
S gfaA0 c1 c1 a20 t21 "1"           "1"      "ev21" "ev21"      "0"         
S gfaA0 c1 c1 a21 t22 "1"           "0"      "ev22" "0"         "0"         
S gfaA0 c1 c1 a22 t23 "1"           "0"      "ev23" "0"         "0"         
f gfaA0 c1 c1 a0 t1   "0"           "0"      "0"    "0"         "1"     

Interaction model

Code: Select all

> ev.res[[6]]$design.matrix # S(~winter:ev)f(~1)
                      S:(Intercept) S:winter:ev f:(Intercept)
S gfaA0 c1 c1 a0 t1   "1"           "ev1"       "0"         
S gfaA0 c1 c1 a1 t2   "1"           "0"         "0"         
S gfaA0 c1 c1 a4 t5   "1"           "ev5"       "0"         
S gfaA0 c1 c1 a8 t9   "1"           "ev9"       "0"         
S gfaA0 c1 c1 a12 t13 "1"           "ev13"      "0"         
S gfaA0 c1 c1 a16 t17 "1"           "ev17"      "0"         
S gfaA0 c1 c1 a20 t21 "1"           "ev21"      "0"         
f gfaA0 c1 c1 a0 t1   "0"           "0"         "1"   

Why doesn't the "winter:ev" column in the interaction model look the same as in the fully-crossed model (i.e., why is there a zero in row 2, column 2 of the interaction model?)?

You don't say what winter is but from the behavior I'm presuming that winter is a 0/1 numeric variable and if that is the case the interaction ev:winter does what it is supposed to. Interactions with numeric variables is simply multiplication which is unlike interactions of 2 factor variables. In this case, it fits only the ev interaction term for winter=1. The ev main effect is excluded here. Thus this model restricts the affect of ev to winter only and for winter=0, it would use the intercept.

--jeff

Wow you are fast Jeff, thank you!

Let's see... the winter term is design data that I added to ddl$S. It is a seasonal model ("Brownie", four seasons/year), and winter happens once a year. I guess it is numeric, should it be a factor?

I don't want it to have the main effect for "ev", which is the number of "severe winter events" (<-1°C for 4 consecutive days) each winter. I only expect "ev" to affect survival from winter to spring.

Anyway, it sounds like it is doing what it should be, but I expected the DM to look like this, without the second row and that pesky zero:

Code: Select all

> ev.res[[6]]$design.matrix # winter:ev
                      S:(Intercept) S:winter:ev f:(Intercept)
S gfaA0 c1 c1 a0 t1   "1"           "ev1"       "0"         
S gfaA0 c1 c1 a4 t5   "1"           "ev5"       "0"         
S gfaA0 c1 c1 a8 t9   "1"           "ev9"       "0"         
S gfaA0 c1 c1 a12 t13 "1"           "ev13"      "0"         
S gfaA0 c1 c1 a16 t17 "1"           "ev17"      "0"         
S gfaA0 c1 c1 a20 t21 "1"           "ev21"      "0"         
f gfaA0 c1 c1 a0 t1   "0"           "0"         "1"   

Thanks again Jeff!

From what you are describing, you want winter=0/1. If you made it a factor then you couldn't restrict the ev covariate to winter. If I understand what you are doing then you can't get rid of the 0 row because it will be for non-winter periods. If I'm misunderstanding something contact me offlist.

--jeff

The following example may help explain the difference between interactions of numeric and factor variables versus interaction of two factor variables. To make it less confusing I'm not using an individual covariate like in your example.

Create a dataframe with 4 records with each combination of 0:1 and A,B. First w=0,2 is a numeric and then a factor variable. I used 0,2 to make use of numeric clear, otherwise the DM looks the same as a factor variable.

Code: Select all

> df=data.frame(w=rep(c(0,2),2),x=factor(c("A","B","B","A")))
> df
  w x
1 0 A
2 2 B
3 0 B
4 2 A

Full interaction of w and x. First column is intercept (value when w=0), second column is slope for w when x=A, third column is effect of x=B relative to A, and fourth column is difference in slope for w when x=B relative to slope for w when x=A. This is a model in which both slope and intercept can differ for A and B.

Code: Select all

> model.matrix(~w*x,df)
  (Intercept) w xB w:xB
1           1 0  0    0
2           1 2  1    2
3           1 0  1    0
4           1 2  0    0
attr(,"assign")
[1] 0 1 2 3
attr(,"contrasts")
attr(,"contrasts")$x
[1] "contr.treatment"

Interaction only of x and w. Firs column is intercept (value when w=0 for both A and B), the second is the slope for A and the third is the slope for B. These slopes are no longer relative to each other. This would be a linear model with constant intercept but different slopes.

Code: Select all

> model.matrix(~w:x,df)
  (Intercept) w:xA w:xB
1           1    0    0
2           1    0    2
3           1    0    0
4           1    2    0
attr(,"assign")
[1] 0 1 1
attr(,"contrasts")
attr(,"contrasts")$x
[1] "contr.treatment"

Now treat w as a factor variable with 2 levels of 0 and 2. When you interact 2 factor variables with 2 levels each you should get 2*2=4 parameters - one for each of the combinations. We no longer have slopes we simply get means - for each combination.

Below is the full interaction model. First column is the overall intercept, second column is the effect of w=2 value (no longer twice the effect), third column is the effect of B, and fourth column is the effect of B for w=2.

Code: Select all

> model.matrix(~as.factor(w)*x,df)
  (Intercept) as.factor(w)2 xB as.factor(w)2:xB
1           1             0  0                0
2           1             1  1                1
3           1             0  1                0
4           1             1  0                0
attr(,"assign")
[1] 0 1 2 3
attr(,"contrasts")
attr(,"contrasts")$`as.factor(w)`
[1] "contr.treatment"

attr(,"contrasts")$x
[1] "contr.treatment"


Now if you did the interaction only it gives you 5 columns which is too many parameters. The intercept is not needed because all you want is the 4 unique values which gives an identity design matrix if you rearrange columns 2-5:

Code: Select all

> model.matrix(~as.factor(w):x,df)
  (Intercept) as.factor(w)0:xA as.factor(w)2:xA as.factor(w)0:xB as.factor(w)2:xB
1           1                1                0                0                0
2           1                0                0                0                1
3           1                0                0                1                0
4           1                0                1                0                0
attr(,"assign")
[1] 0 1 1 1 1
attr(,"contrasts")
attr(,"contrasts")$`as.factor(w)`
[1] "contr.treatment"

attr(,"contrasts")$x
[1] "contr.treatment"


To remove the intercept for this type of DM, just include -1 which removes the intercept.

Code: Select all

> model.matrix(~-1+as.factor(w):x,df)
  as.factor(w)0:xA as.factor(w)2:xA as.factor(w)0:xB as.factor(w)2:xB
1                1                0                0                0
2                0                0                0                1
3                0                0                1                0
4                0                1                0                0
attr(,"assign")
[1] 1 1 1 1
attr(,"contrasts")
attr(,"contrasts")$`as.factor(w)`
[1] "contr.treatment"

attr(,"contrasts")$x
[1] "contr.treatment"

>
$x
[1] "contr.treatment"

>

The confusing aspect here is that if your numeric variable is 0/1 like winter, the numeric variable is being used simply to restrict the effect of an individual covariate to particular times. You only use the -1 when you are interacting two factor variables. And as it explains in the documentation that approach will not work if you have something like ~-1+y+x:w where x,y and w are factor variables. In that case you need to use the remove.intercept argument. To convince yourself of that create a dataframe with 3 factor variables and try model.matrix with the formula.

Thank you Jeff!