spaMM (version 4.5.0)

X.GCA: Fixed-effect terms for dyadic interactions

Description

X.GCA and X.antisym are functions which, when called in a model formula, stand for terms designed to represent the effet of symmetric interactions between pairs of individuals (order of individuals in the pair does not matter). antisym likewise represents anti-symmetric interactions (the effect of reciprocal ordered pairs on the outcome are opposite, as in the so-called Bradley-Terry models). These constructs all account for multiple membership, i.e., the fact that the same individual may act as the first or the second individual among different pairs, or even within one pair if this makes sense).

The outcome of an interaction between a pair \(i,j\) of agents is subject to a symmetric overall effect \(a_{ij}\) when the effect “on” individual \(i\) (or viewed from the perspective of individual \(i\)) equals the effect on \(j\): \(a_{ij}=a_{ji}\). This may result from the additive effect of individual effects \(a_{i}\) and \(a_{j}\): \(a_{ij}=a_i+a_j\). A X.GCA call represents such symmetric additive effects. Conversely, antisymmetry is characterized by \(a_{ij}=a_i-a_j=-a_{ji}\) and is represented by a X.antisym call. See the diallel documentation for similar constructs for random effects, for additional comments on semantics (e.g. about “GCA”), and for further references.

If individual-level factors ID1 + ID2 were included in a formula for dyadic interactions, this would result in different coefficients being fitted for the same level in each factor. By contrast, the present constucts ensure that a single coefficient is fitted for the same-named levels of factors ID1 and ID2.

Usage

X.GCA(term, contr="contr.treatment", ...) 
X.antisym(term, contr="contr.treatment", ...)

Value

The functions return design matrices with additional class "factorS" and attributes "call" and "spec_levs".

Arguments

term

an expression of the form <.>:<.> where each <.> represents a factor (or a variable that will automaticall be converted to a factor) present in the data of the fitting function call.

contr

The contrasts used. Only the default and "contr.sum" are implemented.

...

For programming purposes, not documented.

Details

The fixed-effect terms (GCA(Par1,Par2), etc) from the lmDiallel package (Onofri & Terzaroli, 2021), defined by functions returning design matrices for usage with stats:lm, work in a formula for a spaMM fit, and users can define use such functions as templates for additional functions that will work similarly. However, not all post-fit functions will handle terms defined in this way well: checking robustness of predict on small and permuted newdata, as shown in the Examples, is a good test. Such problems happen because the formula-handling machinery of R handles terms represented by either a matrix or a factor, while both the model matrix and the factor information used to construct dyadic-interaction terms are needed to correctly predict, with new data, from models including such terms.

The presently designed functions are defined to solve this issue. By using such functions as template, users can define additional functions with the same return format (as further explained in the documented source code of X.antisym), which will allow them to perform correct predictions from fitted models.

References

Onofri A., Terzaroli N. (2021) lmDiallel: Linear fixed effects models for diallel crosses. Version 0.9.4. https://cran.r-project.org/package=lmDiallel.

Examples

#### Simulate dyadic data

set.seed(123)
nind <- 10       # Beware data grow as O(nind^2)
x <- runif(nind^2) 
id12 <- expand.grid(id1=seq(nind),id2=seq(nind))
id1 <- id12$id1
id2 <- id12$id2
u <-  rnorm(nind,mean = 0, sd=0.5)

## additive individual effects:
y <-  0.1 + 1*x + u[id1] +  u[id2] + rnorm(nind^2,sd=0.2)

## anti-smmetric individual effects:
t <-  0.1 + 1*x + u[id1] - u[id2] + rnorm(nind^2,sd=0.2)

dyaddf <- data.frame(x=x, y=y, t=t, id1=id1,id2=id2)
# : note that this contains two rows per dyad, which avoids identifiability issues.

# Enforce that interactions are between distinct individuals (not essential for the fit):
dyaddf <- dyaddf[- seq.int(1L,nind^2,nind+1L),] 


# Fits:

(addfit <- fitme(y ~x +X.GCA(id1:id2), data=dyaddf))

(antifit <- fitme(t ~x +X.antisym(id1:id2), data=dyaddf))

if (FALSE) { #### check of correct handling by predict():

  # First scramble the data so that input factors are in no particular order
  set.seed(123)
  dyaddf <- dyaddf[sample(nrow(dyaddf)),]
  
  addfiti <- fitme(y ~x +X.GCA(id1:id2), data=dyaddf)
  foo <- rev(2:4)
  p1 <- predict(addfiti)[foo]
  p2 <- predict(addfiti, newdata=dyaddf[foo,])
  diff(range(p1-p2))<1e-10    # must be TRUE
}