加载中…用R程序搞定多年多点，一年多点，一点多年数据BLUP
(2012-12-26 11:34:53)
标签：
杂谈 |
分类： 关联分析答疑 |
## Set Working Directory for Mac
setwd("/Users/heathermerk/Documents/TBRT2011")
## Set Working Director for PC
setwd("C:/your_working_directory")
## Read in Brix dataset
qualdat = read.csv("TBRTQuality.csv", header=T)
## Check to ensure data imported correctly
str(qualdat)
head(qualdat)
tail(qualdat)
## Attach dataset
attach(qualdat)
## Examine distribution of brix data
hist(Brix, col="gold")
boxplot(Brix~Loc, xlab="Location", ylab="Degrees Brix", main="Degrees Brix by Location", col="pink")
# Rename variables for ease of use
BRIX = as.numeric(Brix)
LINE = as.factor(Line)
LOC = as.factor(Loc)
YEAR = as.factor(Year)
REP = as.factor(Rep)
## Calculate variance components
# requires lme4 package
library(lme4)
# Linear Model with random effects for variance components
brixvarcomp = lmer(BRIX~ (1|LINE) + (1|LOC) + (1|YEAR) + (1|REP%in%LOC:YEAR) + (1|LINE:LOC) + (1|LINE:YEAR))
# Extract variance components
summary(brixvarcomp) （此处得出的结果可以用来计算遗传力）
## BLUPS（多年多点）
# fit the model
brixmodel = lmer(BRIX ~ (1|LINE) + (1|LOC) + (1|YEAR) + (1|REP%in%LOC:YEAR) + (1|LINE:LOC) + (1|LINE:YEAR))
## BLUPS（一年多点）
# fit the model
brixmodel = lmer(BRIX ~ (1|LINE) + (1|LOC)+ (1|REP%in%LOC) +
(1|LINE:LOC))
setwd("/Users/heathermerk/Documents/TBRT2011")
## Set Working Director for PC
setwd("C:/your_working_directory")
## Read in Brix dataset
qualdat = read.csv("TBRTQuality.csv", header=T)
## Check to ensure data imported correctly
str(qualdat)
head(qualdat)
tail(qualdat)
## Attach dataset
attach(qualdat)
## Examine distribution of brix data
hist(Brix, col="gold")
boxplot(Brix~Loc, xlab="Location", ylab="Degrees Brix", main="Degrees Brix by Location", col="pink")
# Rename variables for ease of use
BRIX = as.numeric(Brix)
LINE = as.factor(Line)
LOC = as.factor(Loc)
YEAR = as.factor(Year)
REP = as.factor(Rep)
## Calculate variance components
# requires lme4 package
library(lme4)
# Linear Model with random effects for variance components
brixvarcomp = lmer(BRIX~ (1|LINE) + (1|LOC) + (1|YEAR) + (1|REP%in%LOC:YEAR) + (1|LINE:LOC) + (1|LINE:YEAR))
# Extract variance components
summary(brixvarcomp) （此处得出的结果可以用来计算遗传力）
## BLUPS（多年多点）
# fit the model
brixmodel = lmer(BRIX ~ (1|LINE) + (1|LOC) + (1|YEAR) + (1|REP%in%LOC:YEAR) + (1|LINE:LOC) + (1|LINE:YEAR))
## BLUPS（一年多点）
# fit the model
brixmodel = lmer(BRIX ~ (1|LINE) + (1|LOC)
##
BLUPS（一点多年）
# fit the model
# fit the model
brixmodel = lmer(BRIX ~ (1|LINE) +
(1|YEAR) + (1|REP%in%YEAR) +
(1|LINE:YEAR))
# estimate BLUPS
brixblup = ranef(brixmodel)
# look at output structure
str(brixblup)
# extract blup for line
brixlineblup = brixblup$LINE
# see the structure of the blup for each line
str(brixlineblup)
# save the brixlineblup output to a separate .csv file
write.csv(brixlineblup, file="BrixLineBLUPS.csv")
## Creating plots with the BLUPs
# Create a numeric vector with the BLUP for each line
LINEBLUP = brixlineblup[,1]
# Create a histogram with the BLUP for each line
hist(LINEBLUP, col="brown")
## Compare BLUP to line averages on a scatterplot
lmean = tapply(BRIX, LINE, na.rm=T, mean)
plot(LINEBLUP, lmean, col="blue")
详细操作见，http://www.extension.org/pages/61006/the-solcap-tomato-phenotypic-data:-estimating-heritability-and-blups-for-traits
brixblup = ranef(brixmodel)
# look at output structure
str(brixblup)
# extract blup for line
brixlineblup = brixblup$LINE
# see the structure of the blup for each line
str(brixlineblup)
# save the brixlineblup output to a separate .csv file
write.csv(brixlineblup, file="BrixLineBLUPS.csv")
## Creating plots with the BLUPs
# Create a numeric vector with the BLUP for each line
LINEBLUP = brixlineblup[,1]
# Create a histogram with the BLUP for each line
hist(LINEBLUP, col="brown")
## Compare BLUP to line averages on a scatterplot
lmean = tapply(BRIX, LINE, na.rm=T, mean)
plot(LINEBLUP, lmean, col="blue")
详细操作见，http://www.extension.org/pages/61006/the-solcap-tomato-phenotypic-data:-estimating-heritability-and-blups-for-traits
更新2017-3-13
This error occurs because the model is overparameterized;
there is only one observation for every value of “REP%in%LOC:YEAR”.
That means that estimating both var(residual) and
var(REP%in%LOC:YEAR) can lead to spurious results (i.e. estimated
solutions that represent a good fit, but are not necessarily
reliable). There are at least two ways to enable estimation of
variance components in this situation:
The problem can be avoided by altering the model and removing
(1|REP%in%LOC:YEAR)) so that the model is not
overparamterized.
y = lmer(BRIX~ (1|LINE) + (1|LOC) + (1|YEAR) + (1|LINE:LOC) +
(1|LINE:YEAR))
喜欢
0
赠金笔
前一篇：答疑——显性标记数据格式
后一篇：关联分析答疑-SSR标记处理