roots (Программы с лекций (2013)), страница 2

TheMODEL statement includes Photoperiod, BAP, and their interactions in the model of the linear predictor.The DIST= option fits a zero-inflated Poisson model. The ZEROMODEL statement uses the default logitmodel to model the probability of a zero count and uses the variable Photoperiod as a linear predictor inthe model. The OUTPUT statement saves the predicted values and the estimated conditional zero-inflationprobabilities in the data set Zip. The predicted values and the zero-inflation probabilities are used later togenerate graphical displays that help assess the model’s goodness-of-fit.

The ODS OUTPUT statement savesthe goodness-of-fit statistics to the data set Fit so that a formal test for overdispersion can be performed. Ifthere is overdispersion, then the model is misspecified and the standard errors of the model parameters arebiased downwards.Output 1 displays the fit criteria for the ZIP model.Output 1 ZIP Model of Roots DataCriteria For Assessing Goodness Of FitCriterionDevianceScaled DeviancePearson Chi-SquareScaled Pearson X2Log LikelihoodFull Log LikelihoodAIC (smaller is better)AICC (smaller is better)BIC (smaller is better)DF260260Value1244.45661244.4566330.6476330.64761137.1695-622.22831264.45661265.30601300.4408Value/DF1.27171.2717Most of the criteria are useful only for comparing the model fit among given alternative models.

However, thePearson statistic can be used to determine if there is any evidence of overdispersion. If the model is correctlyspecified and there is no overdispersion, the Pearson chi-square statistic divided by the degrees-of-freedomhas an expected value of 1. The obvious question is whether the observed value of 1.2717 is significantlydifferent from 1, and thus an indication of overdispersion. As indicated in the section “Analysis” on page 1,the scaled Pearson statistic for generalized linear models has a limiting chi-square distribution under certainregularity conditions with degrees of freedom equal to the number of observations minus the number ofExample: Trajan Data Set F 9estimated parameters.

For Poisson and negative binomial models, the scale is fixed at 1, so there is nodifference between the scaled and unscaled versions of the statistic. Therefore, a formal one-sided test foroverdispersion is performed by computing the probability of observing a larger value of the statistic. Thefollowing SAS statements compute the p-value for such a test:data fit;set fit(where=(criterion="Scaled Pearson X2"));format pvalue pvalue6.4;pvalue=1-probchi(value,df);run;proc print data=fit noobs;var criterion value df pvalue;run;Output 2 reveals a p-value of 0.002 indicating rejection of the null hypothesis of no overdispersion at themost commonly used confidence levels.Output 2 Pearson Chi-Square StatisticCriterionScaled Pearson X2ValueDFpvalue330.64762600.0020Output 3 presents the parameter estimates for the ZIP model.

Because of the evidence of overdispersion,inferences based on these estimates are suspect; the standard errors are likely to be biased downwards.Nevertheless, the results as presented indicate that Photoperiod and BAP are significant determinants of theexpected value, as are three of the four interactions. Also as expected, Photoperiod is a significant predictorof the probability of a zero count.10 FOutput 3 ZIP Model Parameter EstimatesThe GENMOD ProcedureAnalysis Of Maximum Likelihood Parameter EstimatesParameterInterceptbapbapbapbapphotoperiodphotoperiodbap*photoperiodbap*photoperiodbap*photoperiodbap*photoperiodbap*photoperiodbap*photoperiodbap*photoperiodbap*photoperiodScale2.24.48.817.68162.22.24.44.48.88.817.617.6816816816816DFEstimateStandardError1111010101010000-2.15810.63220.52090.40580.00000.48570.0000-0.59740.0000-0.19980.0000-0.40740.00000.00000.00001.00000.10330.14490.15210.14680.00000.11930.00000.17390.00000.17600.00000.16860.00000.00000.00000.0000Wald 95% ConfidenceLimits-2.36070.34830.22280.11820.00000.25190.0000-0.93830.0000-0.54480.0000-0.73780.00000.00000.00001.0000WaldChi-SquarePr > ChiSq436.1419.0411.737.65.16.58.11.80.1.29.5.84...<.0001<.00010.00060.0057.<.0001.0.0006.0.2562.0.0157...-1.95560.91620.81910.69350.00000.71950.0000-0.25650.00000.14510.0000-0.07690.00000.00000.00001.0000NOTE: The scale parameter was held fixed.Analysis Of Maximum Likelihood Zero Inflation Parameter EstimatesParameterInterceptphotoperiodphotoperiod816DFEstimateStandardError110-0.1033-4.16980.00000.17660.76110.0000Wald 95% ConfidenceLimits-0.4495-5.66170.00000.2429-2.67800.0000WaldChi-SquarePr > ChiSq0.3430.01.0.5587<.0001.Another method for assessing the goodness-of-fit of the model is to compare the observed relative frequenciesof the various counts to the maximum likelihood estimates of their respective probabilities.

The followingSAS statements demonstrate one method of computing the estimated probabilities and generating twocomparative plots.The first step is to observe the value of the largest count and the sample size and save them into macrovariables.proc means data=Trajan noprint;var roots;output out=maxcount max=max N=N;run;data _null_;set maxcount;call symput('N',N);call symput('max',max);run;Example: Trajan Data Set F 11%let max=%sysfunc(strip(&max));Next, you use the model predictions and the estimated zero-inflation probabilities that are stored in the outputdata set Zip to compute the conditional probabilities Pr.yij D i jxij /. These are the variables ep0–ep&maxin the following DATA step. You also generate an indicator variable for each count i , i D 0; 1; : : : ; &max,where each observation is assigned a value of 1 if count i is observed, and 0 otherwise.

These are the variablesc0–c&max.data zip(drop= i);set zip;lambda=pred/(1-pzero);array ep{0:&max} ep0-ep&max;array c{0:&max} c0-c&max;do i = 0 to &max;if i=0 then ep{i}= pzero + (1-pzero)*pdf('POISSON',i,lambda);elseep{i}=(1-pzero)*pdf('POISSON',i,lambda);c{i}=ifn(roots=i,1,0);end;run;Now you can use PROC MEANS to compute the means of the variables ep0, : : : , ep&max and c0, : : : ,c&max. The means of ep0, : : : , ep&max are the maximum likelihood estimates of Pr.y D i /. The means ofc0, : : : , c&max are the observed relative frequencies.proc means data=zip noprint;var ep0 - ep&max c0-c&max;output out=ep(drop=_TYPE_ _FREQ_) mean(ep0-ep&max)=ep0-ep&max;output out=p(drop=_TYPE_ _FREQ_) mean(c0-c&max)=p0-p&max;run;The output data sets from PROC MEANS are in what is commonly referred to as wide form. That is, thereis one observation for each variable.

In order to generate comparative plots, the data need to be in whatis referred to as long form. Ultimately, you need four variables, one whose observations are an index ofthe values of the counts, a second whose observations are the observed relative frequencies, a third whoseobservations contain the ZIP model estimates of the probabilities Pr.y D i /, and a fourth whose observationscontain the difference between the observed relative frequencies and the estimated probabilities.The following SAS statements transpose the two output data sets so that they are in long form. Then, thetwo data sets are merged and the variables that index the count values and record the difference between theobserved relative frequencies and the estimated probabilities are generated.proc transpose data=ep out=ep(rename=(col1=zip) drop=_NAME_);run;proc transpose data=p out=p(rename=(col1=p) drop=_NAME_);run;12 Fdata zipprob;merge ep p;zipdiff=p-zip;roots=_N_ -1;label zip='ZIP Probabilities'p='Relative Frequencies'zipdiff='Observed minus Predicted';run;Now you can use the SGPLOT procedure to produce the comparative plots.proc sgplot data=zipprob;scatter x=roots y=p /markerattrs=(symbol=CircleFilled size=5px color=blue);scatter x=roots y=zip /markerattrs=(symbol=TriangleFilled size=5px color=red);xaxis type=discrete;run;proc sgplot data=zipprob;series x=roots y=zipdiff /lineattrs=(pattern=ShortDash color=blue)markers markerattrs=(symbol=CircleFilled size=5px color=blue);refline 0/ axis=y;xaxis type=discrete;run;Figure 5 Comparison of ZIP Probabilities to Observed Relative FrequenciesZIP Probabilities versus Relative FrequenciesObserved Relative Frequencies Minus ZIP ProbabilitiesFigure 5 shows that the ZIP model accounts for the excess zeros quite well and that the ZIP distributionreasonably captures the shape of the distribution of the relative frequencies.Clearly, a zero-inflated model can account for the excess zeros.

However, because the Pearson statisticindicates that there is evidence of model misspecification, with overdispersion being the most likely culprit,inference based upon the ZIP model estimates are suspect. If overdispersion is the culprit, then fitting azero-inflated negative binomial (ZINB) might be a solution because it can account for the excess zeros aswell as the ZIP model did and it provides a more flexible estimator for the variance of the response variable.The following SAS statements fit a ZINB model to the response variable Roots.

The model specification isthe same as before except that the DIST= option in the MODEL statement now specifies a ZINB distribution.Example: Trajan Data Set F 13proc genmod data=Trajan;class bap photoperiod;model roots = bap|photoperiod / dist=zinb offset=lshoot;zeromodel photoperiod;output out=zinb predicted=pred pzero=pzero;ods output ParameterEstimates=zinbparms;ods output Modelfit=fit;run;Output 4 displays the fit criteria for the ZINB model. The Pearson chi-square statistic divided by itsdegrees-of-freedom is 1.0313, which is much closer to 1 compared to the ZIP model.Output 4 ZINB Model of Roots DataCriteria For Assessing Goodness Of FitCriterionDFDevianceScaled DeviancePearson Chi-SquareScaled Pearson X2Log LikelihoodFull Log LikelihoodAIC (smaller is better)AICC (smaller is better)BIC (smaller is better)ValueValue/DF1232.45091232.4509268.1486268.1486-616.2255-616.22551254.45091255.47421294.03362602601.03131.0313The following SAS statements perform the same formal test that was used for the ZIP model:data fit;set fit(where=(criterion="Scaled Pearson X2"));format pvalue pvalue6.4;pvalue=1-probchi(value,df);run;proc print data=fit noobs;var criterion value df pvalue;run;Output 5 reveals a p-value of 0.3509, which indicates that you would fail to reject the null hypothesis of nooverdispersion at the most commonly used confidence levels.Output 5 Pearson Chi-Square StatisticCriterionScaled Pearson X2ValueDFpvalue268.14862600.350914 FTable 2 provides a side-by-side comparison of the other fit criteria for the two models.

All of the criteriafavor the ZINB over the ZIP model.Table 2 Comparison of ZIP and ZINB Model Fit CriteriaCriterionZIPZINBFull Log LikelihoodAICAICCBIC–622.22831264.45661265.30601300.4408–616.22551254.45091255.47421294.0336Output 6 displays the ZINB model’s parameter estimates. Compared to the ZIP model, most (but not all)of the ZINB model parameters are slightly smaller in magnitude and the standard errors are larger. Thereis effectively no change in any inference you would make regarding any of the parameters. The negativebinomial dispersion parameter has an estimated value of 0.0649, and the Wald 95% confidence intervalindicates that the estimate is significantly different from 0.Output 6 ZINB Model Parameter EstimatesThe GENMOD ProcedureAnalysis Of Maximum Likelihood Parameter EstimatesParameterInterceptbapbapbapbapphotoperiodphotoperiodbap*photoperiodbap*photoperiodbap*photoperiodbap*photoperiodbap*photoperiodbap*photoperiodbap*photoperiodbap*photoperiodDispersion2.24.48.817.68162.22.24.44.48.88.817.617.6816816816816DFEstimateStandardError1111010101010001-2.16630.63710.52350.40950.00000.48750.0000-0.59600.0000-0.19620.0000-0.40480.00000.00000.00000.06490.11880.17020.17770.16970.00000.13970.00000.20550.00000.20840.00000.19790.00000.00000.00000.0240Wald 95% ConfidenceLimits-2.39920.30360.17530.07690.00000.21370.0000-0.99880.0000-0.60470.0000-0.79270.00000.00000.00000.0314WaldChi-SquarePr > ChiSq332.2214.028.685.82.12.17.8.41.0.89.4.18...<.00010.00020.00320.0158.0.0005.0.0037.0.3466.0.0408...-1.93330.97060.87170.74210.00000.76140.0000-0.19310.00000.21230.0000-0.01690.00000.00000.00000.1340NOTE: The negative binomial dispersion parameter was estimated by maximum likelihood.Analysis Of Maximum Likelihood Zero Inflation Parameter EstimatesParameterInterceptphotoperiodphotoperiod816DFEstimateStandardError110-0.1150-4.29240.00000.17790.86940.0000Wald 95% ConfidenceLimits-0.4636-5.99630.00000.2336-2.58850.0000WaldChi-SquarePr > ChiSq0.4224.38.0.5179<.0001.You can generate comparative plots the same way as before with a few minor alterations.