data milkcs1;              format type $3.;

infile 'c:\Documents and Settings\uzivatel\Dokumenty\genetickeparametry2003\mlkcs1/meanvariances2.milk'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 PE;
  set  milkcs1;
  if type='PE '; RUN;

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



Proc IML;
*matice - vstup;

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

*ADDITIVE Genetic (co)variance matrix;

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


                             /* begin IML session */

use milkcs1_B;


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};
celkem1={0,0,0,0,0,0,0,0,0,0,0,0};
zero={0,0,0,0};


*****************************;
   *PRVNI x druha LAKTACE;
*****************************;


create odhady12 var{cov12,var1,var2};

   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) ;

  *    xv(1,n)=const*dsqrt(0.5d0)
      xv(2,n)=const*dsqrt(1.5d0)*age
      xv(3,n)=const*dsqrt(2.5d0)*(1.5*age**2-0.5)
      xv(4,n)=const*dsqrt(3.5d0)*(2.5*age**3-1.5*age);




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



end;                                    /* end do loop */

*prvni laktace;

cov12=celkem`*G*celkem1;
var1=celkem`*G*celkem;
var2=celkem1`*G*celkem1;


append ;


*****************************;
   *PRVNI x treti LAKTACE;
*****************************;
celkem={0,0,0,0,0,0,0,0,0,0,0,0};
celkem1={0,0,0,0,0,0,0,0,0,0,0,0};


create odhady13 var{cov13,var3};

   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;
     y=zero//zero//x1//x2//x3//x4;
celkem1=celkem1+y;



end;                                    /* end do loop */

*prvni laktace;

cov13=celkem`*G*celkem1;

var3=celkem1`*G*celkem1;


append ;

    *****************************;
   *DRUHA x TRETI LAKTACE;
*****************************;
celkem={0,0,0,0,0,0,0,0,0,0,0,0};
celkem1={0,0,0,0,0,0,0,0,0,0,0,0};


create odhady23 var{cov23};

   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;
     y=zero//zero//x1//x2//x3//x4;
celkem1=celkem1+y;



end;                                    /* end do loop */



cov23=celkem`*G*celkem1;




append ;


QUIT; *IML;







****************************************;
*ADITIVNI GENETICKE korelace mezi laktacemi;
*************************************;





data corrmilk_b;
  merge odhady12(in=in_1)
        odhady13(in=in_2)
        odhady23(in=in_3)
        ;
rg23=cov23/sqrt(var3*var2);
rg12=cov12/sqrt(var1*var2);
rg13=cov13/sqrt(var1*var3);
run;

proc print;run;