Statistical Procedures

Programming the statistical procedures from SAS
BookmarkSubscribeRSS Feed
☑ This topic is solved. Need further help from the community? Please sign in and ask a new question.
Abdus-Salaam
Obsidian | Level 7

Dear SAS-Community,
i want to calculate coefficients of variance (CV) for grain dry matter data with 12 variants and 4 replicates. I want then to test differences of CV-Values between the treatments and the fertiliser levels. Here is a table containing the variants with the corresponding treatment and fertiliser level:

Variant

treatment

fertiliser level

1

1

2

2

3

3

4

4

5

10 

6

11 

7

12 

7

13 

8

14 

8


When i run my code i get the following warning in the Log:


The final Hessian matrix is not positive definite, and therefore the estimated covariance matrix is not full rank and may be unreliable. The variance of some parameter estimates is zero or some parameters are linearly related to other parameters.

 

Could someone of you help me with finding a solution to this problem?

Here comes my code with sample data (similar to my own data):

data Total_Grain_DM_Sample;
input Block Variant Total_Grain_DM_kg_ha;
datalines;
2	1	4016.79
1	1	3904.58
3	1	3857.5
3	1	3962.11
4	2	3777.43
2	2	3857.63
1	2	4053.94
3	2	3902.49
4	3	4494.53
3	3	4607.35
1	3	4616.06
2	3	4660.09
2	4	4580.31
4	4	5225.37
1	4	4666.99
4	4	4782.53
1	5	5271.92
4	5	5302.86
4	5	5056.13
2	5	5033.14
2	6	5213.04
3	6	5477.23
4	6	5173.18
1	6	5482.44
3	7	5389.63
2	7	5482.13
1	7	5764.17
4	7	4824.41
1	8	5662.74
3	8	5906.35
3	8	5904.86
2	8	5693.09
2	9	5821.93
3	9	5989.53
1	9	6351.82
4	9	5813.85
4	10	6210.08
1	10	6257.52
3	10	6125.38
2	10	5959
1	11	6225.87
2	11	6236.02
3	11	6328.76
4	11	6395.25
3	12	6677.88
4	12	6447.28
1	12	6274.82
2	12	6341.02
;run;
;

/* Dummy-Variables for NLMixed */
data Total_Grain_DM_Sample;set Total_Grain_DM_Sample;
if Variant=1 then trt1=1;else trt1=0;
trt2=0;if Variant=2 then trt2=1;
trt3=0;if Variant=3 then trt3=1;
trt4=0;if Variant=4 then trt4=1;
trt5=0;if Variant=5 then trt5=1;
trt6=0;if Variant=6 then trt6=1;
trt7=0;if Variant=7 then trt7=1;
trt8=0;if Variant=8 then trt8=1;
trt9=0;if Variant=9 then trt9=1;
trt10=0;if Variant=10 then trt10=1;
trt11=0;if Variant=11 then trt11=1;
trt12=0;if Variant=12 then trt12=1;
run;
;

/* Berechnung und Vergleich der Variationskoeffizienten */
Proc nlmixed data=Total_Grain_DM_Sample;
parameters t1=3200 t2=4500.00 t3=3700.00 t4=5800 t5=4200.00 t6=4000.00 t7=3250.00 t8=4600.99 t9=4500.00 t10=3200.00 t11=3800.00 t12=4800.00 v1=26000.00 v2=27000.00 v3=262133.76 v4=2120000.00 v5=275000 v6=140500.00 v7=261000.00 v8=280000.00 v9=262500.00 v10=250000.00 v11=108000.00 v12=111000.00;

mu=trt1*t1+trt2*t2+trt3*t3+trt4*t4+trt5*t5+trt6*t6+trt7*t7+trt8*t8+trt9*t9+trt10*t10+trt11*t11+trt12*t12+u1*trt1+u2*trt2+u3*trt3+u4*trt4+u5*trt5+u6*trt6+u7*trt7+u8*trt8+u9*trt9+u10*trt10+u11*trt11+u12*trt12;

model Total_Grain_DM_kg_ha ~ normal (mu,0.0000001);
random u1 u2 u3 u4 u5 u6 u7 u8 u9 u10 u11 u12 ~ normal([0,0,0,0,0,0,0,0,0,0,0,0],
  [v1,
   0, v2,
   0, 0, v3,
   0, 0, 0, v4,
   0, 0, 0, 0, v5,
   0, 0, 0, 0, 0, v6,
   0, 0, 0, 0, 0, 0, v7,
   0, 0, 0, 0, 0, 0, 0, v8,
   0, 0, 0, 0, 0, 0, 0, 0, v9,
   0, 0, 0, 0, 0, 0, 0, 0, 0, v10,
   0, 0, 0, 0, 0, 0, 0, 0, 0, 0, v11,
   0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, v12])   SUBJECT = block OUT = ausgabe;

/* Calculation of CV-values */
predict ((v1)**(1/2)/(t1))*100 out=cv_t1;
predict ((v2)**(1/2)/(t2))*100 out=cv_t2;
predict ((v3)**(1/2)/(t3))*100 out=cv_t3;
predict ((v4)**(1/2)/(t4))*100 out=cv_t4;
predict ((v5)**(1/2)/(t5))*100 out=cv_t5;
predict ((v6)**(1/2)/(t6))*100 out=cv_t6;
predict ((v7)**(1/2)/(t7))*100 out=cv_t7;
predict ((v8)**(1/2)/(t8))*100 out=cv_t8;
predict ((v9)**(1/2)/(t9))*100 out=cv_t9;
predict ((v10)**(1/2)/(t10))*100 out=cv_t10;
predict ((v11)**(1/2)/(t11))*100 out=cv_t11;
predict ((v12)**(1/2)/(t12))*100 out=cv_t12;

/* Comparison of the CV-values */

/* Between the treatments */
/* with fertiliser level 0 */
predict ((v1)**(1/2)/(t1))*100-((v3)**(1/2)/(t3))*100 out=cv_t1_vs_cv_t3;
predict ((v5)**(1/2)/(t5))*100-((v1)**(1/2)/(t1))*100 out=cv_t5_vs_cv_t1;
predict ((v5)**(1/2)/(t5))*100-((v3)**(1/2)/(t3))*100 out=cv_t5_vs_cv_t3;
predict ((v5)**(1/2)/(t5))*100-((v7)**(1/2)/(t7))*100 out=cv_t5_vs_cv_t7;
predict ((v5)**(1/2)/(t5))*100-((v9)**(1/2)/(t9))*100 out=cv_t5_vs_cv_t9;
predict ((v5)**(1/2)/(t5))*100-((v10)**(1/2)/(t10))*100 out=cv_t5_vs_cv_t10;
predict ((v5)**(1/2)/(t5))*100-((v11)**(1/2)/(t11))*100 out=cv_t5_vs_cv_t11;
predict ((v7)**(1/2)/(t7))*100-((v1)**(1/2)/(t1))*100 out=cv_t7_vs_cv_t1;
predict ((v7)**(1/2)/(t7))*100-((v9)**(1/2)/(t9))*100 out=cv_t7_vs_cv_t9;
predict ((v7)**(1/2)/(t7))*100-((v10)**(1/2)/(t10))*100 out=cv_t7_vs_cv_t10;
predict ((v7)**(1/2)/(t7))*100-((v11)**(1/2)/(t11))*100 out=cv_t7_vs_cv_t11;
predict ((v9)**(1/2)/(t9))*100-((v1)**(1/2)/(t1))*100 out=cv_t9_vs_cv_t1;
predict ((v9)**(1/2)/(t9))*100-((v3)**(1/2)/(t3))*100 out=cv_t9_vs_cv_t3;
predict ((v9)**(1/2)/(t9))*100-((v10)**(1/2)/(t10))*100 out=cv_t9_vs_cv_t10;
predict ((v9)**(1/2)/(t9))*100-((v11)**(1/2)/(t11))*100 out=cv_t9_vs_cv_t11;
predict ((v10)**(1/2)/(t10))*100-((v1)**(1/2)/(t1))*100 out=cv_t10_vs_cv_t1;
predict ((v10)**(1/2)/(t10))*100-((v3)**(1/2)/(t3))*100 out=cv_t10_vs_cv_t3;
predict ((v10)**(1/2)/(t10))*100-((v11)**(1/2)/(t11))*100 out=cv_t10_vs_cv_t11;
/* with fertiliser 1evel 1 */
predict ((v2)**(1/2)/(t2))*100-((v4)**(1/2)/(t4))*100 out=cv_t2_vs_cv_t4;
predict ((v6)**(1/2)/(t6))*100-((v2)**(1/2)/(t2))*100 out=cv_t6_vs_cv_t2;
predict ((v6)**(1/2)/(t6))*100-((v4)**(1/2)/(t4))*100 out=cv_t6_vs_cv_t4;
predict ((v6)**(1/2)/(t6))*100-((v8)**(1/2)/(t8))*100 out=cv_t6_vs_cv_t8;
predict ((v6)**(1/2)/(t6))*100-((v12)**(1/2)/(t12))*100 out=cv_t6_vs_cv_t12;
predict ((v8)**(1/2)/(t8))*100-((v2)**(1/2)/(t2))*100 out=cv_t8_vs_cv_t2;
predict ((v8)**(1/2)/(t8))*100-((v4)**(1/2)/(t4))*100 out=cv_t8_vs_cv_t4;
predict ((v8)**(1/2)/(t8))*100-((v12)**(1/2)/(t12))*100 out=cv_t8_vs_cv_t12;
predict ((v12)**(1/2)/(t12))*100-((v2)**(1/2)/(t2))*100 out=cv_t12_vs_cv_t2;
predict ((v12)**(1/2)/(t12))*100-((v4)**(1/2)/(t4))*100 out=cv_t12_vs_cv_t4;

/* Between the fertiliser levels */
predict ((v2)**(1/2)/(t2))*100-((v1)**(1/2)/(t1))*100 out=cv_t2_vs_cv_t1;
predict ((v4)**(1/2)/(t4))*100-((v3)**(1/2)/(t3))*100 out=cv_t4_vs_cv_t3;
predict ((v6)**(1/2)/(t6))*100-((v5)**(1/2)/(t5))*100 out=cv_t6_vs_cv_t5;
predict ((v8)**(1/2)/(t8))*100-((v7)**(1/2)/(t7))*100 out=cv_t8_vs_cv_t7;
predict ((v12)**(1/2)/(t12))*100-((v11)**(1/2)/(t11))*100 out=cv_t12_vs_cv_t11;
run;
;



 

 

1 ACCEPTED SOLUTION

Accepted Solutions
Abdus-Salaam
Obsidian | Level 7

My supervisor has found a simple solution for my problem:

First he has created new start for the parameters of the model in proc NLmixed with proc mixed (means of the variants and the corresponding (standard errors)²) without the repeated /group=block value that i've had in the proc mixed code for creating start values for nlmixed. That he did since the block effect in my model was set 0, what i haven't recognised. Additionally he divided the grain dry matter values by 1000 wich doesn*t effect the CV since it's calculated with a quotient with two values that are equally effected by the transformation of the dry matter values.

Thank you again Mr. Denham and KSharp for your efforts to help me 

View solution in original post

7 REPLIES 7
SteveDenham
Jade | Level 19

This behavior may be due to misspecifying the "block" variable for Variants=1 and 8 (Block 3 appears twice, with no Block 4) and Variants=4 and 5 (Block 4 appears twice, with no Block 3). Another complicating factor is the imbalance in the data, where treatments 5 and 6 only were measured at fertilizer level = 0.  Combining these may make your mapping to the trtXX variable such that you have overlap where there are zero estimates or confounded estimates.

 

Have you tried fitting a model with only Variants 1 and 2? That may shed some light on what is going on.

 

SteveDenham

Abdus-Salaam
Obsidian | Level 7

Hello Mr. Denham,

thank you for your answer. I've tried it again with just variants 1 and 2 with the same warning in the log as before.

data Total_Grain_DM_Sample;
input Block Variant Total_Grain_DM_kg_ha;
datalines;
2	1	4016.79
1	1	3904.58
3	1	3857.5
3	1	3962.11
4	2	3777.43
2	2	3857.63
1	2	4053.94
3	2	3902.49
4	3	4494.53
3	3	4607.35
1	3	4616.06
2	3	4660.09
2	4	4580.31
4	4	5225.37
1	4	4666.99
4	4	4782.53
1	5	5271.92
4	5	5302.86
4	5	5056.13
2	5	5033.14
2	6	5213.04
3	6	5477.23
4	6	5173.18
1	6	5482.44
3	7	5389.63
2	7	5482.13
1	7	5764.17
4	7	4824.41
1	8	5662.74
3	8	5906.35
3	8	5904.86
2	8	5693.09
2	9	5821.93
3	9	5989.53
1	9	6351.82
4	9	5813.85
4	10	6210.08
1	10	6257.52
3	10	6125.38
2	10	5959
1	11	6225.87
2	11	6236.02
3	11	6328.76
4	11	6395.25
3	12	6677.88
4	12	6447.28
1	12	6274.82
2	12	6341.02
;run;
;

/* Dummy-Variables for NLMixed */
data Total_Grain_DM_Sample;set Total_Grain_DM_Sample;
if Variant=1 then trt1=1;else trt1=0;
trt2=0;if Variant=2 then trt2=1;
run;
;

/* Berechnung und Vergleich der Variationskoeffizienten */
Proc nlmixed data=Total_Grain_DM_Sample;
parameters t1=3200 t2=4500.00 t3=3700.00 t4=5800 t5=4200.00 t6=4000.00 t7=3250.00 t8=4600.99 t9=4500.00 t10=3200.00 t11=3800.00 t12=4800.00 v1=26000.00 v2=27000.00 v3=262133.76 v4=2120000.00 v5=275000 v6=140500.00 v7=261000.00 v8=280000.00 v9=262500.00 v10=250000.00 v11=108000.00 v12=111000.00;

mu=trt1*t1+trt2*t2+u1*trt1+u2*trt2;

model Total_Grain_DM_kg_ha ~ normal (mu,0.0000001);
random u1 u2 ~ normal([0,0],
  [v1,
   0, v2])   SUBJECT = block OUT = ausgabe;

/* Calculation of CV-values */
predict ((v1)**(1/2)/(t1))*100 out=cv_t1;
predict ((v2)**(1/2)/(t2))*100 out=cv_t2;

/* Comparison of the CV-values */

/* Between the fertiliser levels */
predict ((v2)**(1/2)/(t2))*100-((v1)**(1/2)/(t1))*100 out=cv_t2_vs_cv_t1;
run;
;
Abdus-Salaam
Obsidian | Level 7

Hello KSharp,

thank you for the document. It helped me to understand the warning better.

Abdus-Salaam
Obsidian | Level 7

I'm sorry, i've forgotten to mention an additional note that appears in the log and that may be relevant:


NOTE: Moore-Penrose inverse is used in covariance matrix.

Abdus-Salaam
Obsidian | Level 7

I think the covariance matrix of the parameter estimates from the original data is helpful to understand the problem. I guess the model is overparameterised, or what do you think ?

The issue is that i need all the variants in the model to compare their CV's.

covariance matrix of the parameter estimates
  t1 t2 t3 t4 t5 t6 t7 t8 t9 t10 t11 t12 v1 v2 v3 v4 v5 v6 v7 v8 v9 v10 v11 v12
t1 88381 0 -5.14E-07 0 6.28E-07 -1.60E-07 0 -1.34E-07 -1.29E-07 0 4.10E-07 2.64E-07 35.243 . . 0 0 . . 0 . . 0 .
t2 0 68045 -1.63E-07 -1.96E-07 4.11E-07 -2.38E-07 1.92E-07 -3.53E-07 1.42E-07 -1.59E-07 2.30E-07 -1.31E-07 0 . . 0 0 . . 0 . . 0 .
t3 -5.14E-07 -1.63E-07 65533 -1.82E-06 0 -2.20E-07 -2.08E-07 -1.34E-07 -3.16E-07 -2.10E-07 -1.12E-07 0 0 . . 0 0 . . 0 . . 0 .
t4 0 -1.96E-07 -1.82E-06 708844 -4.06E-07 -2.56E-06 -2.50E-07 -9.31E-07 -5.56E-07 3.82E-07 0 1.11E-06 8.51E-08 . . 65.6625 8.12E-08 . . 0 . . 2.09E-07 .
t5 6.28E-07 4.11E-07 0 -4.06E-07 91854 0 2.59E-07 -6.47E-07 6.06E-07 1.91E-07 3.69E-07 2.12E-07 0 . . 0 19.1913 . . 0 . . 0 .
t6 -1.60E-07 -2.38E-07 -2.20E-07 -2.56E-06 0 35216 -3.98E-07 1.74E-07 0 0 1.20E-07 -8.15E-08 0 . . 0 0 . . 0 . . 0 .
t7 0 1.92E-07 -2.08E-07 -2.50E-07 2.59E-07 -3.98E-07 65382 0 1.39E-07 0 0 9.81E-08 0 . . 0 0 . . 0 . . 0 .
t8 -1.34E-07 -3.53E-07 -1.34E-07 -9.31E-07 -6.47E-07 1.74E-07 0 94267 -3.08E-07 -1.16E-07 0 2.21E-07 0 . . 0 0 . . 74.6835 . . 0 .
t9 -1.29E-07 1.42E-07 -3.16E-07 -5.56E-07 6.06E-07 0 1.39E-07 -3.08E-07 65556 0 2.76E-07 1.11E-07 0 . . 0 0 . . 0 . . 0 .
t10 0 -1.59E-07 -2.10E-07 3.82E-07 1.91E-07 0 0 -1.16E-07 0 62535 0 1.01E-07 0 . . 0 0 . . 0 . . 0 .
t11 4.10E-07 2.30E-07 -1.12E-07 0 3.69E-07 1.20E-07 0 0 2.76E-07 0 36036 0 0 . . 0 0 . . 0 . . 31.9055 .
t12 2.64E-07 -1.31E-07 0 1.11E-06 2.12E-07 -8.15E-08 9.81E-08 2.21E-07 1.11E-07 1.01E-07 0 27782 0 . . 0 0 . . 0 . . 0 .
v1 35.243 0 0 8.51E-08 0 0 0 0 0 0 0 0 0.01405 . . 0 0 . . 0 . . 0 .
v2 . . . . . . . . . . . . . . . . . . . . . . . .
v3 . . . . . . . . . . . . . . . . . . . . . . . .
v4 0 0 0 65.6625 0 0 0 0 0 0 0 0 0 . . 0.006083 0 . . 0 . . 0 .
v5 0 0 0 8.12E-08 19.1913 0 0 0 0 0 0 0 0 . . 0 0.00401 . . 0 . . 0 .
v6 . . . . . . . . . . . . . . . . . . . . . . . .
v7 . . . . . . . . . . . . . . . . . . . . . . . .
v8 0 0 0 0 0 0 0 74.6835 0 0 0 0 0 . . 0 0 . . 0.05917 . . 0 .
v9 . . . . . . . . . . . . . . . . . . . . . . . .
v10 . . . . . . . . . . . . . . . . . . . . . . . .
v11 0 0 0 2.09E-07 0 0 0 0 0 0 31.9055 0 0 . . 0 0 . . 0 . . 0.02825 .
v12 . . . . . . . . . . . . . . . . . . . . . . . .

 

Abdus-Salaam
Obsidian | Level 7

My supervisor has found a simple solution for my problem:

First he has created new start for the parameters of the model in proc NLmixed with proc mixed (means of the variants and the corresponding (standard errors)²) without the repeated /group=block value that i've had in the proc mixed code for creating start values for nlmixed. That he did since the block effect in my model was set 0, what i haven't recognised. Additionally he divided the grain dry matter values by 1000 wich doesn*t effect the CV since it's calculated with a quotient with two values that are equally effected by the transformation of the dry matter values.

Thank you again Mr. Denham and KSharp for your efforts to help me 

sas-innovate-white.png

Register Today!

Join us for SAS Innovate 2025, our biggest and most exciting global event of the year, in Orlando, FL, from May 6-9.

 

Save $200 when you sign up by March 14!

Register now!

What is ANOVA?

ANOVA, or Analysis Of Variance, is used to compare the averages or means of two or more populations to better understand how they differ. Watch this tutorial for more.

Find more tutorials on the SAS Users YouTube channel.

Discussion stats
  • 7 replies
  • 810 views
  • 2 likes
  • 3 in conversation