data fatcs1;              format type $3.;

infile 'c:\Documents and Settings\uzivatel\Dokumenty\genetickeparametry2003\fatcs1\meanvariances2.fat' 
lrecl=160 missover;

INPUT
v1
V2;
poradi=_n_;
if poradi<79 then type='G';
if poradi>156 then type='REZ';
if poradi ge 79 and poradi le 156 then type='PE';

run;


title ' ';


************************************************;
*ADDITIVe Genetic variances pro celou laktaci;
************************************************;




***********************************************;
*PRIPRAVA DAT;
***********************************************;

  data GE;
  set  fatcs1;
  if type='G '; RUN;


  data PE;
  set  fatcs1;
  if type='PE '; RUN;

  data RE;
  set  fatcs1;
  if type='REZ';



Proc IML;
*matice - vstup;

***********************************************;

*ADDITIVE Genetic (co)variance matrix;

***********************************************;


                             /* begin IML session */

use ge;


v=j(1);

G=j(12);
i=0;



   DO j=1 to 12 ;                     /* begin do loop */
   DO k=j to 12 ;
        read next into v;
      if j = k then do;

             G[j,j]=V[1,1];
             end;
      else do;
              G[j,k]=V[1,1];
              G[k,j]=V[1,1];

       end;


   end;
   end;
*print G;


***********************************************;

*PERMANENT environmental (co)variance matrix;

***********************************************;
use PE;
v=j(1);
P=j(12);
i=0;
   DO j=1 to 12 ;                     /* begin do loop */
   DO k=j to 12 ;
        read next into v;
      if j = k then do;

             P[j,j]=V[1,1];
             end;
      else do;
              P[j,k]=V[1,1];
              P[k,j]=V[1,1];

       end;

   end;
   end;


*print P;


***********************************************;

*RESIUUM variance matrix;

***********************************************;

                             /* begin IML session */

use RE;


v=j(1);

R=j(12,1);
i=0;



   DO j=1 to 12 ;                     /* begin do loop */
        read next into v;
             R[j]=V[1,1];
   END;


*print R;



************************************************;
*ADDITIVe Genetic variances pro celou laktaci;
************************************************;

celkem={0,0,0,0,0,0,0,0,0,0,0,0};
zero={0,0,0,0};

   *PRVNI LAKTACE;

create odhady1 var{var1};

   DO i=6 to 305;                     /* begin do loop */

     age=(2*(i-6)*inv(299))-1;

      x1= sqrt(2)*sqrt(0.50);

      x2= sqrt(2)*sqrt(1.50)*age ;

      x3= sqrt(2)*sqrt(2.50)*(1.5*age**2-0.5);

      x4= sqrt(2)*sqrt(3.50)*(2.5*age**3-1.5*age) ;


     z=x1//x2//x3//x4//zero//zero;
celkem=celkem+z;


end;                                    /* end do loop */

*prvni laktace;

var1=celkem`*G*celkem;

print var1;

append ;

********************;
*DRUHA LAKTACE;
********************;

create odhady2 var{var2};

celkem={0,0,0,0,0,0,0,0,0,0,0,0};


     DO i=6 to 305;                     /* begin do loop */

     age=(2*(i-6)*inv(299))-1;

      x1= sqrt(2)*sqrt(0.50);

      x2= sqrt(2)*sqrt(1.50)*age ;

      x3= sqrt(2)*sqrt(2.50)*(1.5*age**2-0.5);

      x4= sqrt(2)*sqrt(3.50)*(2.5*age**3-1.5*age) ;


     z=zero//x1//x2//x3//x4//zero;
celkem=celkem+z;
end;

var2=celkem`*G*celkem;

print var2;

append ;


********************;
*TRETI LAKTACE;
********************;

create odhady3 var{var3};
   ;
celkem={0,0,0,0,0,0,0,0,0,0,0,0};



DO i=6 to 305;                     /* begin do loop */

     age=(2*(i-6)*inv(299))-1;

      x1= sqrt(2)*sqrt(0.50);

      x2= sqrt(2)*sqrt(1.50)*age ;

      x3= sqrt(2)*sqrt(2.50)*(1.5*age**2-0.5);

      x4= sqrt(2)*sqrt(3.50)*(2.5*age**3-1.5*age) ;


     z=zero//zero//x1//x2//x3//x4;
celkem=celkem+ z;
end;



   var3=celkem`*G*celkem;
        *print var3;

append ;

************************************************;
*PE variances pro celou laktaci;
************************************************;

celkem={0,0,0,0,0,0,0,0,0,0,0,0};
zero={0,0,0,0};


**********************;
   *PRVNI LAKTACE;
**********************;

create perhady1 var{per1,re1};

   DO i=6 to 305;                     /* begin do loop */

     age=(2*(i-6)*inv(299))-1;

      x1= sqrt(2)*sqrt(0.50);

      x2= sqrt(2)*sqrt(1.50)*age ;

      x3= sqrt(2)*sqrt(2.50)*(1.5*age**2-0.5);

      x4= sqrt(2)*sqrt(3.50)*(2.5*age**3-1.5*age) ;


     z=x1//x2//x3//x4//zero//zero;
celkem=celkem+z;


end;                                    /* end do loop */

*prvni laktace;

per1=celkem`*P*celkem;

*print per1;

***********************************;
Vypocet residua;
************************************;

RM=Block(r[1],r[2],r[3],r[4]);

*print RM;

z={40,70,150,60};


RE1=z`*RM*z;


*print RE1;

append ;

********************;
*DRUHA LAKTACE;
********************;

create perhady2 var{per2,re2};

celkem={0,0,0,0,0,0,0,0,0,0,0,0};


     DO i=6 to 305;                     /* begin do loop */

     age=(2*(i-6)*inv(299))-1;

      x1= sqrt(2)*sqrt(0.50);

      x2= sqrt(2)*sqrt(1.50)*age ;

      x3= sqrt(2)*sqrt(2.50)*(1.5*age**2-0.5);

      x4= sqrt(2)*sqrt(3.50)*(2.5*age**3-1.5*age) ;


     z=zero//x1//x2//x3//x4//zero;
celkem=celkem+z;
end;

per2=celkem`*P*celkem;

*print per2;

***********************************;
*Vypocet residua;
************************************;

RM=Block(r[5],r[6],r[7],r[8]);



z={40,70,150,60};

*print RM;

RE2=z`*RM*z;


*print RE2;

append ;


********************;
*TRETI LAKTACE;
********************;

create perhady3 var{per3,re3};
   ;
celkem={0,0,0,0,0,0,0,0,0,0,0,0};



DO i=6 to 305;                     /* begin do loop */

     age=(2*(i-6)*inv(299))-1;

      x1= sqrt(2)*sqrt(0.50);

      x2= sqrt(2)*sqrt(1.50)*age ;

      x3= sqrt(2)*sqrt(2.50)*(1.5*age**2-0.5);

      x4= sqrt(2)*sqrt(3.50)*(2.5*age**3-1.5*age) ;


     z=zero//zero//x1//x2//x3//x4;
celkem=celkem+ z;
end;



   per3=celkem`*P*celkem;
        *print per3;



***********************************;
*Vypocet residua;
************************************;

RM=Block(r[9],r[10],r[11],r[12]);



z={40,70,150,60};

*print RM;

RE3=z`*RM*z;


*print RE3;

append ;


QUIT;


data fat;
  merge odhady1(in=in_1)
        odhady2(in=in_2)
        odhady3(in=in_3)
        perhady1(in=in_1)
        perhady2(in=in_2)
        perhady3(in=in_3);

        ;  run;







data fat;
set fat;

h1=var1/(var1+per1+re1);
h2=var2/(var2+per2+re2);
h3=var3/(var3+per3+re3);run;


proc print data=fat;run;