Nothing
#----------------------------------------------------------------
# QGglmm R package
# Functions to estimate QG parameters from GLMM estimates
#----------------------------------------------------------------
# Pierre de Villemereuil & Michael Morrissey (2015)
#----------------------------------------------------------------
# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
# You should have received a copy of the GNU General Public License
# along with this program. If not, see <http://www.gnu.org/licenses/>.
##-------------------------- General functions ---------------------- ##
#Calculating the observed/expected scale mean
QGmean <- function(mu = NULL,
var,
link.inv,
predict = NULL,
width = 10)
{
if(length(mu) > 1 | length(var) > 1) {
stop("The parameters mu and var must be of length 1,
please check your input.")
}
#If no fixed effects were included in the model
if (is.null(predict)) {
if(is.null(mu)) {
stop("Please provide either mu or predict.")
} else {
predict <- mu
}
}
mean(
sapply(predict,
function(pred_i) {
integrate(f = function(x) {
link.inv(x) * dnorm(x, pred_i, sqrt(var))
},
lower = pred_i - width * sqrt(var),
upper = pred_i + width * sqrt(var)
)$value
}
)
)
}
#Calculating the expected scale variance
QGvar.exp <- function(mu = NULL,
var,
link.inv,
obs.mean = NULL,
predict = NULL,
width = 10)
{
if(length(mu) > 1 | length(var) > 1) {
stop("The parameters mu and var must be of length 1,
please check your input.")
}
#If no fixed effects were included in the model
if (is.null(predict)) {
if(is.null(mu)) {
stop("Please provide either mu or predict.")
} else {
predict <- mu
}
}
#If not provided, compute the obsereved mean
if (is.null(obs.mean)) {
obs.mean <- QGmean(mu = mu,
var = var,
link.inv = link.inv,
width = width,
predict = predict)
}
#Note: Using Koenig's formula
mean(
sapply(predict,
function(pred_i) {
integrate(f = function(x) {
(link.inv(x)^2) *
dnorm(x, pred_i, sqrt(var))
},
lower = pred_i - width * sqrt(var),
upper = pred_i + width * sqrt(var)
)$value
}
)
) - (obs.mean^2)
}
#Calculating the "distribution" variance
QGvar.dist <- function(mu = NULL,
var,
var.func,
predict = NULL,
width = 10)
{
if(length(mu) > 1 | length(var) > 1) {
stop("The parameters mu and var must be of length 1,
please check your input.")
}
#If no fixed effects were included in the model
if (is.null(predict)) {
if(is.null(mu)) {
stop("Please provide either mu or predict.")
} else {
predict <- mu
}
}
mean(
sapply(predict,
function(pred_i) {
integrate(f = function(x) {
var.func(x) * dnorm(x, pred_i, sqrt(var))
},
lower = pred_i - width * sqrt(var),
upper = pred_i + width * sqrt(var)
)$value
}
)
)
}
#Calculating "psi" for the observed additive genetic variance computation
QGpsi <- function(mu = NULL,
var,
d.link.inv,
predict = NULL,
width = 10)
{
if(length(mu) > 1 | length(var) > 1) {
stop("The parameters mu and var must be of length 1,
please check your input.")
}
#If no fixed effects were included in the model
if (is.null(predict)) {
if(is.null(mu)) {
stop("Please provide either mu or predict.")
} else {
predict <- mu
}
}
mean(
sapply(predict,
function(pred_i) {
integrate(f = function(x) {
d.link.inv(x) * dnorm(x, pred_i, sqrt(var))
},
lower = pred_i - width * sqrt(var),
upper = pred_i + width * sqrt(var)
)$value
}
)
)
}
## ------------------ Dictionary functions ----------------- ##
#Function creating the needed functions according to the link name
QGlink.funcs <- function(name,
n.obs = NULL,
theta = NULL)
{
if (name == "Gaussian") {
inv.link <- function(x) {x}
var.func <- function(x) {0}
d.inv.link <- function(x) {1}
} else if (name == "binom1.probit") {
inv.link <- function(x) {pnorm(x)}
var.func <- function(x) {pnorm(x) * (1 - pnorm(x))}
d.inv.link <- function(x) {dnorm(x)}
} else if (name == "binomN.probit") {
if (is.null(n.obs)) {
stop("binomN model used,
but no observation number (n.obs) defined.")
}
inv.link <- function(x) {n.obs * pnorm(x)}
var.func <- function(x) {n.obs * pnorm(x) * (1 - pnorm(x))}
d.inv.link <- function(x) {n.obs * dnorm(x)}
} else if (name == "binom1.logit") {
inv.link <- function(x) {plogis(x)}
var.func <- function(x) {plogis(x) * (1 - plogis(x))}
d.inv.link <- function(x) {dlogis(x)}
} else if (name == "binomN.logit") {
if (is.null(n.obs)) {
stop("binomN model used,
but no observation number (n.obs) defined.")
}
inv.link <- function(x) {n.obs * plogis(x)}
var.func <- function(x) {n.obs * plogis(x) * (1 - plogis(x))}
d.inv.link <- function(x) {n.obs * dlogis(x)}
} else if (name == "ordinal") {
stop("ordinal models are particular,
please use the QGparams function only")
} else if (name == "Poisson.log") {
inv.link <- function(x) {exp(x)}
var.func <- function(x) {exp(x)}
d.inv.link <- function(x) {exp(x)}
} else if (name == "Poisson.sqrt") {
inv.link <- function(x) {x^2}
var.func <- function(x) {x^2}
d.inv.link <- function(x) {2 * x}
} else if (name == "negbin.log") {
if (is.null(theta)) {stop("negbin model used, but theta not defined.")}
inv.link <- function(x) {exp(x)}
var.func <- function(x) {exp(x) + (exp(2 * x) / theta)}
d.inv.link <- function(x) {exp(x)}
} else if (name == "negbin.sqrt") {
if (is.null(theta)) {
stop("negbin model used, but theta not defined.")
}
inv.link <- function(x) {x^2}
var.func <- function(x) {(x^2) + ((x^4) / theta)}
d.inv.link <- function(x) {2 * x}
} else {
stop("Invalid model name.
Use a valid model name or enter a custom model specification.")
}
list(inv.link = inv.link, var.func = var.func, d.inv.link = d.inv.link)
}
## -----------Special functions for known analytical solutions-------------- ##
qg.Gaussian <- function(mu = NULL,
var.a,
var.p,
predict = NULL)
{
if(length(mu) > 1 | length(var.a) != 1 | length(var.p) != 1) {
stop("The parameters mu and var must be of length 1,
please check your input.")
}
#If no fixed effects were included in the model
if (is.null(predict)) {
if(is.null(mu)) {
stop("Please provide either mu or predict.")
} else {
predict <- mu
}
}
#Nothing to be done, except averaging over predict
if (length(predict) == 1) {
var_fixed <- 0
} else {
var_fixed <- var(predict)
}
data.frame(mean.obs = mean(predict),
var.obs = var.p + var_fixed,
var.a.obs = var.a,
h2.obs = var.a / (var.p + var_fixed))
}
qg.binom1.probit <- function(mu = NULL,
var.a,
var.p,
predict = NULL)
{
if(length(mu) > 1 | length(var.a) != 1 | length(var.p) != 1) {
stop("The parameters mu and var must be of length 1,
please check your input.")
}
#If no fixed effects were included in the model
if (is.null(predict)) {
if(is.null(mu)) {
stop("Please provide either mu or predict.")
} else {
predict <- mu
}
}
#Observed mean
p <- mean(1 - pnorm(0, predict, sqrt(var.p + 1)))
#Observed variance
var_obs <- p * (1 - p)
#Psi
Psi <- mean(dnorm(0, (predict), sqrt(var.p + 1)))
data.frame(mean.obs = p,
var.obs = var_obs,
var.a.obs = (Psi^2) * var.a,
h2.obs = ((Psi^2) * var.a) / var_obs)
}
qg.binomN.probit <- function(mu = NULL,
var.a,
var.p,
n.obs,
predict = NULL,
width = 10)
{
if(length(mu) > 1 | length(var.a) != 1 | length(var.p) != 1) {
stop("The parameters mu and var must be of length 1,
please check your input.")
}
#If no fixed effects were included in the model
if (is.null(predict)) {
if(is.null(mu)) {
stop("Please provide either mu or predict.")
} else {
predict <- mu
}
}
#Observed mean
p <- n.obs * mean(1 - pnorm(0, predict, sqrt(var.p + 1)))
#Observed variance
prob.sq.int <-
mean(
sapply(predict,
function(pred_i) {
integrate(f = function(x) {
(pnorm(x)^2) *
dnorm(x, pred_i, sqrt(var.p))
},
lower = pred_i - width * sqrt(var.p),
upper = pred_i + width * sqrt(var.p)
)$value
}
)
)
var_obs <- ((n.obs^2) - n.obs) * prob.sq.int - p^2 +p
#Psi
Psi <- n.obs * mean(dnorm(0, (predict), sqrt(var.p + 1)))
data.frame(mean.obs = p,
var.obs = var_obs,
var.a.obs = (Psi^2) * var.a,
h2.obs = ((Psi^2) * var.a) / var_obs)
}
qg.Poisson.log <- function(mu = NULL,
var.a,
var.p,
predict = NULL)
{
if(length(mu) > 1 | length(var.a) != 1 | length(var.p) != 1) {
stop("The parameters mu and var must be of length 1,
please check your input.")
}
#If no fixed effects were included in the model
if (is.null(predict)) {
if(is.null(mu)) {
stop("Please provide either mu or predict.")
} else {
predict <- mu
}
}
#Observed mean
lambda <- mean(exp(predict + (var.p / 2)))
#Mean of lambda square, needed for the following
lambda_sq <- mean(exp(2 * (predict + var.p / 2)))
#Observed variance
var_obs <- lambda_sq * exp(var.p) - lambda^2 + lambda
data.frame(mean.obs = lambda,
var.obs = var_obs,
var.a.obs = (lambda^2) * var.a,
h2.obs = ((lambda^2) * var.a) / var_obs)
}
qg.Poisson.sqrt <- function(mu = NULL,
var.a,
var.p,
predict = NULL)
{
if(length(mu) > 1 | length(var.a) != 1 | length(var.p) != 1) {
stop("The parameters mu and var must be of length 1,
please check your input.")
}
#If no fixed effects were included in the model
if (is.null(predict)) {
if(is.null(mu)) {
stop("Please provide either mu or predict.")
} else {
predict <- mu
}
}
#Observed mean
lambda <- mean((predict)^2 + var.p)
#Observed variance
var_obs <- mean((predict)^4 +
6 * var.p * ((predict)^2) +
3 * (var.p^2)) -
lambda^2 + lambda
#Psi
Psi <- mean(2 * (predict))
data.frame(mean.obs = lambda,
var.obs = var_obs,
var.a.obs = (Psi^2) * var.a,
h2.obs = ((Psi^2) * var.a) / var_obs)
}
qg.negbin.log <- function(mu = NULL,
var.a,
var.p,
theta,
predict = NULL)
{
if(length(mu) > 1 | length(var.a) != 1 | length(var.p) != 1) {
stop("The parameters mu and var must be of length 1,
please check your input.")
}
#If no fixed effects were included in the model
if (is.null(predict)) {
if(is.null(mu)) {
stop("Please provide either mu or predict.")
} else {
predict <- mu
}
}
#Observed mean
lambda <- mean(exp(predict + (var.p / 2)))
#Mean of lambda square, needed for the following
lambda_sq <- mean(exp(2 * (predict + var.p / 2)))
#Observed variance
var_obs <- lambda_sq * exp(var.p) -
lambda^2 +
lambda +
mean(exp(2 * (predict + var.p))) / theta
data.frame(mean.obs = lambda,
var.obs = var_obs,
var.a.obs = (lambda^2) * var.a,
h2.obs = ((lambda^2) * var.a) / var_obs)
}
qg.negbin.sqrt <- function(mu = NULL,
var.a,
var.p,
theta,
predict = NULL)
{
if(length(mu) > 1 | length(var.a) != 1 | length(var.p) != 1) {
stop("The parameters mu and var must be of length 1,
please check your input.")
}
#If no fixed effects were included in the model
if (is.null(predict)) {
if(is.null(mu)) {
stop("Please provide either mu or predict.")
} else {
predict <- mu
}
}
#Observed mean
lambda <- mean((predict)^2 + var.p)
#Observed variance
var_obs <- mean((predict)^4 +
6 * var.p * ((predict)^2) +
3 * (var.p^2)) -
lambda^2 + lambda +
(mean((predict)^4 +
6 * var.p * ((predict)^2) +
3 * (var.p^2)) / theta)
#Psi
Psi <- mean(2 * (predict))
data.frame(mean.obs = lambda,
var.obs = var_obs,
var.a.obs = (Psi^2) * var.a,
h2.obs = ((Psi^2) * var.a) / var_obs)
}
##------------------Meta-function for general calculation------------------- ##
QGparams <- function(mu = NULL,
var.a,
var.p,
model = "",
width = 10,
predict = NULL,
closed.form = TRUE,
custom.model = NULL,
n.obs = NULL,
cut.points = NULL,
theta = NULL,
verbose = TRUE)
{
if(length(mu) > 1 | length(var.a) != 1 | length(var.p) != 1) {
stop("The parameters mu and var must be of length 1,
please check your input.")
}
#If no fixed effects were included in the model
if (is.null(predict)) {
if(is.null(mu)) {
stop("Please provide either mu or predict.")
} else {
predict <- mu
}
}
## Using analytical solutions if possible (and asked for: see closed.form)
if (model == "Gaussian" & closed.form) {
#Gaussian model with identity link (e.g. LMM)
if (verbose) {
print("Using the closed forms for a Gaussian model
with identity link (e.g. LMM).")
}
qg.Gaussian(mu = mu,
var.a = var.a,
var.p = var.p,
predict = predict)
} else if (model == "binom1.probit" & closed.form) {
#Binary.probit model
if (verbose) {
print("Using the closed forms for a Binomial1 - probit model.")
}
qg.binom1.probit(mu = mu,
var.a = var.a,
var.p = var.p,
predict = predict)
} else if (model == "threshold") {
#Binary.probit model
if (verbose) {
print("Using the closed forms for a threshold model
(e.g. for MCMCglmm package).
Closed form argument is ignored...")
}
qg.binom1.probit(mu = mu,
var.a = var.a,
var.p = var.p - 1,
predict = predict)
} else if (model == "binomN.probit" & closed.form) {
#Binomial - not - binary model
if (is.null(n.obs)) {
stop("binomN.probit model used,
but no observation number (n.obs) defined.")
}
if (verbose) {
print("Using a semi - closed form for a BinomialN - probit model.")
}
warning("Some component (var.obs) are not totally
from a closed form solution: an integral is computed")
qg.binomN.probit(mu = mu,
var.a = var.a,
var.p = var.p,
predict = predict,
n.obs = n.obs,
width = width)
} else if (model == "Poisson.log" & closed.form){
#Poisson - log model
if(verbose) {
print("Using the closed forms for a Poisson - log model.")
}
qg.Poisson.log(mu = mu,
var.a = var.a,
var.p = var.p,
predict = predict)
} else if (model == "Poisson.sqrt" & closed.form){
#Poisson - sqrt model
if(verbose) {
print("Using the closed forms for a Poisson - sqrt model.")
}
qg.Poisson.sqrt(mu = mu,
var.a = var.a,
var.p = var.p,
predict = predict)
} else if (model == "negbin.log" & closed.form){
#NegBin - log model
if (is.null(theta)) {
stop("negbin model used, but theta not defined.")
}
if(verbose) {
print("Using the closed forms for a NegativeBinomial - log model.")
}
qg.negbin.log(mu = mu,
var.a = var.a,
var.p = var.p,
predict = predict,
theta = theta)
} else if (model == "negbin.sqrt" & closed.form){
#NegBin - sqrt model
if (is.null(theta)) {
stop("negbin model used, but theta not defined.")
}
if(verbose) {
print("Using the closed forms for a NegativeBinomial - sqrt model.")
}
qg.negbin.sqrt(mu = mu,
var.a = var.a,
var.p = var.p,
predict = predict,
theta = theta)
} else if (model == "ordinal"){
#Ordinal model
if (is.null(cut.points)) {
stop("cut points must be specified to use the ordinal model.")
}
if(verbose) {
print("Using the closed forms for an ordinal model
(ignoring the closed.form argument)")
}
qg.ordinal(mu = mu,
var.a = var.a,
var.p = var.p,
predict = predict,
cut.points = cut.points)
} else {
##Else, use the general integral equations
#Use a custom model if defined, otherwise look into the "Dictionary"
if (is.null(custom.model)) {
if (model == "") {
stop("The function requires either model or custom.model.")
} else {
funcs <- QGlink.funcs(model, n.obs = n.obs, theta = theta)
}
} else {
funcs <- custom.model
}
#Observed mean computation
if (verbose) {
print("Computing observed mean...")
}
z_bar <- QGmean(mu = mu,
var = var.p,
link.inv = funcs$inv.link,
width = width,
predict = predict)
#Variances computation
if (verbose) {
print("Computing variances...")
}
var_exp <- QGvar.exp(mu = mu,
var = var.p,
link.inv = funcs$inv.link,
obs.mean = z_bar,
width = width,
predict = predict)
var_dist <- QGvar.dist(mu = mu,
var = var.p,
var.func = funcs$var.func,
width = width,
predict = predict)
var_obs <- var_exp + var_dist
#Psi computation (for the observed additive genetic variance)
if (verbose) {
print("Computing Psi...")
}
Psi <- QGpsi(mu = mu,
var = var.p,
d.link.inv = funcs$d.inv.link,
width = width,
predict = predict)
#Return a data.frame with the calculated components
data.frame(mean.obs = z_bar,
var.obs = var_obs,
var.a.obs = (Psi^2) * var.a,
h2.obs = ((Psi^2) * var.a) / var_obs)
}
}
## -----------Function to calculate the evolutive prediction--------------- ##
QGpred <- function(mu = NULL,
var.a,
var.p,
fit.func,
d.fit.func,
width = 10,
predict = NULL,
verbose = TRUE)
{
if(length(mu) > 1 | length(var.a) != 1 | length(var.p) != 1) {
stop("The parameters mu, var.a and var.p must be of length 1,
please check your input.")
}
if (is.null(predict)) {
if(is.null(mu)) {
stop("Please provide either mu or predict.")
} else {
predict <- mu
}
}
#Calculating the latent mean fitness
if (verbose) {
print("Computing mean fitness...")
}
Wbar <-
mean(
sapply(predict,
function(pred_i) {
integrate(f = function(x) {
fit.func(x) *
dnorm(x, pred_i, sqrt(var.p))
},
lower = pred_i - width * sqrt(var.p),
upper = pred_i + width * sqrt(var.p)
)$value
}
)
)
#Calculating the covariance between latent trait and latent fitness
if (verbose) {
print("Computing the latent selection and response...")
}
#Computing the derivative of fitness
dW <-
mean(
sapply(predict,
function(pred_i) {
integrate(f = function(x) {
d.fit.func(x) *
dnorm(x, pred_i, sqrt(var.p))
},
lower = pred_i - width * sqrt(var.p),
upper = pred_i + width * sqrt(var.p)
)$value
}
)
)
#Computing the selection
if (length(predict) > 1) {
sel <- (var.p + var(predict)) * dW / Wbar
} else {
sel <- var.p * dW / Wbar
}
#Computing the evolutionary response
resp <- var.a * dW / Wbar
#Returning the results on the latent scale
data.frame(mean.lat.fit = Wbar,
lat.grad = dW / Wbar,
lat.sel = sel,
lat.resp = resp)
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.