/
fit1_pvl.Rd
61 lines (55 loc) · 2.16 KB
/
fit1_pvl.Rd
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
% Generated by roxygen2: do not edit by hand
% Please edit documentation in R/fit1_pvl.R
\name{fit1_pvl}
\alias{fit1_pvl}
\title{Fit a model for a specified d-tuple of markers}
\usage{
fit1_pvl(indices, start_snp, probs, addcovar, inv_S, S, pheno)
}
\arguments{
\item{indices}{a vector of indices for extracting elements of `probs` array}
\item{start_snp}{an integer to specify the index of the marker where the scan - in call to scan_pvl - starts. This argument is needed because `mytab` has only relative indices (relative to the `start_snp` marker)}
\item{probs}{founder allele probabilities matrix}
\item{addcovar}{additive covariates matrix}
\item{inv_S}{inverse covariance matrix for the vectorized phenotype}
\item{S}{covariance matrix for the vectorized phenotype, ie, the inverse of inv_S}
\item{pheno}{a n by d phenotypes matrix}
}
\value{
a number, the log-likelihood for the specified model
}
\description{
`fit1_pvl` uses several functions in the package qtl2pleio to fit the
linear mixed effects model for a single d-tuple of markers.
Creation of `fit1_pvl` - from code that originally resided in `scan_pvl`, enabled
parallelization via the `parallel` R package.
}
\examples{
# read data
iron <- qtl2::read_cross2(system.file("extdata", "iron.zip", package="qtl2"))
# insert pseudomarkers into map
map <- qtl2::insert_pseudomarkers(iron$gmap, step=1)
# calculate genotype probabilities
probs <- qtl2::calc_genoprob(iron, map, error_prob=0.002)
# grab phenotypes and covariates; ensure that covariates have names attribute
pheno <- iron$pheno
# leave-one-chromosome-out kinship matrices
kinship <- qtl2::calc_kinship(probs, "loco")
# get founder allele probabilites
probs <- qtl2::genoprob_to_alleleprob(probs)
cc_out <- calc_covs(pheno, kinship$`1`, max_iter = 10000, max_prec = 1e-07, covariates = NULL)
Vg <- cc_out$Vg
Ve <- cc_out$Ve
# define Sigma
Sigma <- calc_Sigma(Vg, Ve, kinship$`1`)
# define Sigma_inv
Sigma_inv <- solve(Sigma)
# prepare table of marker indices for each call of scan_pvl
mytab <- prep_mytab(d_size = 2, n_snp = 10)
# set up parallel analysis
fit1_pvl(mytab[1, ], start_snp = 1,
probs = probs$`1`, addcovar = NULL, inv_S = Sigma_inv,
S = Sigma,
pheno = pheno
)
}