RepeatABEL (version 1.1)

preFitModel: Fits a linear mixed model (without fixed SNP effects) and computes the fitted variance-covariance matrix for later use in the rGLS function.

Description

Uses a GenABEL object and phenotype data as input. The model is fitted using the hglm function in the hglm package.

Usage

preFitModel(fixed = y ~ 1, random = ~1 | id, id.name = "id", genabel.data, phenotype.data, corStruc = NULL, GRM = NULL, Neighbor.Matrix = NULL)

Arguments

fixed

A formula including the response and fixed effects

random

A formula for the random effects

id.name

The column name of the IDs in phen.data

genabel.data

An GenABEL object including marker information. This object has one observation per individual.

phenotype.data

A data frame including the repeated observations and IDs.

corStruc

A list specifying the correlation structure for each random effect. The options are: "Ind" for iid random effects, "GRM" for a correlation structure given by a genetic relationship matrix, or "CAR" for a spatial correlation structure given by a Conditional Autoregressive model specified by a neighborhood matrix.

GRM

A genetic relationship matrix. If not specified whilst the "GRM" option is given for corStruc then the GRM is computed internally within the function.

Neighbor.Matrix

A neighborhood matrix having non-zero value for an element (i,j) where the observations i and j come from neighboring locations. The diagonal elements should be zero.

Value

Returns a list including the fitted hglm object fitted.hglm, the variance-covariance matrix V and the ratios between estimated variance components for the random effects divided by the residual variance, ratio.

Examples

  
####### FIRST EXAMPLE USING GRM #############
 data(Phen.Data) #Phenotype data with repeated observations
 data(gen.data) #GenABEL object including IDs and marker genotypes
 GWAS1 <- rGLS(y ~ age + sex, genabel.data = gen.data, phenotype.data = Phen.Data)
 plot(GWAS1, main="")
 summary(GWAS1)
 #Summary for variance component estimation without SNP effects
 summary(GWAS1@call$hglm)
 #The same results can be computed using the preFitModel as follows
 fixed = y ~ age + sex
 Mod1 <- preFitModel(fixed, random=~1|id, genabel.data = gen.data,
                     phenotype.data = Phen.Data, corStruc=list( id=list("GRM","Ind") ))
 GWAS1b <- rGLS(fixed, genabel.data = gen.data,
                phenotype.data = Phen.Data, V = Mod1$V)
 plot(GWAS1b, main="Results using the preFitModel function")
 ####### SECOND EXAMPLE USING CAR #############
 # Add a fake nest variable to the data just to run the example
 #In this example there are 6 nests and 60 observations per nest
 Phen.Data$nest <- rep(1:6, each=60)
 #A model including polygenic effects, permanent environmental effects,
 #and nest effect as random
 Mod2 <- preFitModel(fixed, random=~1|id + 1|nest, genabel.data = gen.data,
          phenotype.data = Phen.Data, corStruc=list( id=list("GRM","Ind"), nest=list("Ind")) )
 GWAS2 <- rGLS(fixed, genabel.data = gen.data, phenotype.data = Phen.Data, V = Mod2$V)
 plot(GWAS2)
 #Similar to previous plot because the nest effect variance component is almost 0.
 ###################
 #Construct a fake nighbourhood matrix
 D = matrix(0,6,6)
 D[1,2] = D[2,1] = 1
 D[5,6] = D[6,5] = 1
 D[2,4] = D[4,2] = 1
 D[3,5] = D[5,3] = 1
 D[1,6] = D[6,1] = 1
 D[3,4] = D[4,3] = 1
 #The matrix shows which pair of nests that can be considered as neighbours
 image(Matrix(D), main="Neighbourhood matrix")
 Mod3 <- preFitModel(fixed, random=~1|id + 1|nest, genabel.data = gen.data,
          phenotype.data = Phen.Data, corStruc=list( id=list("GRM","Ind"),
                                     nest=list("CAR")), Neighbor.Matrix=D )
 GWAS2b <- rGLS(fixed, genabel.data = gen.data,
                phenotype.data = Phen.Data, V = Mod3$V)
 plot(GWAS2b)