- Mark as New
- Bookmark
- Subscribe
- Mute
- RSS Feed
- Permalink
- Report Inappropriate Content
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 |
0 |
2 |
1 |
1 |
3 |
2 |
0 |
4 |
2 |
1 |
5 |
3 |
0 |
6 |
3 |
1 |
7 |
4 |
0 |
8 |
4 |
1 |
9 |
5 |
0 |
10 |
6 |
0 |
11 |
7 |
0 |
12 |
7 |
1 |
13 |
8 |
0 |
14 |
8 |
1 |
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; ;
Accepted Solutions
- Mark as New
- Bookmark
- Subscribe
- Mute
- RSS Feed
- Permalink
- Report Inappropriate Content
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
- Mark as New
- Bookmark
- Subscribe
- Mute
- RSS Feed
- Permalink
- Report Inappropriate Content
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
- Mark as New
- Bookmark
- Subscribe
- Mute
- RSS Feed
- Permalink
- Report Inappropriate Content
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; ;
- Mark as New
- Bookmark
- Subscribe
- Mute
- RSS Feed
- Permalink
- Report Inappropriate Content
- Mark as New
- Bookmark
- Subscribe
- Mute
- RSS Feed
- Permalink
- Report Inappropriate Content
Hello KSharp,
thank you for the document. It helped me to understand the warning better.
- Mark as New
- Bookmark
- Subscribe
- Mute
- RSS Feed
- Permalink
- Report Inappropriate Content
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.
- Mark as New
- Bookmark
- Subscribe
- Mute
- RSS Feed
- Permalink
- Report Inappropriate Content
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 | . | . | . | . | . | . | . | . | . | . | . | . | . | . | . | . | . | . | . | . | . | . | . | . |
- Mark as New
- Bookmark
- Subscribe
- Mute
- RSS Feed
- Permalink
- Report Inappropriate Content
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