Genomic dominant relationship matrix for the Jersey dataset.

Description

This matrix was calculated by using the dominance incidence matrix derived from 33,267 Single Nucleotide Polymorphisms (SNPs) information on 297 individually cows,

D=\frac{X_d X_d'}{2 \sum_{j=1}^p (p_j^2+q_j^2) p_j q_j},

where

• X_d is the design matrix for allele substitution effects for dominance.

• p_j is the frecuency of the second allele at locus j and q_j=1-p_j.

Source

University of Wisconsin at Madison, USA.

Genomic additive relationship matrix for the Jersey dataset.

Description

A matrix, similar to this was used in Gianola et al. (2011) for predicting milk, fat and protein production in Jersey cows. In this software version we do not center the incidence matrix for the additive effects.

G=\frac{X_a X_a'}{2\sum_{j=1}^p p_j (1-p_j)},

where

• X_a is the design matrix for allele substitution effects for additivity.

• p_j is the frecuency of the second allele at locus j and q_j=1-p_j.

Source

University of Wisconsin at Madison, USA.

References

Gianola, D. Okut, H., Weigel, K. and Rosa, G. 2011. "Predicting complex quantitative traits with Bayesian neural networks: a case study with Jersey cows and wheat". BMC Genetics, 12,87.

Genomic additive relationship matrix for the GLS dataset.

Description

Genomic relationship matrix for the GLS dataset. The matrix was derived from 46347 markers for the 278 individuals. The matrix was calculated as follows G=MM'/p, where M is the matrix of markers centered and standarized and p is the number of markers.

Source

International Maize and Wheat Improvement Center (CIMMYT), Mexico.

References

Montesinos-Lopez, O., A. Montesinos-Lopez, J. Crossa, J. Burgueno and K. Eskridge. 2015. "Genomic-Enabled Prediction of Ordinal Data with Bayesian Logistic Ordinal Regression". G3: Genes | Genomes | Genetics, 5, 2113-2126.

brnn

Description

The brnn function fits a two layer neural network as described in MacKay (1992) and Foresee and Hagan (1997). It uses the Nguyen and Widrow algorithm (1990) to assign initial weights and the Gauss-Newton algorithm to perform the optimization. This function implements the functionality of the function trainbr in Matlab 2010b.

Usage

  brnn(x, ...)
  
  ## S3 method for class 'formula'
brnn(formula, data, contrasts=NULL,...)

  ## Default S3 method:
brnn(x,y,neurons=2,normalize=TRUE,epochs=1000,mu=0.005,mu_dec=0.1, 
       mu_inc=10,mu_max=1e10,min_grad=1e-10,change = 0.001,cores=1,
       verbose=FALSE,Monte_Carlo = FALSE,tol = 1e-06, samples = 40,...)

Arguments

Details

The software fits a two layer network as described in MacKay (1992) and Foresee and Hagan (1997). The model is given by:

y_i=g(\boldsymbol{x}_i)+e_i = \sum_{k=1}^s w_k g_k (b_k + \sum_{j=1}^p x_{ij} \beta_j^{[k]}) + e_i, i=1,...,n

where:

• e_i \sim N(0,\sigma_e^2).

• s is the number of neurons.

• w_k is the weight of the k-th neuron, k=1,...,s.

• b_k is a bias for the k-th neuron, k=1,...,s.

• \beta_j^{[k]} is the weight of the j-th input to the net, j=1,...,p.

• g_k(\cdot) is the activation function, in this implementation g_k(x)=\frac{\exp(2x)-1}{\exp(2x)+1}.

The software will minimize

F=\beta E_D + \alpha E_W

where

• E_D=\sum_{i=1}^n (y_i-\hat y_i)^2, i.e. the error sum of squares.

• E_W is the sum of squares of network parameters (weights and biases).

• \beta=\frac{1}{2\sigma^2_e}.

• \alpha=\frac{1}{2\sigma_\theta^2}, \sigma_\theta^2 is a dispersion parameter for weights and biases.

Value

object of class "brnn" or "brnn.formula". Mostly internal structure, but it is a list containing:

References

Bai, Z. J., M. Fahey and G. Golub. 1996. "Some large-scale matrix computation problems." Journal of Computational and Applied Mathematics 74(1-2), 71-89.

Foresee, F. D., and M. T. Hagan. 1997. "Gauss-Newton approximation to Bayesian regularization", Proceedings of the 1997 International Joint Conference on Neural Networks.

Gianola, D. Okut, H., Weigel, K. and Rosa, G. 2011. "Predicting complex quantitative traits with Bayesian neural networks: a case study with Jersey cows and wheat". BMC Genetics, 12,87.

MacKay, D. J. C. 1992. "Bayesian interpolation", Neural Computation, 4(3), 415-447.

Nguyen, D. and Widrow, B. 1990. "Improving the learning speed of 2-layer neural networks by choosing initial values of the adaptive weights", Proceedings of the IJCNN, 3, 21-26.

Paciorek, C. J. and Schervish, M. J. 2004. "Nonstationary Covariance Functions for Gaussian Process Regression". In Thrun, S., Saul, L., and Scholkopf, B., editors, Advances in Neural Information Processing Systems 16. MIT Press, Cambridge, MA.

See Also

predict.brnn

Examples

## Not run: 

#Load the library
library(brnn)

###############################################################
#Example 1 
#Noise triangle wave function, similar to example 1 in Foresee and Hagan (1997)

#Generating the data
x1=seq(0,0.23,length.out=25)
y1=4*x1+rnorm(25,sd=0.1)
x2=seq(0.25,0.75,length.out=50)
y2=2-4*x2+rnorm(50,sd=0.1)
x3=seq(0.77,1,length.out=25)
y3=4*x3-4+rnorm(25,sd=0.1)
x=c(x1,x2,x3)
y=c(y1,y2,y3)


#With the formula interface
out=brnn(y~x,neurons=2)

#With the default S3 method the call is
#out=brnn(y=y,x=as.matrix(x),neurons=2)

plot(x,y,xlim=c(0,1),ylim=c(-1.5,1.5),
     main="Bayesian Regularization for ANN 1-2-1")
lines(x,predict(out),col="blue",lty=2)
legend("topright",legend="Fitted model",col="blue",lty=2,bty="n")

###############################################################
#Example 2
#sin wave function, example in the Matlab 2010b demo.

x = seq(-1,0.5,length.out=100)
y = sin(2*pi*x)+rnorm(length(x),sd=0.1)

#With the formula interface
out=brnn(y~x,neurons=3)

#With the default method the call is
#out=brnn(y=y,x=as.matrix(x),neurons=3)

plot(x,y)
lines(x,predict(out),col="blue",lty=2)

legend("bottomright",legend="Fitted model",col="blue",lty=2,bty="n")

###############################################################
#Example 3
#2 Inputs and 1 output
#the data used in Paciorek and
#Schervish (2004). The data is from a two input one output function with Gaussian noise
#with mean zero and standard deviation 0.25

data(twoinput)

#Formula interface
out=brnn(y~x1+x2,data=twoinput,neurons=10)

#With the default S3 method
#out=brnn(y=as.vector(twoinput$y),x=as.matrix(cbind(twoinput$x1,twoinput$x2)),neurons=10)

f=function(x1,x2) predict(out,cbind(x1,x2))
x1=seq(min(twoinput$x1),max(twoinput$x1),length.out=50)
x2=seq(min(twoinput$x2),max(twoinput$x2),length.out=50)
z=outer(x1,x2,f) # calculating the density values

transformation_matrix=persp(x1, x2, z,
                            main="Fitted model",
                            sub=expression(y==italic(g)~(bold(x))+e),
                            col="lightgreen",theta=30, phi=20,r=50, 
                            d=0.1,expand=0.5,ltheta=90, lphi=180,
                            shade=0.75, ticktype="detailed",nticks=5)
points(trans3d(twoinput$x1,twoinput$x2, f(twoinput$x1,twoinput$x2), 
               transformation_matrix), col = "red")

###############################################################
#Example 4
#Gianola et al. (2011).
#Warning, it will take a while

#Load the Jersey dataset
data(Jersey)

#Fit the model with the FULL DATA
#Formula interface
out=brnn(pheno$yield_devMilk~G,neurons=2,verbose=TRUE)

#Obtain predictions and plot them against fitted values
plot(pheno$yield_devMilk,predict(out))

#Predictive power of the model using the SECOND set for 10 fold CROSS-VALIDATION
data=pheno
data$X=G
data$partitions=partitions

#Fit the model for the TESTING DATA
out=brnn(yield_devMilk~X,
         data=subset(data,partitions!=2),neurons=2,verbose=TRUE)

#Plot the results
#Predicted vs observed values for the training set
par(mfrow=c(2,1))
plot(out$y,predict(out),xlab=expression(hat(y)),ylab="y")
cor(out$y,predict(out))

#Predicted vs observed values for the testing set
yhat_R_testing=predict(out,newdata=subset(data,partitions==2))
ytesting=pheno$yield_devMilk[partitions==2]
plot(ytesting,yhat_R_testing,xlab=expression(hat(y)),ylab="y")
cor(ytesting,yhat_R_testing)


## End(Not run)

brnn_extended

Description

The brnn_extended function fits a two layer neural network as described in MacKay (1992) and Foresee and Hagan (1997). It uses the Nguyen and Widrow algorithm (1990) to assign initial weights and the Gauss-Newton algorithm to perform the optimization. The hidden layer contains two groups of neurons that allow us to assign different prior distributions for two groups of input variables.

Usage

  brnn_extended(x, ...)

  ## S3 method for class 'formula'
brnn_extended(formula, data, contrastsx=NULL,contrastsz=NULL,...)

  ## Default S3 method:
brnn_extended(x,y,z,neurons1,neurons2,normalize=TRUE,epochs=1000,
              mu=0.005,mu_dec=0.1, mu_inc=10,mu_max=1e10,min_grad=1e-10,
              change = 0.001, cores=1,verbose =FALSE,...)

Arguments

Details

The software fits a two layer network as described in MacKay (1992) and Foresee and Hagan (1997). The model is given by:

y_i= \sum_{k=1}^{s_1} w_k^{1} g_k (b_k^{1} + \sum_{j=1}^p x_{ij} \beta_j^{1[k]}) + \sum_{k=1}^{s_2} w_k^{2} g_k (b_k^{2} + \sum_{j=1}^q z_{ij} \beta_j^{2[k]})\,\,e_i, i=1,...,n

• e_i \sim N(0,\sigma_e^2).

• g_k(\cdot) is the activation function, in this implementation g_k(x)=\frac{\exp(2x)-1}{\exp(2x)+1}.

The software will minimize

F=\beta E_D + \alpha \theta_1' \theta_1 +\delta \theta_2' \theta_2

where

• E_D=\sum_{i=1}^n (y_i-\hat y_i)^2, i.e. the sum of squared errors.

• \beta=\frac{1}{2\sigma^2_e}.

• \alpha=\frac{1}{2\sigma_{\theta_1}^2}, \sigma_{\theta_1}^2 is a dispersion parameter for weights and biases for the associated to the first group of neurons.

• \delta=\frac{1}{2\sigma_{\theta_2}^2}, \sigma_{\theta_2}^2 is a dispersion parameter for weights and biases for the associated to the second group of neurons.

Value

object of class "brnn_extended" or "brnn_extended.formula". Mostly internal structure, but it is a list containing:

References

Foresee, F. D., and M. T. Hagan. 1997. "Gauss-Newton approximation to Bayesian regularization", Proceedings of the 1997 International Joint Conference on Neural Networks.

MacKay, D. J. C. 1992. "Bayesian interpolation", Neural Computation, 4(3), 415-447.

Nguyen, D. and Widrow, B. 1990. "Improving the learning speed of 2-layer neural networks by choosing initial values of the adaptive weights", Proceedings of the IJCNN, 3, 21-26.

See Also

predict.brnn_extended

Examples

## Not run: 

#Example 5
#Warning, it will take a while

#Load the Jersey dataset
data(Jersey)

#Predictive power of the model using the SECOND set for 10 fold CROSS-VALIDATION
data=pheno
data$G=G
data$D=D
data$partitions=partitions

#Fit the model for the TESTING DATA for Additive + Dominant
out=brnn_extended(yield_devMilk ~ G | D,
                                  data=subset(data,partitions!=2),
                                  neurons1=2,neurons2=2,epochs=100,verbose=TRUE)

#Plot the results
#Predicted vs observed values for the training set
par(mfrow=c(2,1))
yhat_R_training=predict(out)
plot(out$y,yhat_R_training,xlab=expression(hat(y)),ylab="y")
cor(out$y,yhat_R_training)

#Predicted vs observed values for the testing set
newdata=subset(data,partitions==2,select=c(D,G))
ytesting=pheno$yield_devMilk[partitions==2]
yhat_R_testing=predict(out,newdata=newdata)
plot(ytesting,yhat_R_testing,xlab=expression(hat(y)),ylab="y")
cor(ytesting,yhat_R_testing)
  

## End(Not run)
 

brnn_ordinal

Description

The brnn_ordinal function fits a Bayesian Regularized Neural Network for Ordinal data.

Usage

  brnn_ordinal(x, ...)
  
  ## S3 method for class 'formula'
brnn_ordinal(formula, data, contrasts=NULL,...)
  
  ## Default S3 method:
brnn_ordinal(x,
               y,
               neurons=2,
               normalize=TRUE,
               epochs=1000,
               mu=0.005,
               mu_dec=0.1,
               mu_inc=10,
               mu_max=1e10,
               min_grad=1e-10,
               change_F=0.01,
               change_par=0.01,
               iter_EM=1000,
               verbose=FALSE,
               ...)

Arguments

Details

The software fits a Bayesian Regularized Neural Network for Ordinal data. The model is an extension of the two layer network as described in MacKay (1992); Foresee and Hagan (1997), and Gianola et al. (2011). We use the latent variable approach described in Albert and Chib (1993) to model ordinal data, the Expectation maximization (EM) and Levenberg-Marquardt algorithm (Levenberg, 1944; Marquardt, 1963) to fit the model.

Following Albert and Chib (1993), suppose that Y_1,...,Y_n are observed and Y_i can take values on L ordered values. We are interested in modelling the probability p_{ij}=P(Y_i=j) using the covariates x_{i1},...,x_{ip}. Let

g(\boldsymbol{x}_i)= \sum_{k=1}^s w_k g_k (b_k + \sum_{j=1}^p x_{ij} \beta_j^{[k]}),

where:

• s is the number of neurons.

• w_k is the weight of the k-th neuron, k=1,...,s.

• b_k is a bias for the k-th neuron, k=1,...,s.

• \beta_j^{[k]} is the weight of the j-th input to the net, j=1,...,p.

• g_k(\cdot) is the activation function, in this implementation g_k(x)=\frac{\exp(2x)-1}{\exp(2x)+1}.

Let

Z_i=g(\boldsymbol{x}_i)+e_i,

where:

• e_i \sim N(0,1).

• Z_i is an unobserved (latent variable).

The output from the model for latent variable is related to observed data using the approach employed in the probit and logit ordered models, that is Y_i=j if \lambda_{j-1}<Z_i<\lambda_{j}, where \lambda_j are a set of unknown thresholds. We assign prior distributions to all unknown quantities (see Albert and Chib, 1993; Gianola et al., 2011) for further details. The Expectation maximization (EM) and Levenberg-Marquardt algorithm (Levenberg, 1944; Marquardt, 1963) to fit the model.

Value

object of class "brnn_ordinal". Mostly internal structure, but it is a list containing:

References

Albert J, and S. Chib. 1993. Bayesian Analysis of Binary and Polychotomus Response Data. JASA, 88, 669-679.

Foresee, F. D., and M. T. Hagan. 1997. "Gauss-Newton approximation to Bayesian regularization", Proceedings of the 1997 International Joint Conference on Neural Networks.

Gianola, D. Okut, H., Weigel, K. and Rosa, G. 2011. "Predicting complex quantitative traits with Bayesian neural networks: a case study with Jersey cows and wheat". BMC Genetics, 12,87.

Levenberg, K. 1944. "A method for the solution of certain problems in least squares", Quart. Applied Math., 2, 164-168.

MacKay, D. J. C. 1992. "Bayesian interpolation", Neural Computation, 4(3), 415-447.

Marquardt, D. W. 1963. "An algorithm for least-squares estimation of non-linear parameters". SIAM Journal on Applied Mathematics, 11(2), 431-441.

See Also

predict.brnn_ordinal

Examples

## Not run: 
#Load the library
library(brnn)

#Load the dataset
data(GLS)

#Subset of data for location Harare
HarareOrd=subset(phenoOrd,Loc=="Harare")

#Eigen value decomposition for GOrdm keep those 
#eigen vectors whose corresponding eigen-vectors are bigger than 1e-10
#and then compute principal components

evd=eigen(GOrd)
evd$vectors=evd$vectors[,evd$value>1e-10]
evd$values=evd$values[evd$values>1e-10]
PC=evd$vectors%*%sqrt(diag(evd$values))
rownames(PC)=rownames(GOrd)

#Response variable
y=phenoOrd$rating
gid=as.character(phenoOrd$Stock)

Z=model.matrix(~gid-1)
colnames(Z)=gsub("gid","",colnames(Z))

if(any(colnames(Z)!=rownames(PC))) stop("Ordering problem\n")

#Matrix of predictors for Neural net
X=Z%*%PC

#Cross-validation
set.seed(1)
testing=sample(1:length(y),size=as.integer(0.10*length(y)),replace=FALSE)
isNa=(1:length(y)%in%testing)
yTrain=y[!isNa]
XTrain=X[!isNa,]
nTest=sum(isNa)

neurons=2
	
fmOrd=brnn_ordinal(XTrain,yTrain,neurons=neurons,verbose=FALSE)

#Predictions for testing set
XTest=X[isNa,]
predictions=predict(fmOrd,XTest)
predictions



## End(Not run)

estimate.trace

Description

The estimate.trace function estimates the trace of the inverse of a possitive definite and symmetric matrix using the algorithm developed by Bai et al. (1996). It is specially useful when the matrix is huge.

Usage

  estimate.trace(A,tol=1E-6,samples=40,cores=1)

Arguments

References

Bai, Z. J., M. Fahey and G. Golub. 1996. "Some large-scale matrix computation problems." Journal of Computational and Applied Mathematics, 74(1-2), 71-89.

Examples

## Not run: 
library(brnn)
data(Jersey)

#Estimate the trace of the inverse of G matrix
estimate.trace(G)

#The TRUE value
sum(diag(solve(G)))


## End(Not run)

Initialize networks weights and biases

Description

Function to initialize the weights and biases in a neural network. It uses the Nguyen-Widrow (1990) algorithm.

Usage

     initnw(neurons,p,n,npar)

Arguments

Details

The algorithm is described in Nguyen-Widrow (1990) and in other books, see for example Sivanandam and Sumathi (2005). The algorithm is briefly described below.

• 1.-Compute the scaling factor \theta=0.7 p^{1/n}.

• 2.- Initialize the weight and biases for each neuron at random, for example generating random numbers from U(-0.5,0.5).

• 3.- For each neuron:

  • compute \eta_k=\sqrt{\sum_{j=1}^p (\beta_j^{(k)})^2},

  • update (\beta_1^{(k)},...,\beta_p^{(k)})',

    \beta_j^{(k)}=\frac{\theta \beta_j^{(k)}}{\eta_k}, j=1,...,p,

  • Update the bias (b_k) generating a random number from U(-\theta,\theta).

Value

A list containing initial values for weights and biases. The first s components of the list contains vectors with the initial values for the weights and biases of the k-th neuron, i.e. (\omega_k, b_k, \beta_1^{(k)},...,\beta_p^{(k)})'.

References

Nguyen, D. and Widrow, B. 1990. "Improving the learning speed of 2-layer neural networks by choosing initial values of the adaptive weights", Proceedings of the IJCNN, 3, 21-26.

Sivanandam, S.N. and Sumathi, S. 2005. Introduction to Neural Networks Using MATLAB 6.0. Ed. McGraw Hill, First edition.

Examples

## Not run: 
#Load the library
library(brnn)

#Set parameters
neurons=3
p=4
n=10
npar=neurons*(1+1+p)+1
initnw(neurons=neurons,p=p,n=n,npar=npar)


## End(Not run)

Jacobian

Description

Internal function for the calculation of the Jacobian.

normalize

Description

Internal function for normalizing the data. This function makes a linear transformation of the inputs such that the values lie between -1 and 1.

Usage

   normalize(x,base,spread)

Arguments

Details

z=2*(x-base)/spread - 1

Value

A vector or matrix with the resulting normalized values.

Partitions for cross validation (CV)

Description

Is a vector (297 \times 1) that assigns observations to 10 disjoint sets; the assignment was generated at random. This is used later to conduct a 10-fold CV.

Source

University of Wisconsin at Madison, USA.

Phenotypic information for Jersey

Description

The format of the phenotype dataframe is animal ID, herd, number of lactations, average milk yield, average fat yield, average protein yield, yield deviation for milk, yield deviation for fat, and yield deviation for protein. Averages are adjusted for age, days in milk, and milking frequency. Yield deviations are adjusted further, for herd-year-season effect. You may wish to use yield deviations, because there are relatively few cows per herd (farmers don't pay to genotype all of their cows, just the good ones).

Source

University of Wisconsin at Madison, USA.

References

Gianola, D. Okut, H., Weigel, K. and Rosa, G. 2011. "Predicting complex quantitative traits with Bayesian neural networks: a case study with Jersey cows and wheat". BMC Genetics, 12,87.

Phenotypic information for GLS (ordinal trait)

Description

The format of the phenotype dataframe is Location (Loc), Replicate (Rep), Genotype ID (Stock) and rating (ordinal score). Gray Leaf Spot (GLS) is a disease caused by the fungus Cercospora zeae-maydis. This dataset consists of genotypic and phenotypic information for 278 maize lines from the Drought Tolerance Maize (DTMA) project of CIMMYT's Global Maize Program. The dataset includes information on disease severity measured on an ordinal scale with 5 points: 1= no disease, 2= low infection, 3=moderate infection, 4=high infection and 5=totally infected.

Source

International Maize and Wheat Improvement Center (CIMMYT), Mexico.

References

Montesinos-Lopez, O., A. Montesinos-Lopez, J. Crossa, J. Burgueno and K. Eskridge. 2015. "Genomic-Enabled Prediction of Ordinal Data with Bayesian Logistic Ordinal Regression". G3: Genes | Genomes | Genetics, 5, 2113-2126.

predict.brnn

Description

The function produces the predictions for a two-layer feed-forward neural network.

Usage

   ## S3 method for class 'brnn'
predict(object,newdata,...)

Arguments

Details

This function is a method for the generic function predict() for class "brnn". It can be invoked by calling predict(x) for an object x of the appropriate class, or directly by calling predict.brnn(x) regardless of the class of the object.

Value

A vector containing the predictions

predict.brnn_extended

Description

The function produces the predictions for a two-layer feed-forward neural network.

Usage

   ## S3 method for class 'brnn_extended'
predict(object,newdata,...)

Arguments

Details

This function is a method for the generic function predict() for class "brnn_extended". It can be invoked by calling predict(x) for an object x of the appropriate class, or directly by calling predict.brnn(x) regardless of the class of the object.

Value

A vector containig the predictions

predict.brnn_ordinal

Description

The function produces the predictions for a two-layer feed-forward neural network for ordinal data.

Usage

   ## S3 method for class 'brnn_ordinal'
predict(object,newdata,...)

Arguments

Details

This function is a method for the generic function predict() for class "brnn_ordinal". It can be invoked by calling predict(x) for an object x of the appropriate class, or directly by calling predict.brnn_ordinal(x) regardless of the class of the object.

Value

A list with components:

2 Inputs and 1 output.

Description

The data used in Paciorek and Schervish (2004). This is a data.frame with 3 columns, columns 1 and 2 corresponds to the predictors and column 3 corresponds to the target.

References

Paciorek, C. J. and Schervish, M. J. 2004. "Nonstationary Covariance Functions for Gaussian Process Regression". In Thrun, S., Saul, L., and Scholkopf, B., editors, Advances in Neural Information Processing Systems 16. MIT Press, Cambridge, MA.

un_normalize

Description

Internal function for going back to the original scale.

Usage

   un_normalize(z,base,spread)

Arguments

Details

x=base+0.5*spread*(z+1)

Value

A vector or matrix with the resulting un normalized values.