Работа с отдельными полами в летнем пакете - PullRequest
Новая
       6

Работа с отдельными полами в летнем пакете

0 голосов
/ 05 ноября 2018

Я пытаюсь выяснить, как указать конкретную смешанную модель в пакете sommer для R.

Я измерил две черты на группе мужчин и две черты на их родственниках. Цель оценить генетическую (со) дисперсию внутри и между этими 4 признаками. Однако, поскольку люди не могут быть обоего пола, я хочу подогнать модель так, чтобы она оценивала остаточную ковариацию между чертой женщины 1 и чертой женщины 2, а также остаточную ковариацию между чертой мужчины 1 и чертой 2 мужчины, но НЕ остаточные ковариации между любыми мужскими и женскими чертами (поскольку в данных не должно быть никакой информации об этих ковариациях). В MCMCglmm это можно сделать с помощью кода, подобного следующему (предполагается, что две женские черты находятся в столбцах 1 и 2, а две мужские черты - в столбцах 3 и 4 матрицы переменных ответа): 

rcov = ~us(at.level(trait, 1:2)):units + us(at.level(trait, 3:4)):units

Но в sommer, похоже, нет эквивалентной функции: я получаю сообщение об ошибке 

Error: On the meantime the only rcov structures available are:
     'rcov=~units' or 'rcov=~at(.):units'.

Затем я попытался назначить каждому человеку номер и подобрать случайный эффект индивидуального уровня, например: 

library(sommer)

# Generate some fake data: 
# 100 males and 100 females
# Two traits are measured on each male, and two traits on each female
# 20 individuals per sex are measured for each of 5 different genotypes 
df <- data.frame(
  sex = rep(c("female", "male"), each = 100),
  female_trait_1 = c(rnorm(100), rep(NA, 100)),
  female_trait_2 = c(rnorm(100), rep(NA, 100)),
  male_trait_1 = c(rep(NA, 100), rnorm(100)),
  male_trait_2 = c(rep(NA, 100), rnorm(100)),
  genotype = rep(rep(1:5, each = 20), 2),
  individual = 1:200
)
df$genotype <- as.factor(df$genotype)
df$individual <- as.factor(df$individual)

sommer_test <- mmer2(

  # four traits as multivariate response
  cbind(female_trait_1, 
        female_trait_2,
        male_trait_1, 
        male_trait_2) ~ 1,

  # Fit the random effect of genotype (to estimate genetic covariance within and between sexes)
  # Try to fit US covariance matrices to specific levels of 'trait' (does not work)
  random =~ 
    us(trait):genotype +
    us(at.levels(trait, c("female_trait_1", "female_trait_2"))):individual +
    us(at.levels(trait, c("male_trait_1", "male_trait_2"))):individual,

  data = df
)
summary(sommer_test)

Однако, последний также не работает - он работает, но матрица США для индивидуума только что подобрана дважды, для всех комбинаций черт (включая комбинации мужчины и женщины, которых я пытался избежать, используя at.levels). Таким образом, похоже, что at.levels не работает так же, как at.level в MCMCglmm, поскольку, похоже, ничего не делает, поскольку я использовал его здесь.

Я что-то упустил, или эта функция отсутствует в sommer в настоящее время?

Большое спасибо!

1 Ответ

0 голосов
/ 06 ноября 2018

Для того, чтобы иметь дело с отдельными полами в летнее время, нужно установить летнее время> = 3.7. Эта модель напрямую манипулирует аргументом Gtc (в функции vs), который обрабатывает ограничения, применяемые для компонентов дисперсии и ковариации. Значения 1,2,3 соответствуют положительным, неограниченным, фиксированным.

Учитывая сказанное, структура данных проясняет, что ковариация среди всех признаков может быть оценена только для случайного эффекта "генотип": 

> unsm(4)
     [,1] [,2] [,3] [,4]
[1,]    1    2    2    2
[2,]    0    1    2    2
[3,]    0    0    1    2
[4,]    0    0    0    1

, тогда как для «индивида» это не полная неструктурированная черта, а выглядит так: 

> mm <- adiag1(unsm(2),unsm(2));mm
     [,1] [,2] [,3] [,4]
[1,]    1    2    0    0
[2,]    0    1    0    0
[3,]    0    0    1    2
[4,]    0    0    0    1

Короче говоря, модель выглядит следующим образом: 

# Generate some fake data: 
# 100 males and 100 females
# Two traits are measured on each male, and two traits on each female
# 20 individuals per sex are measured for each of 5 different genotypes 
set.seed(3434)
df <- data.frame(
  sex = rep(c("female", "male"), each = 100),
  female_trait_1 = c(rnorm(100), rep(NA, 100)),
  female_trait_2 = c(rnorm(100), rep(NA, 100)),
  male_trait_1 = c(rep(NA, 100), rnorm(100)),
  male_trait_2 = c(rep(NA, 100), rnorm(100)),
  genotype = rep(rep(1:5, each = 20), 2),
  individual = 1:200
)
df$genotype <- as.factor(df$genotype)
df$individual <- as.factor(df$individual)

library(sommer)
mm <- adiag1(unsm(2),unsm(2));mm
mix <- mmer(cbind(female_trait_1, 
                  female_trait_2,
                  male_trait_1,
                  male_trait_2) ~ 1,
     random=~vs(genotype,Gtc=unsm(4)) + vs(individual,Gtc=mm),
     rcov=~vs(units), na.method.Y = "include",
     data=df)

mix$sigma
cov2cor(mix$sigma$genotype)
cov2cor(mix$sigma$individual)

Возможно, вы захотите использовать реальные данные, чтобы получить значимые результаты, но это показывает путь. Приветствия.

...