millironx/beefblup
1
0
Fork 0
mirror of https://github.com/MillironX/beefblup.git synced 2024-12-22 09:08:16 +00:00
Code Issues Releases Packages Activity

Made display prettier

Browse source
This commit is contained in:
Thomas A. Christensen II 2018-11-13 12:37:04 -07:00
parent dfddabde76
commit da5397e4f9
1 changed files with 73 additions and 6 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

77
MATLAB/beefblup.m
View file

@ -10,12 +10,43 @@ clear
clc
close all
% Import data from a suitable spreadsheet
%% Display stuff
disp('beefblup v. 0.0.0.1')
disp('(C) 2018 Thomas A. Christensen II')
disp('https://github.com/millironx/beefblup')
disp(' ')
%% Prompt User
% Ask for an input spreadsheet
[name, path] = uigetfile('*.xlsx','Select a beefblup worksheet');
% Ask for an ouput text file
[savename, savepath, ~] = uiputfile('*.txt', 'Save your beefblup results', 'results');
% Ask for heritability
h2 = input('What is the heritability for this trait? >> ');
%% Import input file
tic
disp(' ')
disp('Importing Excel file...')
% Import data from a suitable spreadsheet
fullname = [path name];
clear name path
[~, ~, data] = xlsread(fullname);
disp('Importing Excel file... Done!')
toc
disp(' ')
%% Process input file
tic
disp(' ')
disp('Processing and formatting data...')
disp(' ')
% Extract the headers into a separate array
headers = data(1,:);
data(1,:) = [];
@ -97,6 +128,15 @@ disp('If no animal matching this description exists, the results may appear')
disp('outlandish, but are still as correct as the accuracy suggests')
disp(' ')
disp('Processing and formatting data... Done!')
toc
disp(' ')
%% Create the fixed-effect matrix
tic
disp(' ')
disp('Creating the fixed-effect matrix...')
% Form the fixed effect matrix
X = zeros(numanimals, sum(numgroups)-length(numgroups)+1);
X(:,1) = ones(1, numanimals);
@ -129,6 +169,15 @@ for i = 1:length(normal)
end
end
disp('Creating the fixed-effect matrix... Done!')
toc
disp(' ')
%% Additive relationship matrix
tic
disp(' ')
disp('Creating the additive relationship matrix...')
% Create an empty matrix for the additive relationship matrix
A = zeros(numanimals, numanimals);
@ -181,6 +230,15 @@ for i = 1:numanimals
end
disp('Creating the additive relationship matrix... Done!')
toc
disp(' ')
%% Perform BLUP
tic
disp(' ')
disp('Solving the mixed-model equations')
% Extract the observed data
Y = cell2mat(data(:, 5));
@ -191,8 +249,7 @@ Z = eye(numanimals, numanimals);
nullobs = find(isnan(Y));
Z(nullobs, nullobs) = 0;
% Prompt for heritability
h2 = input('What is the heritability for this trait? >> ');
% Calculate heritability
lambda = (1-h2)/h2;
% Use the mixed-model equations
@ -202,8 +259,14 @@ solutions = [X'*X X'*Z; Z'*X (Z'*Z)+(inv(A).*lambda)]\[X'*Y; Z'*Y];
diaginv = diag(inv([X'*X X'*Z; Z'*X (Z'*Z)+(inv(A).*lambda)]));
reliability = 1 - diaginv.*lambda;
% Ask the user for where they would like the file saved
[savename, savepath, ~] = uiputfile('*.txt', 'Save your beefblup results', 'results');
disp('Solving the mixed-model equations... Done!')
toc
disp(' ')
%% Output the results
tic
disp(' ')
disp('Saving results...')
% Find how many traits we found BLUE for
numgroups = numgroups - 1;
@ -275,3 +338,7 @@ end
fprintf(fileID, '\n - END REPORT -');
fclose(fileID);
disp('Saving results... Done!')
toc
disp(' ')