处理 sommer 包中的不同性别
Dealing with separate sexes in the sommer package
我正在尝试弄清楚如何在 R 的 sommer
包中指定特定的混合模型。
我测量了一群男性个体的两个特征,以及一群女性亲属的两个特征。目的是估计这 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
的工作方式与 MCMCglmm
中的 at.level
不同,因为我在这里使用它似乎什么都不做。
我是不是漏掉了一些关键的东西,或者 sommer
目前缺少这个功能?
非常感谢!
为了在 sommer 中处理不同的性别,你需要安装 sommer >=3.7。该模型直接操作 Gtc 参数(在 vs 函数中),该参数处理应用于方差-协方差分量的约束。值 1、2、3 对应于 positive、unconstrained、fixed。
鉴于此,数据的结构清楚地表明所有性状之间的协方差只能针对 "genotype" 随机效应进行估计:
> 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
而 "individual" 不是一个完全非结构化的特征,而是看起来像这样:
> 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)
您可能希望使用真实数据来获得一些有意义的结果,但这说明了方法。干杯。
我正在尝试弄清楚如何在 R 的 sommer
包中指定特定的混合模型。
我测量了一群男性个体的两个特征,以及一群女性亲属的两个特征。目的是估计这 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
的工作方式与 MCMCglmm
中的 at.level
不同,因为我在这里使用它似乎什么都不做。
我是不是漏掉了一些关键的东西,或者 sommer
目前缺少这个功能?
非常感谢!
为了在 sommer 中处理不同的性别,你需要安装 sommer >=3.7。该模型直接操作 Gtc 参数(在 vs 函数中),该参数处理应用于方差-协方差分量的约束。值 1、2、3 对应于 positive、unconstrained、fixed。
鉴于此,数据的结构清楚地表明所有性状之间的协方差只能针对 "genotype" 随机效应进行估计:
> 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
而 "individual" 不是一个完全非结构化的特征,而是看起来像这样:
> 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)
您可能希望使用真实数据来获得一些有意义的结果,但这说明了方法。干杯。