Matched Pair Tests and Comparing Data to Standards
Compare Data to a Standard Using a Matched Pair Test
Example 1: Use the cen_paired function
to determine if arsenic concentrations in groundwater exceed the
drinking water standard of 10 ug/L standard for the Example1
dataset.
data(Example1) # From NADA2 package
head(Example1,5L)
#> Arsenic NDis1 NDisTRUE
#> 1 4.00000 1 TRUE
#> 2 4.20000 0 FALSE
#> 3 0.61606 0 FALSE
#> 4 5.27628 0 FALSE
#> 5 3.00000 1 TRUE
with(Example1,cen_paired(Arsenic, NDisTRUE, 10, alt = "greater"))
#> Censored paired test for mean(Arsenic) equals 10
#> alternative hypothesis: true mean Arsenic exceeds 10.
#>
#> n = 21 Z= -20.4157 p-value = 1
#> Mean Arsenic = 2.252
The mean arsenic concentration does not exceed 10 ug/L.
Test for Differences in Two Paired Columns of Data
Example 2: Test whether atrazine concentrations were
the same in June versus September groundwaters on the same dates in a
variety of wells (rows – paired data). Test both for differences in the
mean as well as differences in the cdfs and the medians.
data(Atra); # From NADA package
head(Atra,5L)
#> June Sept JuneCen SeptCen
#> 1 0.38 2.66 FALSE FALSE
#> 2 0.04 0.63 FALSE FALSE
#> 3 0.01 0.59 TRUE FALSE
#> 4 0.03 0.05 FALSE FALSE
#> 5 0.03 0.84 FALSE FALSE
with(Atra,cen_paired(June, JuneCen, Sept, SeptCen))
#> Censored paired test for mean(June - Sept) equals 0.
#> alternative hypothesis: true mean difference does not equal 0.
#>
#> n = 24 Z= -1.0924 p-value = 0.2747
#> Mean difference = -3.927
The p-value is well above 0.05. Do not reject that the mean
difference in concentration for the two months could be 0.
# test for the median difference = 0 using the sign test.
with(Atra,cen_signtest(June, JuneCen, Sept, SeptCen))
#> Censored sign test for median(x:June - y:Sept) equals 0
#> alternative hypothesis: true median difference is not = 0
#> n = 24 n+ = 3 n- = 16 ties: 5
#>
#> No correction for ties: p-value = 0.004425
#> Fong correction for ties: p-value = 0.08956
Because it is important to correct for the numbers of tied values
within a pair, the p-value of 0.089 results in the conclusion to not
reject that the median difference in concentration between the two
months could be 0.
# test for a difference in the cdfs of the two months using the signed-rank
with(Atra,cen_signedranktest(June, JuneCen, Sept, SeptCen))
#> Censored signed-rank test for (x:June - y:Sept) equals 0
#> alternative hypothesis: true difference June - Sept does not equal 0
#>
#> Pratt correction for ties
#> n = 24, Z = -3.319, p.value = 0.0009033
The signed-rank test has more power to see differences than did the
sign test. It also is comparing the cdfs, the entire set of percentiles,
between the two months. It finds a difference because the upper end of
the distribution is quite a bit higher in the Sept data.
Comparing Data to Standards Using an Upper Confidence Limit
Using the Example 1 data, compute the UCL95 for censored data.
Step 1. Sample size. There are 21 observations.
Since it is on the borderline for deciding whether to use a
distributional or nonparametric method, both will be demonstrated
here.
Step 2. Distributional Method
Draw the boxplot for “censored data” (data with nondetects).
with(Example1,
cboxplot(Arsenic, NDisTRUE, Ylab="Arsenic Conc", show = TRUE))
Note that the highest detection limit is drawn as the horizontal
dashed line at 4 ug/L. Everything below that includes values estimated
using a lognormal ROS. Three “outliers” (not ‘bad data’) lie above the
estimated whisker, showing that the data are skewed.
Decide which of three distributions best fits the data using the
cenCompareCdfs command. Choose the distribution with the
smallest BIC.
with(Example1,
cenCompareCdfs (Arsenic, NDisTRUE, Yname = "Arsenic concentration in ug/L"))
The gamma distribution has the smallest BIC.
Note that the curve representing the normal distribution dips below
zero (x=0) at about the 10th percentile. A distribution of
concentrations with 10% negative numbers is not realistic, which results
in a higher BIC statistic.
Use the best-fit distribution (gamma) from 2b to compute the
UCL95.
egam <- with(Example1,
egammaAltCensored(Arsenic, NDisTRUE,
ci=TRUE, ci.type = "upper",
ci.method = "normal.approx"))
print(egam)
#>
#> Results of Distribution Parameter Estimation
#> Based on Type I Censored Data
#> --------------------------------------------
#>
#> Assumed Distribution: Gamma
#>
#> Censoring Side: left
#>
#> Censoring Level(s): 0.5 2.0 3.0 4.0
#>
#> Estimated Parameter(s): mean = 1.8399269
#> cv = 0.9131572
#>
#> Estimation Method: MLE
#>
#> Data: Arsenic
#>
#> Censoring Variable: NDisTRUE
#>
#> Sample Size: 21
#>
#> Percent Censored: 66.66667%
#>
#> Confidence Interval for: mean
#>
#> Confidence Interval Method: Normal Approximation
#>
#> Confidence Interval Type: upper
#>
#> Confidence Level: 95%
#>
#> Confidence Interval: LCL = -Inf
#> UCL = 2.575537
Use the print statement to get the “table format” for the output from
this EnvStats function. The UCL95 equals 2.57 assuming a gamma
distribution. Because this is lower than the 10 ug/L standard, the null
hypothesis of non-compliance is rejected, and the site from which these
data came is found to be in compliance.
Step 3. Nonparametric Method
Example 1:
There are multiple detection limits for this arsenic data. Compute
the Kaplan-Meier estimate of the mean and percentile bootstrap UCL95,
the latter because of the high percent of nondetects (66.67%) in the
data.
arsenic.out <- with(Example1,
enparCensored(Arsenic, NDisTRUE,
ci=TRUE, ci.method="bootstrap", ci.type="upper",
n.bootstraps=5000))
print(arsenic.out)
#>
#> Results of Distribution Parameter Estimation
#> Based on Type I Censored Data
#> --------------------------------------------
#>
#> Assumed Distribution: None
#>
#> Censoring Side: left
#>
#> Censoring Level(s): 0.5 2.0 3.0 4.0
#>
#> Estimated Parameter(s): mean = 1.7169702
#> sd = 1.5928374
#> se.mean = 0.1159666
#>
#> Estimation Method: Kaplan-Meier
#>
#> Data: Arsenic
#>
#> Censoring Variable: NDisTRUE
#>
#> Sample Size: 21
#>
#> Percent Censored: 66.66667%
#>
#> Confidence Interval for: mean
#>
#> Assumed Sample Size: 21
#>
#> Confidence Interval Method: Bootstrap
#>
#> Number of Bootstraps: 5000
#>
#> Number of Bootstrap Samples
#> With No Censored Values: 0
#>
#> Number of Times Bootstrap
#> Repeated Because Too Few
#> Uncensored Observations: 6
#>
#> Confidence Interval Type: upper
#>
#> Confidence Level: 95%
#>
#> Confidence Interval: Pct.LCL = 0.000000
#> Pct.UCL = 2.559268
#> BCa.LCL = 0.000000
#> BCa.UCL = 2.511469
#> t.LCL = 0.000000
#> t.UCL = 3.593110
The percentile bootstrap estimate of the UCL95 will be near to 2.53,
with the output varying slightly each time a bootstrap estimate is run.
This is essentially the same estimate as that for the gamma
distribution, with the identical result – the site is found to be in
compliance.
Example 2: Computation of a UCL95 for data with both
detected and non-detected values, DL unknown.
Data: Methyl Isobutyl Ketone (MIBK) in air above a medium-sized US
city.
A: Computation of the mean and UCL95
The MIBK concentrations are given as reported in column 1 – no
detection limit was provided. Nondetects were designated only as ND. The
lowest detected value in the data equals 0.1229. Assuming all ND values
are lower than this, all NDs were changed to <0.1229 as shown in the
MIBK and MIBKcen columns.
This results in only one reporting limit in the data, so the
Kaplan-Meier estimate will be biased a bit high. Instead, use the robust
ROS method with bootstrapping:
mibk.ucl95 <- with(Example2,
elnormAltCensored (MIBK, MIBKcen, method = "rROS",
ci=TRUE, ci.method = "bootstrap",
ci.type = "upper", n.bootstraps = 5000))
print(mibk.ucl95)
#>
#> Results of Distribution Parameter Estimation
#> Based on Type I Censored Data
#> --------------------------------------------
#>
#> Assumed Distribution: Lognormal
#>
#> Censoring Side: left
#>
#> Censoring Level(s): 0.1229
#>
#> Estimated Parameter(s): mean = 0.2160198
#> cv = 0.9338747
#>
#> Estimation Method: Imputation with
#> Q-Q Regression (rROS)
#>
#> Data: MIBK
#>
#> Censoring Variable: MIBKcen
#>
#> Sample Size: 31
#>
#> Percent Censored: 48.3871%
#>
#> Confidence Interval for: mean
#>
#> Confidence Interval Method: Bootstrap
#>
#> Number of Bootstraps: 5000
#>
#> Number of Bootstrap Samples
#> With No Censored Values: 0
#>
#> Number of Times Bootstrap
#> Repeated Because Too Few
#> Uncensored Observations: 0
#>
#> Confidence Interval Type: upper
#>
#> Confidence Level: 95%
#>
#> Confidence Interval: Pct.LCL = 0.0000000
#> Pct.UCL = 0.2883277
#> BCa.LCL = 0.0000000
#> BCa.UCL = 0.2666083
The percentile bootstrap UCL95 based on the robust ROS mean equals
0.290 (the Kaplan-Meier estimate with the slight bias would have equaled
0.293). Remember that your bootstrap result will slightly differ from
the one here. To decrease differences between runs, increase the number
of bootstraps, say to 10,000.
B: What if the detection limit had been known?
If a reporting limit of 0.029 had been provided by the laboratory,
the data would be as given in the MIBK2 and MIBK2cen columns. Using the
same procedure gives slightly lower results for both mean and UCL95:
mibk2.out <- with(Example2,
elnormAltCensored (MIBK2, MIBK2cen, method = "rROS",
ci=TRUE, ci.method = "bootstrap",
ci.type = "upper", n.bootstraps = 5000))
print(mibk2.out)
#>
#> Results of Distribution Parameter Estimation
#> Based on Type I Censored Data
#> --------------------------------------------
#>
#> Assumed Distribution: Lognormal
#>
#> Censoring Side: left
#>
#> Censoring Level(s): 0.029
#>
#> Estimated Parameter(s): mean = 0.2146941
#> cv = 0.9436391
#>
#> Estimation Method: Imputation with
#> Q-Q Regression (rROS)
#>
#> Data: MIBK2
#>
#> Censoring Variable: MIBK2cen
#>
#> Sample Size: 31
#>
#> Percent Censored: 48.3871%
#>
#> Confidence Interval for: mean
#>
#> Confidence Interval Method: Bootstrap
#>
#> Number of Bootstraps: 5000
#>
#> Number of Bootstrap Samples
#> With No Censored Values: 0
#>
#> Number of Times Bootstrap
#> Repeated Because Too Few
#> Uncensored Observations: 0
#>
#> Confidence Interval Type: upper
#>
#> Confidence Level: 95%
#>
#> Confidence Interval: Pct.LCL = 0.0000000
#> Pct.UCL = 0.2836123
#> BCa.LCL = 0.0000000
#> BCa.UCL = 0.2766389
The percentile bootstrap UCL95 based on the robust ROS mean will be
near to 0.290 with this known detection limit. It is always better to
use the laboratory documented limit, but not having one should not stop
the user from computing estimates using the lowest detected observation
as the limit.
Example 3: Computation of the expected percent of
observations exceeding a health advisory when all data are NDs. More
details of this method are found in Chapter 8 of Statistics for Censored
Environmental Data Using Minitab and R (Helsel, 2012).
All detection limits used are below the 10 ppb drinking water MCL for
arsenic. Therefore we know that 0 out of 14 current observations exceed
the MCL of 10 ppb. What is the range of percent of observations in the
aquifer that might exceed the MCL (with 95% probability)? Use the
binomial test command, entering the number of observations in the
dataset that exceed the MCL (0) and the number of total observations
(14). The alternative =”less” option states that this is a
one-sided confidence interval – we are looking only for possible
exceedances, nothing on the low end.
binom.test(0,14,alternative="less")
#>
#> Exact binomial test
#>
#> data: 0 and 14
#> number of successes = 0, number of trials = 14, p-value = 6.104e-05
#> alternative hypothesis: true probability of success is less than 0.5
#> 95 percent confidence interval:
#> 0.0000000 0.1926362
#> sample estimates:
#> probability of success
#> 0
Most of what is returned concerns a test for whether the proportion
of observations above the MCL differs from 50%, but this test is of no
interest here. What is of interest is the confidence interval on the
proportion of observations in the population that could be above the
MCL, based on the 14 samples observed. The UCL95 of the proportion
equals 0.192. Therefore we can say with 95% probability that there are
no more than 19.2% of concentrations in the aquifer exceeding the MCL –
we expect that there are fewer because the MCL of 10 is considerably
above the highest detection limit of 4 ppb, and this interval is
actually the probability of exceeding 4 ppb. Taking this conservative
approach that the probability of values falling above 4 ppb is the same
probability of falling above 10 ppb, the expected percent of samples at
this location above the MCL of 10 ppb is no more than 19.2%. This range
could be tightened by taking more samples, of course. For other
questions that can be answered when all values are nondetects, see
Chapter 8 in Helsel (2012).
Two-Group Tests
The MLE version of a “t-test” for censored data uses an MLE
regression with one X variable, a 0/1 group indicator. Assuming a normal
distribution:
with(CuZn,cen2means(Zn,ZnCen,Zone,LOG=FALSE))
#> MLE 't-test' of mean CensData: Zn by Factor: Zone
#> Assuming normal distribution of residuals around group means
#> mean of AlluvialFan = 11.49 mean of BasinTrough = 18.13
#> Chisq = 0.2928 on 1 degrees of freedom p = 0.588
#> Warning in cen2means(Zn, ZnCen, Zone, LOG = FALSE): NOTE: Data with nondetects may be projected below 0 with MLE normal distribution.
#> If so, p-values will be unreliable (often too small). Use perm test instead.
No difference between group means can be seen. But do the residuals
follow a normal distribution, as required? The corresponding Q-Q plot of
regression residuals and corresponding Shapiro-Francia test show that
the data do not follow a normal distribution.
MLE also builds a model of the two groups after estimating their mean
and standard deviations. When assuming a normal distribution with data
close to zero it is easy for the model to project data down below zero.
Besides being unrealistic, this can lead to a false separation between
the groups and p-values that are too low. To avoid this, either assume a
lognormal distribution (the default, or use LOG=TRUE) or use a
permutation test instead.
Use the cenperm2 function to perform a two-group permutation test.
This avoids an assumption of a normal distribution while testing for
differences in means:
with(CuZn,cenperm2(Zn,ZnCen,Zone))
#> Permutation test of mean CensData: Zn by Factor: Zone
#> 9999 Permutations alternative = two.sided
#> mean of AlluvialFan = 21.22 to 23.51 mean of BasinTrough = 21.28 to 21.94
#> Mean (AlluvialFan - BasinTrough) = -0.05612 to 1.567 p = 0.9999 to 0.9989
#>
The two groups do not have significantly different means. These are
p-values that you can believe, as they do not depend on the normal
assumption and do not project data values below zero.
Using cen2means to perform the test using the default lognormal
distribution, the difference in geometric means (not arithmetic means)
is performed below.
with(CuZn,cen2means(Zn,ZnCen,Zone))
#> MLE 't-test' of mean natural logs of CensData: Zn by Factor: Zone
#> Assuming lognormal distribution of residuals around group geometric means
#> geometric mean of AlluvialFan = 11.78 geometric mean of BasinTrough = 15.24
#> Chisq = 2.547 on 1 degrees of freedom p = 0.11
#>
No significant difference between geometric means. And the
corresponding Q-Q plot of regression residuals show that the data do not
follow a lognormal distribution either (one large outlier), but it is
the better fit of the two distributions. A permutation test using
cenperm2 after computing the logarithms would be a better test for
differences in geometric means.
“A nonparametric test between the two groups will test differences
between their cdfs (often considered as a test for difference in
medians) instead of testing the difference between means. The Peto-Peto
test is run using the cen1way function. It reports the Kaplan-Meier
medians in each of the groups:
with(CuZn,cen1way(Zn,ZnCen,Zone))
#> grp N PctND KMmean KMsd LCLmean
#> 1 AlluvialFan 67 23.88 22.70 74.03 4.796
#> 2 BasinTrough 50 8.00 21.61 18.77 16.490
#>
#> Oneway Peto-Peto test of CensData: Zn by Factor: Zone
#> Chisq = 5.183 on 1 degrees of freedom p = 0.0228
The two group medians (10 vs 17) are found different at p = 0.0228,
without assuming normality or substituting anything for the nondetects
censored at multiple (in this case, two) reporting limits. The cdfs for
the two groups are also shown, drawn using the cen_ecdf function.
with(CuZn,cen_ecdf (Zn,ZnCen,Zone,
Ylab="Zinc Concentration, in ug/L"))
Above approximately the 40th percentile or so, the Basin Trough
(dashed line) data are higher than (to the right of) the Alluvial Fan
data. Boxplots might show the group differences better than do the cdfs,
especially for non-statisticians.
with(CuZn,cboxplot(Zn,ZnCen,Zone,
Ylab="Zinc Concentration, in ug/L"))
Now for the two simpler tests, where all values below the highest DL
must be re-censored.
Contingency tables are a test to determine whether the proportions of
data in categories are the same in two or more groups. With nondetects
the cutoff level is again the highest reporting limit. Use the
TCE2 dataset with the highest of 4 DLs at 5 ug/L. The test
determines whether the proportions of data above versus below the cutoff
are the same in each group. It differs from a rank-sum in that all data
above the limit are simply in the same ‘Above’ group rather than ranked
individually. This loses some information as compared to the rank-sum
test. Density is a text variable (a factor) representing the groups,
while Below5Cens has a 1 for data below 5 ug/L and a 0 for data at or
above 5. First we combine them into a matrix using the ftable command,
then compute a ‘cross-tabulation’ with the xtabs(...)
command that the chisq.test(...) command expects. Finally
the chisquare test is computed. Its null hypothesis is that there is no
difference in the percent of data above the cutoff of 5 ug/L in the two
groups.
data(TCE2)
with(TCE2,ftable(Density~Below5Cens))
#> Density High Medium
#> Below5Cens
#> 0 18 12
#> 1 74 118
tab=with(TCE2,xtabs(~Below5Cens+Density))
chisq.test(tab)
#>
#> Pearson's Chi-squared test with Yates' continuity correction
#>
#> data: tab
#> X-squared = 4.0785, df = 1, p-value = 0.04343
The contingency table finds a difference in the proportions. It is
less powerful than the rank-sum test (which will also find a difference)
when there are data that can be ranked above the highest reporting
limit.
The nonparametric Wilcoxon rank-sum test can be calculated using the
wilcox.test(...) command. First you must have or create a
column that contains the concentrations for all detected values at the
maxDL and above, plus a single number (I use -1) below the max DL for
all values below the maxDL.
TCE2$Below5[TCE2$Below5Cens== 1] <- -1 # all <5s are now a -1
wilcox.test(Below5~Density,TCE2)
#>
#> Wilcoxon rank sum test with continuity correction
#>
#> data: Below5 by Density
#> W = 6599.5, p-value = 0.02713
#> alternative hypothesis: true location shift is not equal to 0
The two-sided p-value is significant at 0.027. The smaller p-value
than the contingency table test reflects the additional information in
the individual values at and above 5 that the rank-sum tests uses. The
slightly smaller p-value for the Peto-Peto test shows that it is the
most appropriate test when there are multiple detection limits.
Is re-censoring at the highest DL and running the rank-sum test
really better than the typical method of running a t-test on data with
one-half DL subbed for nondetects?
t.test(Half.DL~Density,TCE2)
#>
#> Welch Two Sample t-test
#>
#> data: Half.DL by Density
#> t = -0.065623, df = 201.88, p-value = 0.9477
#> alternative hypothesis: true difference in means between group High and group Medium is not equal to 0
#> 95 percent confidence interval:
#> -8.073324 7.553257
#> sample estimates:
#> mean in group High mean in group Medium
#> 7.763043 8.023077
No significant difference found. The arbitrariness of the
substitution process, not to mention that this puts the same number in
many times, and so likely decreases the standard deviation artificially,
should get you to quickly use a nonparametric test instead.
Three or more groups
The Golden data present lead concentrations in organs of
herons after exposing them to lead. There are four dosage groups
(“Dosage” or “Group” columns), zero plus three amounts of lead. The
objective was to determine if feathers or another non-destructive part
of the birds could be used as an environmental indicator, so it would
not be necessary to sacrifice a bird in order to measure their lead
concentrations. A censored boxplot of the Golden liver lead data shows
that the third and fourth groups have generally higher values.
data(Golden)
with(Golden,cboxplot(Liver, LiverCen, Dosage,
Ylab = "Lead concentrations in liver, in ppb"))
The skewness of the 0.05 group and the outlier of the 0.25 group
indicates that logs should end up being the better set of units to
use.
Kruskal-Wallis test
First the groups can be compared using a Kruskal-Wallis test, setting
all values below the highest detection limit of 0.04 as tied. Note that
there are detected observations below 0.04, so either the data had a
second and lower detection limit with no nondetects below it, or more
likely were reported using “insider censoring” (see Statistics for
Censored Environmental Data Using Minitab and R to find out what
that is and the problem it causes).
Step 1 - Create a variable – call it Below04 – with
zeros (or -1, or any value below the highest DL) for all data below the
highest DL of 0.04. Be careful not just to assign all 0.04s as
nondetects, as some of these could be detected 0.04s. Instead, use two
steps, the first to set all values BELOW 0.04 as a 0 (or -1), and the
second to set all data marked as nondetects (which will include the
<0.04 values) as a 0 (or -1). The result is a variable with an
indicator (-1 recommended) for all data below the highest reporting
limit, and original values for all detected data at and above the
highest reporting limit. The logical operators < (less-than) and ==
(equal to) are used here.
Golden$Below04 <- Golden$Liver
Golden$Below04[Golden$Liver<0.04] <- -1
Golden$Below04[Golden$LiverCen==TRUE] <- -1
Step 2 - run the Kruskal-Wallis test
kruskal.test(Below04~Dosage,Golden)
#>
#> Kruskal-Wallis rank sum test
#>
#> data: Below04 by Dosage
#> Kruskal-Wallis chi-squared = 7.8565, df = 3, p-value = 0.04907
The result shows that there is a difference (p = 0.049) between group
medians using this simple nonparametric test. An ANOVA on data after
substituting one-half DL will not find a difference (trust me on
this).
Peto-Peto test
The nonparametric Peto-Peto test, the multi-DL nonparametric test, is
computed using the cen1way(...) command:
with(Golden,cen1way(Liver,LiverCen,Dosage))
#> Warning: One or more group(s) do not have censored data.
#> Warning: One or more group(s) do not have censored data.
#> grp N PctND KMmean KMsd LCLmean
#> 1 0 7 28.57 0.1020 0.08834 0.03338
#> 2 0.01 7 28.57 0.1384 0.10590 0.07153
#> 3 0.05 6 0.00 12.1100 16.16000 -4.84900
#> 4 0.25 7 0.00 6.8660 17.46000 -9.28600
#>
#> Oneway Peto-Peto test of CensData: Liver by Factor: Dosage
#> Chisq = 7.795 on 3 degrees of freedom p = 0.0504
#>
#> Pairwise comparisons using Peto & Peto test
#>
#> data: CensData and Factor
#>
#> 0 0.01 0.05
#> 0.01 0.887 - -
#> 0.05 0.171 0.321 -
#> 0.25 0.079 0.127 0.887
#>
#> P value adjustment method: BH
The cdfs show that the higher two groups appear to differ in their
percentiles as compared to the lower two groups.
with(Golden,cen_ecdf(Liver,LiverCen,Dosage))
This is more easily seen by plotting the empirical cdfs in log
units:
Golden$lnLiver=log(Golden$Liver)
with(Golden,cen_ecdf(lnLiver,LiverCen,Dosage,
xlim = c(min(lnLiver), max(lnLiver)),
Ylab = "Natural Logs of Lead Concentrations in Liver"))
The 0.05 and 0.025 groups appear to have the higher liver lead
concentrations (are further to the right) than the other two groups.
The MLE “ANOVA”
For the parametric approach, use the cenanova(...)
command to run a censored regression with the groups as ‘factor’
explanatory variables. By default, cenanova assumes the
residuals follow a lognormal distribution, so use the associated Q-Q
plot to see if the residuals in log units appear approximately like a
normal distribution.
with(Golden,cenanova(Liver,LiverCen,Dosage))
#>
#> MLE test of mean natural logs of CensData: Liver by Factor: Dosage
#> Assuming lognormal distribution of CensData
#> Chisq = 10.67 on 3 degrees of freedom p = 0.0137
#>
#> mean(0) mean(0.01) mean(0.05) mean(0.25)
#> -2.989004 -2.701335 0.3031636 -0.7304445
#>
#> Simultaneous Tests for General Linear Hypotheses
#>
#> Multiple Comparisons of Means: Tukey Contrasts
#>
#>
#> Fit: survreg(formula = logCensData ~ Factor, dist = "gaussian")
#>
#> Linear Hypotheses:
#> Estimate Std. Error z value Pr(>|z|)
#> 0.01 - 0 == 0 0.2877 1.0841 0.265 0.9935
#> 0.05 - 0 == 0 3.2922 1.1036 2.983 0.0150 *
#> 0.25 - 0 == 0 2.2586 1.0625 2.126 0.1449
#> 0.05 - 0.01 == 0 3.0045 1.1010 2.729 0.0322 *
#> 0.25 - 0.01 == 0 1.9709 1.0599 1.860 0.2456
#> 0.25 - 0.05 == 0 -1.0336 1.0749 -0.962 0.7712
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> (Adjusted p values reported -- single-step method)
The overall test has a p-value of 0.014. Therefore the four group
mean logarithms (geometric means) differ. p-values for the individual
pairwise tests of differences show differences in two pairs of groups.
The p-values of 0.0150 (0.05 versus 0 groups) and 0.0323 (0.05 versus
0.01 groups) show that the 0.05 group differs from the lowest two
groups, but not from the 0.25 group. The residuals plot shows that log
are a very good set of units to use, as their residuals are close to a
normal distribution:
If instead you had wanted to test differences in the arithmetic
means, use a permutation test as a normal distribution will not fit
these data very well.
with(Golden,cenpermanova(Liver,LiverCen,Dosage))
#> Permutation test of mean CensData: Liver by Factor: Dosage
#> 9999 Permutations
#> Test Statistic = 1211 to 1211 p = 0.135 to 0.1365
#>
#> mean(0) mean(0.01) mean(0.05) mean(0.25)
#> 0.1000 0.1395 12.1100 6.8660
No significant difference in the means was found. This test did not
assume a normal distribution, though it is still influenced by outliers
because it evaluates means, which are influenced by outliers. The
permutation test will not extrapolate data to values below zero as would
MLE.
Note that the means of two groups, and so much of the data in the two
groups, falls below zero when estimated by MLE assuming a normal
distribution:
with(Golden,cenanova(Liver,LiverCen,Dosage,LOG=FALSE))
#>
#> MLE test of mean CensData: Liver by Factor: Dosage
#> Assuming normal distribution of CensData
#> Chisq = 6.889 on 3 degrees of freedom p = 0.0755
#>
#> NOTE: Data with nondetects may be projected below 0 with MLE normal distribution.
#> If so, p-values will be unreliable (often too small). Use perm test instead.
#>
#> mean(0) mean(0.01) mean(0.05) mean(0.25)
#> -2.889326 -2.847255 12.11417 6.865988
#>
#> Simultaneous Tests for General Linear Hypotheses
#>
#> Multiple Comparisons of Means: Tukey Contrasts
#>
#>
#> Fit: survreg(formula = CensData ~ Factor, dist = "gaussian")
#>
#> Linear Hypotheses:
#> Estimate Std. Error z value Pr(>|z|)
#> 0.01 - 0 == 0 0.04207 6.65879 0.006 1.000
#> 0.05 - 0 == 0 15.00349 6.71721 2.234 0.114
#> 0.25 - 0 == 0 9.75531 6.47075 1.508 0.433
#> 0.05 - 0.01 == 0 14.96142 6.71674 2.227 0.116
#> 0.25 - 0.01 == 0 9.71324 6.47026 1.501 0.437
#> 0.25 - 0.05 == 0 -5.24818 6.50033 -0.807 0.851
#> (Adjusted p values reported -- single-step method)
The p-value of 0.0755 is too small because the group differences are
exaggerated by pushing data down below 0. Given that the actual data
cannot go below zero, the cenpermanova p-value of 0.14 is a much more
realistic test result.
Two-factor “ANOVA” by MLE
The above commands cenpermanova, cenanova and cen1way are all “one
way” or one factor tests. They test for the effect of one and only one
factor (group rather than continuous) explanatory variable. A two factor
test is also available with the cen2way command in NADA2.
It can test for the effect of each of the two factors and the
interaction between the two factors, incorporating censored data by
maximum likelihood. The function may also be used on data without
nondetects, resulting in tests that are similar to but not necessarily
identical with the more common least-squares analysis of variance
software. A normal distribution of residuals is assumed by cen2way and
so the default is to log-transform the Y variable. If you decide that
other transformations are better than the log, transform (say using the
square root of Y) and then set the LOG= option to
FALSE. Here I use chromium concentrations in a stream and
test for the effect of a trend with the year (Yr) variable and a
seasonal effect with the (Season) variable. A better trend analysis
procedure will be performed in the Trend Analysis section later. Note
that Season is a text variable and Yr is numeric, so Yr will first be
set to a factor using the as.factor command.
# Load data
data(Gales_Creek)
Gales_Creek <- as.data.frame(Gales_Creek)
# set Year to factor by making a new variable.
Gales_Creek$Yr.f <- as.factor(Gales_Creek$Yr)
with(Gales_Creek,cen2way(TCr, CrND, Yr.f, Season))
#> Two-way Fixed Effects ANOVA by MLE Likelihood Ratio Tests for Censored Data
#> ln(TCr) modeled by Factor 1: Yr.f and Factor 2: Season
#> ANOVA Table with interaction
#> FACTORS chisquare df pval
#> Yr.f 24.03 9 4.254e-03
#> Season 36.03 1 1.943e-09
#> interaction 21.37 9 1.111e-02
#>
#> MODELS loglikelihood delta.lr0x2 AIC BIC
#> Yr.f only -69.28833 22.55 179.5767 221.7543
#> Season only -63.28820 34.55 151.5764 176.4830
#> both factors -61.95872 37.20 148.9174 173.8240
#> both + interaction -51.27145 58.58 145.5429 189.8795
#> Rescaled_Loglik_R2
#> 0.3261
#> 0.4576
#> 0.4834
#> 0.6562
The residuals plot and Shapiro-Francia test show that the log(Y) are
reasonable units. Both Yr and Season are significant with small
p-values. Their interaction is also significant with a p-value of 0.011.
The best model, judged by its lowest AIC of 145.54, is that all three
effects including the interaction, are significant. As mentioned in the
next section on regression, the BIC statistic tends to not find
variables significant that should be, based on simulation studies. Here
BIC would toss out the interaction tern and state that only the two main
effects are to be kept. BIC ‘best models’ tend to be too small and I
generally recommend using AIC instead. As with regression, models should
not be compared using AIC or R2 unless the Y scales are identical. Here
we have used only log(Y) so numerical measures of quality are
appropriate.
Correlation and Regression
First test for high vifs by computing a standard regression equation
that ignores the censoring indicator column, ignoring all results except
for the vifs:
vif(lm(AtraConc ~ Area + Applic + PctCorn + SoilGp + Temp + Precip + Dyplant + Pctl,Recon))
#> Area Applic PctCorn SoilGp Temp Precip Dyplant Pctl
#> 1.101992 2.739602 1.996707 1.480307 2.587299 2.206457 1.068839 1.131753
All of the variables appear uncorrelated with the others (all VIFs
well below 10). Therefore the p-values obtained in regression should be
reliable.
Step 1 - Create the 8-X variable
data.frame and run the regression with all variables.
Decide which scale the Y variable should be used (no transformation,
log, cube-root).
recon.8 <- with(Recon,data.frame(Area, Applic, PctCorn, SoilGp, Temp, Precip, Dyplant, Pctl))
reg.recon.8 <- with(Recon,cencorreg(AtraConc, AtraCen, recon.8))
#> Likelihood R2 = 0.6387 AIC = 804.4707
#> Rescaled Likelihood R2 = 0.6771 BIC = 843.968
#> McFaddens R2 = 0.3547157
#>
summary(reg.recon.8)
#>
#> Call:
#> survreg(formula = "log(AtraConc)", data = "Area+Applic+PctCorn+SoilGp+Temp+Precip+Dyplant+Pctl",
#> dist = "gaussian")
#> Value Std. Error z p
#> (Intercept) -8.76e+00 1.25e+00 -7.00 2.6e-12
#> Area 2.19e-05 1.98e-05 1.11 0.2685
#> Applic -2.75e-02 1.73e-02 -1.59 0.1127
#> PctCorn 5.91e-02 1.88e-02 3.14 0.0017
#> SoilGp 2.35e-01 3.66e-01 0.64 0.5198
#> Temp 6.25e-01 1.15e-01 5.44 5.3e-08
#> Precip -5.18e-03 1.39e-02 -0.37 0.7087
#> Dyplant -1.86e-02 1.51e-03 -12.36 < 2e-16
#> Pctl 4.17e-02 4.45e-03 9.36 < 2e-16
#> Log(scale) 5.96e-01 6.12e-02 9.75 < 2e-16
#>
#> Scale= 1.82
#>
#> Gaussian distribution
#> Loglik(model)= -391.7 Loglik(intercept only)= -607.1
#> Chisq= 430.68 on 8 degrees of freedom, p= 5.1e-88
#> Number of Newton-Raphson Iterations: 6
#> n= 423
The Rescaled likelihood R is fairly high (0.82) and the AIC equals
804.4. The Q-Q plot (below) shows a fairly straight pattern of data and
W = 0.966, so it would be difficult to find a better transformation of
the Y variable than the log. Use log Y.
Step 2 - Decide whether to transform the scale of
each X variable. Run the partplots procedure to see whether curvature in
the Y-X relationship means that a transformation of the X variable
should be taken.
layout(matrix(1:8,4,2))
with(Recon,partplots(AtraConc, AtraCen, recon.8,multiplot = F))
#> Area
#> untransformed
#> Likelihood R2 = 0.6387 AIC = 804.4707
#> cube root
#> Likelihood R2 = 0.6388 AIC = 804.4493
#> log transform
#> Likelihood R2 = 0.6389 AIC = 804.3342
#> Decrease in AIC from transformation of Area = 0.1364942
#>
#> Applic
#> untransformed
#> Likelihood R2 = 0.6387 AIC = 804.4707
#> cube root
#> Likelihood R2 = 0.6378 AIC = 805.5327
#> log transform
#> Likelihood R2 = 0.6371 AIC = 806.3659
#> Decrease in AIC from transformation of Applic = 0
#>
#> PctCorn
#> untransformed
#> Likelihood R2 = 0.6387 AIC = 804.4707
#> cube root
#> Likelihood R2 = 0.641 AIC = 801.7987
#> log transform
#> Likelihood R2 = 0.6423 AIC = 800.281
#> Decrease in AIC from transformation of PctCorn = 4.189691
#>
#> SoilGp
#> untransformed
#> Likelihood R2 = 0.6387 AIC = 804.4707
#> cube root
#> Likelihood R2 = 0.6385 AIC = 804.6955
#> log transform
#> Likelihood R2 = 0.6385 AIC = 804.7786
#> Decrease in AIC from transformation of SoilGp = 0
#>
#> Temp
#> untransformed
#> Likelihood R2 = 0.6387 AIC = 804.4707
#> cube root
#> Likelihood R2 = 0.6384 AIC = 804.8959
#> log transform
#> Likelihood R2 = 0.6381 AIC = 805.2458
#> Decrease in AIC from transformation of Temp = 0
#>
#> Precip
#> untransformed
#> Likelihood R2 = 0.6387 AIC = 804.4707
#> cube root
#> Likelihood R2 = 0.6387 AIC = 804.4702
#> log transform
#> Likelihood R2 = 0.6387 AIC = 804.4687
#> Decrease in AIC from transformation of Precip = 0.002027404
#>
#> Dyplant
#> untransformed
#> Likelihood R2 = 0.6387 AIC = 804.4707
#> cube root
#> Likelihood R2 = 0.6374 AIC = 806.063
#> log transform
#> Cannot take logs of zero or negative values.
#> Decrease in AIC from transformation of Dyplant = 0
#>
#> Pctl
#> untransformed
#> Likelihood R2 = 0.6387 AIC = 804.4707
#> cube root
#> Likelihood R2 = 0.6267 AIC = 818.3161
#> log transform
#> Cannot take logs of zero or negative values.
#> Decrease in AIC from transformation of Pctl = 0
#>
Only PctCorn showed an appreciable drop in AIC with a transformation.
As a percentage it is bounded by 0 and 100 so it’s a little odd to do a
transform. However the partial plot shows an increase in the percent of
“filled circle” detected observations and so a general increase in
atrazine with increasing PctCorn. The relationship may be nonlinear so
I’ll take the cube root of PctCorn.
Recon$cbrtPctCorn <- recon.8$PctCorn^(1/3)
recon.8onecube <- cbind(recon.8[, -3], Recon$cbrtPctCorn)
reg.recon.8onecube = with(Recon,cencorreg(AtraConc, AtraCen, recon.8onecube))
#> Likelihood R2 = 0.641 AIC = 801.7987
#> Rescaled Likelihood R2 = 0.6795 BIC = 841.296
#> McFaddens R2 = 0.3569165
#>
layout(matrix(1:8,4,2))
with(Recon,partplots(AtraConc, AtraCen, recon.8onecube,multiplot = F))
#> Area
#> untransformed
#> Likelihood R2 = 0.641 AIC = 801.7987
#> cube root
#> Likelihood R2 = 0.641 AIC = 801.8123
#> log transform
#> Likelihood R2 = 0.6411 AIC = 801.6904
#> Decrease in AIC from transformation of Area = 0.1083058
#>
#> Applic
#> untransformed
#> Likelihood R2 = 0.641 AIC = 801.7987
#> cube root
#> Likelihood R2 = 0.6404 AIC = 802.5295
#> log transform
#> Likelihood R2 = 0.6397 AIC = 803.316
#> Decrease in AIC from transformation of Applic = 0
#>
#> SoilGp
#> untransformed
#> Likelihood R2 = 0.641 AIC = 801.7987
#> cube root
#> Likelihood R2 = 0.6408 AIC = 802.0643
#> log transform
#> Likelihood R2 = 0.6407 AIC = 802.1689
#> Decrease in AIC from transformation of SoilGp = 0
#>
#> Temp
#> untransformed
#> Likelihood R2 = 0.641 AIC = 801.7987
#> cube root
#> Likelihood R2 = 0.6404 AIC = 802.5673
#> log transform
#> Likelihood R2 = 0.6399 AIC = 803.0917
#> Decrease in AIC from transformation of Temp = 0
#>
#> Precip
#> untransformed
#> Likelihood R2 = 0.641 AIC = 801.7987
#> cube root
#> Likelihood R2 = 0.641 AIC = 801.8317
#> log transform
#> Likelihood R2 = 0.641 AIC = 801.8476
#> Decrease in AIC from transformation of Precip = 0
#>
#> Dyplant
#> untransformed
#> Likelihood R2 = 0.641 AIC = 801.7987
#> cube root
#> Likelihood R2 = 0.6391 AIC = 804.098
#> log transform
#> Cannot take logs of zero or negative values.
#> Decrease in AIC from transformation of Dyplant = 0
#>
#> Pctl
#> untransformed
#> Likelihood R2 = 0.641 AIC = 801.7987
#> cube root
#> Likelihood R2 = 0.6287 AIC = 816.0618
#> log transform
#> Cannot take logs of zero or negative values.
#> Decrease in AIC from transformation of Pctl = 0
#>
#> Recon.cbrtPctCorn
#> untransformed
#> Likelihood R2 = 0.641 AIC = 801.7987
#> cube root
#> Likelihood R2 = 0.6419 AIC = 800.788
#> log transform
#> Likelihood R2 = 0.6423 AIC = 800.281
#> Decrease in AIC from transformation of Recon.cbrtPctCorn = 1.517664
#>
No other variables indicate a further transformation is necessary
after running partplots – the cbrtPctCorn variable is already
transformed so taking the log or cube root of the cube root would not
make much sense.
Step 3a - Can we lower the AIC by dropping
unimportant variables?
summary(reg.recon.8onecube)
#>
#> Call:
#> survreg(formula = "log(AtraConc)", data = "Area+Applic+SoilGp+Temp+Precip+Dyplant+Pctl+Recon.cbrtPctCorn",
#> dist = "gaussian")
#> Value Std. Error z p
#> (Intercept) -1.07e+01 1.58e+00 -6.74 1.6e-11
#> Area 2.04e-05 1.97e-05 1.03 0.30074
#> Applic -2.61e-02 1.62e-02 -1.61 0.10795
#> SoilGp 2.69e-01 3.62e-01 0.74 0.45644
#> Temp 6.41e-01 1.15e-01 5.57 2.6e-08
#> Precip -9.65e-03 1.41e-02 -0.68 0.49348
#> Dyplant -1.86e-02 1.50e-03 -12.44 < 2e-16
#> Pctl 4.16e-02 4.42e-03 9.40 < 2e-16
#> Recon.cbrtPctCorn 1.23e+00 3.50e-01 3.51 0.00046
#> Log(scale) 5.89e-01 6.12e-02 9.64 < 2e-16
#>
#> Scale= 1.8
#>
#> Gaussian distribution
#> Loglik(model)= -390.4 Loglik(intercept only)= -607.1
#> Chisq= 433.35 on 8 degrees of freedom, p= 1.4e-88
#> Number of Newton-Raphson Iterations: 6
#> n= 423
Both Precip and Soil Gp are very non-significant. Soil Group was a
number looked up in a book at the county scale and so wasn’t really
expected to provide much information (your insight as a scientist is
needed when using regression). Delete either Precip or Soil Group and
run the 7-variable model.
# an alternative way to make a data.frame without using data.frame(...)
recon.7 <- Recon[,c("Area", "Applic", "cbrtPctCorn", "Temp", "Precip", "Dyplant", "Pctl")]
reg.recon.7 <- with(Recon,cencorreg(AtraConc, AtraCen, recon.7))
#> Likelihood R2 = 0.6405 AIC = 800.3545
#> Rescaled Likelihood R2 = 0.679 BIC = 835.8021
#> McFaddens R2 = 0.3564587
#>
AIC has decreased so this is better than the 8-variable model. The
residuals plot looks much the same – this is expected in regression. The
normality of residuals is primarily determined by the scale of the Y
variable. Going to a six-variable model is an easy choice: Precip has a
high p-value.
recon.6 <- Recon[,c("Area", "Applic", "cbrtPctCorn", "Temp","Dyplant", "Pctl")]
reg.recon.6 <- with(Recon,cencorreg(AtraConc, AtraCen, recon.6))
#> Likelihood R2 = 0.6404 AIC = 798.5004
#> Rescaled Likelihood R2 = 0.6789 BIC = 829.8982
#> McFaddens R2 = 0.3563385
#>
summary(reg.recon.6)
#>
#> Call:
#> survreg(formula = "log(AtraConc)", data = "Area+Applic+cbrtPctCorn+Temp+Dyplant+Pctl",
#> dist = "gaussian")
#> Value Std. Error z p
#> (Intercept) -1.02e+01 1.35e+00 -7.55 4.3e-14
#> Area 2.27e-05 1.95e-05 1.16 0.24482
#> Applic -2.35e-02 1.59e-02 -1.48 0.13925
#> cbrtPctCorn 1.11e+00 3.19e-01 3.48 0.00051
#> Temp 6.06e-01 9.03e-02 6.71 2.0e-11
#> Dyplant -1.86e-02 1.49e-03 -12.45 < 2e-16
#> Pctl 4.17e-02 4.41e-03 9.45 < 2e-16
#> Log(scale) 5.90e-01 6.12e-02 9.65 < 2e-16
#>
#> Scale= 1.8
#>
#> Gaussian distribution
#> Loglik(model)= -390.8 Loglik(intercept only)= -607.1
#> Chisq= 432.65 on 6 degrees of freedom, p= 2.7e-90
#> Number of Newton-Raphson Iterations: 6
#> n= 423
AIC has decreased by 1.5. The next variable with a high p-value is
Area. The order of deleting these 3 variables likely wouldn’t matter and
you’d get to this 5 variable model even if you dropped them in a
different order.
recon.5 <- Recon[,c("Applic", "cbrtPctCorn", "Temp","Dyplant", "Pctl")]
reg.recon.5 <- with(Recon,cencorreg(AtraConc, AtraCen, recon.5))
#> Likelihood R2 = 0.6393 AIC = 797.8078
#> Rescaled Likelihood R2 = 0.6777 BIC = 825.1559
#> McFaddens R2 = 0.3552617
#>
summary(reg.recon.5)
#>
#> Call:
#> survreg(formula = "log(AtraConc)", data = "Applic+cbrtPctCorn+Temp+Dyplant+Pctl",
#> dist = "gaussian")
#> Value Std. Error z p
#> (Intercept) -9.85413 1.31831 -7.47 7.7e-14
#> Applic -0.02454 0.01594 -1.54 0.12368
#> cbrtPctCorn 1.10879 0.32018 3.46 0.00053
#> Temp 0.58724 0.08900 6.60 4.2e-11
#> Dyplant -0.01862 0.00150 -12.45 < 2e-16
#> Pctl 0.04155 0.00442 9.40 < 2e-16
#> Log(scale) 0.59339 0.06116 9.70 < 2e-16
#>
#> Scale= 1.81
#>
#> Gaussian distribution
#> Loglik(model)= -391.4 Loglik(intercept only)= -607.1
#> Chisq= 431.34 on 5 degrees of freedom, p= 5.2e-91
#> Number of Newton-Raphson Iterations: 6
#> n= 423
AIC has decreased by 0.7, so this is a very slightly better model.
The only other insignificant variable is Applic, the application amounts
of ag chemicals. These are known only on a county level so are estimated
by cutting and pasting county boundaries with watershed boundaries. They
also are voluntary amounts, and may not always provide accurate
information to the Federal government. But probably the main evidence
against the variable is its negative slope – we would expect more
atrazine to wash off with more applications. So we drop this variable to
see its effect.
recon.4 <- Recon[,c( "cbrtPctCorn", "Temp","Dyplant", "Pctl")]
reg.recon.4 <- with(Recon,cencorreg(AtraConc, AtraCen, recon.4))
#> Likelihood R2 = 0.6373 AIC = 798.192
#> Rescaled Likelihood R2 = 0.6756 BIC = 821.4904
#> McFaddens R2 = 0.3532981
#>
summary(reg.recon.4)
#>
#> Call:
#> survreg(formula = "log(AtraConc)", data = "cbrtPctCorn+Temp+Dyplant+Pctl",
#> dist = "gaussian")
#> Value Std. Error z p
#> (Intercept) -8.79657 1.10700 -7.95 1.9e-15
#> cbrtPctCorn 0.81989 0.25697 3.19 0.0014
#> Temp 0.51042 0.07332 6.96 3.4e-12
#> Dyplant -0.01869 0.00150 -12.47 < 2e-16
#> Pctl 0.04050 0.00437 9.26 < 2e-16
#> Log(scale) 0.59752 0.06121 9.76 < 2e-16
#>
#> Scale= 1.82
#>
#> Gaussian distribution
#> Loglik(model)= -392.6 Loglik(intercept only)= -607.1
#> Chisq= 428.96 on 4 degrees of freedom, p= 1.5e-91
#> Number of Newton-Raphson Iterations: 6
#> n= 423
The AIC goes up but only a little. This would be the scientist’s
choice to use either the 5-variable or the 4-variable model. I usually
choose the larger model if all p-values are under 0.10 because AIC and
similar metrics are known to choose too few variables. Here however the
Applic p-value in the 5-varable model is 0.124 and so I’d drop it,
choosing the 4-variable model. What is also behind my decision to drop
Applic is that it is a crude measure of amount of pesticide applied
(county level data cut and pasted) and cbrtPctCorn essentially measures
the same thing. Use your knowledge of the data to make your
decision.
Step 3b - Use the bestaic function to lower the
AIC?
Starting with the full 8 variables, though PctCorn has been
transformed to become cbrtPctCorn, run the bestaic function to see what
models the computer selects:
with(Recon,bestaic(AtraConc, AtraCen, recon.8onecube))
#> Evaluating 255 models and printing the 10 lowest AIC models
#> n.xvars model.xvars aic
#> 5 Applic Temp Dyplant Pctl Recon.cbrtPctCorn 797.8078
#> 4 Temp Dyplant Pctl Recon.cbrtPctCorn 798.1920
#> 6 Area Applic Temp Dyplant Pctl Recon.cbrtPctCorn 798.5004
#> 5 Area Temp Dyplant Pctl Recon.cbrtPctCorn 798.6990
#> 6 Applic SoilGp Temp Dyplant Pctl Recon.cbrtPctCorn 799.4286
#> 6 Applic Temp Precip Dyplant Pctl Recon.cbrtPctCorn 799.6515
#> 5 SoilGp Temp Dyplant Pctl Recon.cbrtPctCorn 800.0487
#> 5 Temp Precip Dyplant Pctl Recon.cbrtPctCorn 800.0798
#> 7 Area Applic SoilGp Temp Dyplant Pctl Recon.cbrtPctCorn 800.2693
#> 7 Area Applic Temp Precip Dyplant Pctl Recon.cbrtPctCorn 800.3545
Many models are very similar in AIC, but the ‘best’ is the 5-variable
model that we named recon.5, above. Second best was the
four variable model we called recon.4. So by deleting
sequentially we did get to the “best” models, but this is a lot quicker.
As you see in the list there are several other models around an AIC of
798, and if it were less expensive to use the variables in one of these,
it would be an excellent substitute for the mathematically lowest AIC
model. For example the fourth model down uses Area instead of Applic,
with all other variables the same. If Applic were expensive
to collect, this model has an AIC only 0.9 units higher.
Finding the best one-variable model. To find the
best 1-variable model (just to compare to the ATS equation), run the
four possible 1-variable models using the variables from the 4-variable
model.
reg.recon.cbrtPctCorn <- with(Recon,cencorreg(AtraConc, AtraCen, cbrtPctCorn))
#> Likelihood R2 = 0.01643 AIC = 1214.138
#> Rescaled Likelihood R2 = 0.01742 BIC = 1225.288
#> McFaddens R2 = 0.005773351
#>
reg.recon.Temp <- with(Recon,cencorreg(AtraConc, AtraCen, Temp))
#> Likelihood R2 = 0.08714 AIC = 1182.583
#> Rescaled Likelihood R2 = 0.09237 BIC = 1193.733
#> McFaddens R2 = 0.03176273
#>
reg.recon.Dyplant <- with(Recon,cencorreg(AtraConc, AtraCen, Dyplant ))
#> Likelihood R2 = 0.4759 AIC = 947.8357
#> Rescaled Likelihood R2 = 0.5045 BIC = 958.9849
#> McFaddens R2 = 0.2251064
#>
reg.recon.Pctl <- with(Recon,cencorreg(AtraConc, AtraCen, Pctl))
#> Likelihood R2 = 0.3455 AIC = 1041.861
#> Rescaled Likelihood R2 = 0.3662 BIC = 1053.011
#> McFaddens R2 = 0.1476647
#>
The Dyplant (days since planting) variable has the lowest AIC. Its
residuals plot (above) shows a linear pattern pulled away from the
normal theory line because of six high outliers. The negative slope is
reasonable: as there are more days since planting of corn, and atrazine
is only applied before planting, the more time it sits on the ground the
lower the amounts available to be washed off into the stream. The
outliers are probably high flows that would be accounted for by the Pctl
variable if it were in the model. Here are the details for the Dyplant
1-variable model:
summary(reg.recon.Dyplant)
#>
#> Call:
#> survreg(formula = "log(AtraConc)", data = "Dyplant", dist = "gaussian")
#> Value Std. Error z p
#> (Intercept) 1.34758 0.22401 6.02 1.8e-09
#> Dyplant -0.03063 0.00229 -13.40 < 2e-16
#> Log(scale) 0.96971 0.06270 15.47 < 2e-16
#>
#> Scale= 2.64
#>
#> Gaussian distribution
#> Loglik(model)= -470.4 Loglik(intercept only)= -607.1
#> Chisq= 273.31 on 1 degrees of freedom, p= 2.2e-61
#> Number of Newton-Raphson Iterations: 6
#> n= 423
Plots in Multiple Regression
To plot the predicted line from a regression model, set
cencorreg’s pred.plot option =
TRUE.
with(Recon,cencorreg(AtraConc, AtraCen, Dyplant,pred.plot=TRUE))
#> Likelihood R2 = 0.4759 AIC = 947.8357
#> Rescaled Likelihood R2 = 0.5045 BIC = 958.9849
#> McFaddens R2 = 0.2251064
#>
#> Call:
#> survreg(formula = "log(AtraConc)", data = "Dyplant", dist = "gaussian")
#>
#> Coefficients:
#> (Intercept) Dyplant
#> 1.3475770 -0.0306273
#>
#> Scale= 2.637173
#>
#> Loglik(model)= -470.4 Loglik(intercept only)= -607.1
#> Chisq= 273.31 on 1 degrees of freedom, p= <2e-16
#> n= 423
For the Dyplant only model the plot shows the relationship between Y
and Dyplant, but ignores the influence of all the other variables that
should be significant and in the model. It is too simple given that we
know there are 4 more variables in the ‘best’ model that affect
predictions which are not shown here. This isn’t much help when a
multiple regression model is the best model.
Plotting the partial plot of a multiple regression model provides a
better picture of the effect of one explanatory variable on Y after
‘subtracting off’ the effect of all the other variables in the model.
The partial plots of five pages earlier are an example. They are the
best way to picture the effect of one X variable within the context of
the multiple regression model. However they do not show what is actually
predicted for the Y variable from the model.
A prediction plot shows predicted values for Y for the full multiple
regression model using one X variable as the horizontal scale. The
predictions will not look like a nice straight line or smooth curve
unless there is only one variable in the model, as effects of the other
explanatory variables will cause the predictions to move up or down in
relation to the effect of the X variable chosen for the horizontal
scale. The cencorreg function uses the first variable in the X
data.frame as the one chosen for the horizontal axis on the plot. If
Dyplant as the ‘best’ X variable is to be used for the horizontal axis,
a new data frame with Dyplant as the first in the list must be
created:
recon.5.alt <- Recon[,c("Dyplant","Applic", "cbrtPctCorn", "Temp", "Pctl")]
plot.reg.recon.5 <- with(Recon,cencorreg(AtraConc, AtraCen, recon.5.alt,pred.plot = TRUE))
#> Likelihood R2 = 0.6393 AIC = 797.8078
#> Rescaled Likelihood R2 = 0.6777 BIC = 825.1559
#> McFaddens R2 = 0.3552617
#>
The pattern of predictions is much more jagged than the partial plot
of Dyplant. But the predicted values also show the effects
of nondetects and the other 4 variables in the model in addition to the
effect of the one X variable used in the plot.
The Nonparametric ATS line: Using Dyplant as the X
variable,
with(Recon,ATS(AtraConc, AtraCen, Dyplant))
#> Akritas-Theil-Sen line for censored data
#>
#> ln(AtraConc) = 3.3637 -0.0364 * Dyplant
#> Kendall's tau = -0.3995 p-value = 0
#>
Seeing this transformed back into the original units will look much
better.
with(Recon,ATS(AtraConc, AtraCen, Dyplant,retrans=TRUE))
#> Akritas-Theil-Sen line for censored data
#>
#> ln(AtraConc) = 3.3637 -0.0364 * Dyplant
#> Kendall's tau = -0.3995 p-value = 0
#>
For this dataset, the maximum likelihood and ATS slopes for Dyplant
are very similar (different by 0.006). The intercepts are similar as
well when you realize that a difference of 2.0 is small when
concentrations go up to 100.
The plot of the relationship of atrazine to flow percentile (Pctl)
shows a clear washoff effect at higher flows:
with(Recon,ATS(AtraConc,AtraCen, Pctl,retrans=TRUE))
#> Akritas-Theil-Sen line for censored data
#>
#> ln(AtraConc) = -5.4878 + 0.0889 * Pctl
#> Kendall's tau = 0.3465 p-value = 0
#>
There is not yet (June 2022) a good nonparametric “multiple
regression” method for censored data. There are ‘robust regression’
methods that perform nonparametric regression but I’ve never seen them
applied to censored data.
Trend Analysis for Censored Data
Load the Gales_Creek dataset:
Nonparametric Methods
ATS (no covariate or seasonal variation) We choose to use
the original units (LOG=FALSE) because the data appear
linear over time with one large outlier, and a nonparametric test will
not be overly influenced by one outlier. Running the ATS function on
concentration versus decimal time, we find strong evidence for a
downtrend (p = 0.006):
with(Gales_Creek,ATS(TCr,CrND,dectime,LOG=FALSE))
#> Akritas-Theil-Sen line for censored data
#>
#> TCr = 181.2617 -0.0896 * dectime
#> Kendall's tau = -0.234 p-value = 0.00648
#>
It isn’t easy to see on the plot, but the detection limits shown as
dashed lines are higher before 2012 as opposed to after 2012. The
methods of this section of the course work well with multiple detection
limits in the data record.
ATS on residuals from a smooth with a covariate Using the
centrend function, we first smooth the chromium – streamflow
relationship, and then test the residuals for trend:
with(Gales_Creek,centrend(TCr,CrND,discharge,dectime))
#> Trend analysis of TCr adjusted for discharge
#> Akritas-Theil-Sen line for censored data
#>
#> TCr residuals = 60.117 -0.0301 * dectime
#> Kendall's tau = -0.0579 p-value = 0.5051
#>
There is no trend in chromium concentration once the effect of
streamflow has been subtracted out. It appears that the evidence for a
downtrend was due to a change in the flow regime over the time period.
There is a strong relationship between flow and chromium
concentrations.
Seasonal Kendall test Perhaps there is a trend in either the
dry season alone, ignoring the effects of high flows on the trend test?
Perform the Seasonal Kendall test using the censeaken function and pay
attention to the individual season results by plotting them using the
seaplots = TRUE option.
with(Gales_Creek,censeaken(dectime,TCr,CrND,Season,seaplots=TRUE))
#>
#> DATA ANALYZED: TCr vs dectime by Season
#> ----------
#> Season N S tau pval Intercept slope
#> 1 Dry 34 -120 -0.214 0.069046 101.24 -0.05001
#> ----------
#> Season N S tau pval Intercept slope
#> 1 Wet 29 -83 -0.204 0.12381 233.15 -0.1151
#> ----------
#> Seasonal Kendall test and Theil-Sen line
#> reps_R N S_SK tau_SK pval intercept slope
#> 1 4999 63 -203 -0.21 0.012 181.26 -0.08965
#> ----------
There is an overall trend once the Seasonal Kendall test has removed
all comparisons between values in different seasons. Also, the dry
season has a p-value of 0.069. The prevailing wisdom in statistics in
2019 is to not get too rigid about an alpha of 0.05. A value of 0.069 is
close to 0.05 and the trend in the dry season graph appears strong. I
would report in this case that there is an overall downtrend and a
downtrend in the dry season. The high flows in the wet season were
preventing the non-seasonal centrend function from seeing
the trend.
Seasonal Kendall test with a covariate The centrendsea
function will perform the seasonal Kendall test on the residuals from a
GAM smooth of the original Y variable verses a covariate, like flow.
This first computes the same covariate adjustment of the centrend
function and using the residuals from the smooth, tests the
“covariate-adjusted trend” using the Seasonal Kendall test. The
variables to be entered are, in order, the Y variable to be tested for
trend, the T/F indicator of censoring for the Y variable, the covariate,
the time variable, and the season category variable.
with(Gales_Creek,centrendsea(TCr,CrND,discharge,dectime,Season))
#>
#> Trend analysis by Season of: TCr adjusted for discharge
#> ----------
#> Season N S tau pval Intercept slope
#> 1 Dry 34 -8 -0.0143 0.91638 34.52 -0.0172
#> ----------
#> Season N S tau pval Intercept slope
#> 1 Wet 29 -58 -0.143 0.28481 86.864 -0.04316
#> ----------
#> Seasonal Kendall test and Akritas-Theil-Sen line on residuals
#> reps_R N S_SK tau_SK pval intercept slope
#> 1 4999 63 -66 -0.0683 0.4292 59.654 -0.02967
#> ----------
There is no trend in flow-adjusted chromium found for this site and
dates. The slight decrease of the ATS line above, which includes the
influence of the below detection limit observations as well as detected
observations, is not significantly different from a zero slope.
Parametric Method
Simple Censored Regression Using the default log
transformation of chromium because we know there is one large outlier
lurking, the cencorreg function shows that the residuals
are not a normal distribution, but the data appear quite straight except
for the one high outlier. There is likely no better scale to work in –
untransformed concentrations would be far worse. Without deleting the
outlier (you should check it to see if there’s been an error, but you
can’t because this isn’t your data!), do not delete the outlier without
cause and work in the log units.
with(Gales_Creek,cencorreg(TCr,CrND,dectime))
#> Likelihood R2 = 0.07102 AIC = 163.4814
#> Rescaled Likelihood R2 = 0.07698 BIC = 168.9581
#> McFaddens R2 = 0.02880396
#>
#> Call:
#> survreg(formula = "log(TCr)", data = "dectime", dist = "gaussian")
#>
#> Coefficients:
#> (Intercept) dectime
#> 192.32546251 -0.09565766
#>
#> Scale= 0.9595278
#>
#> Loglik(model)= -78.2 Loglik(intercept only)= -80.6
#> Chisq= 4.64 on 1 degrees of freedom, p= 0.0312
#> n= 63
The regression p-value of 0.03 says that there is a trend. The slope
of – 0.095 log units per year will be approximately a 10% decrease in
chromium per year. But is this slope a good estimate, given that there
appear to be a confounding effect of streamflow? So perform a censored
multiple regression.
Censored Multiple Regression Create a data frame of the two
X variables, dectime and flow, and try again. This is a better model if
flow explains a lot of the variation in concentration. If that’s the
case the model’s AIC will be lower than the previous AIC of 163.48.
timeflow <- with(Gales_Creek,data.frame(dectime, discharge))
with(Gales_Creek,cencorreg(TCr,CrND,timeflow))
#> Likelihood R2 = 0.5926 AIC = 113.5493
#> Rescaled Likelihood R2 = 0.6424 BIC = 121.1848
#> McFaddens R2 = 0.351119
#>
#> Call:
#> survreg(formula = "log(TCr)", data = "dectime+discharge", dist = "gaussian")
#>
#> Coefficients:
#> (Intercept) dectime discharge
#> 220.206866401 -0.109654346 0.001290593
#>
#> Scale= 0.5499025
#>
#> Loglik(model)= -52.3 Loglik(intercept only)= -80.6
#> Chisq= 56.57 on 2 degrees of freedom, p= 5.19e-13
#> n= 63
The QQ plot looks great, and the residuals do not differ from a
normal distribution. The AIC is considerably lower for the 2-variable
model, so this model that accounts for flow variation should be used
instead of the original model.
Censored Multiple Regression with Seasonal Variables Sounds
like a menu option (‘seasonal vegetables’), doesn’t it? Create the sin
and cos function variables using 2*pi*dectime, and add it
to the stew. See if they add anything.
sinT <- with(Gales_Creek, sin(2*pi*dectime))
cosT <- with(Gales_Creek, cos(2*pi*dectime))
timeflowseas <- with(Gales_Creek,data.frame(dectime, discharge))
timeflowseas <- cbind(timeflowseas,sinT,cosT)
with(Gales_Creek, cencorreg(TCr, CrND, timeflowseas))
#> Likelihood R2 = 0.659 AIC = 106.3479
#> Rescaled Likelihood R2 = 0.7143 BIC = 118.3012
#> McFaddens R2 = 0.4206397
#>
#> Call:
#> survreg(formula = "log(TCr)", data = "dectime+discharge+sinT+cosT",
#> dist = "gaussian")
#>
#> Coefficients:
#> (Intercept) dectime discharge sinT cosT
#> 196.27125199 -0.09773029 0.00104809 0.21495631 0.30433249
#>
#> Scale= 0.5097264
#>
#> Loglik(model)= -46.7 Loglik(intercept only)= -80.6
#> Chisq= 67.77 on 4 degrees of freedom, p= 6.69e-14
#> n= 63
The QQ plot looks good. The sin and cos model has a lower AIC (106.3
versus the 2-variable model’s 113.5) so this is the best model of the
three. The slope of -0.098 per year still maps to around a 10% decrease
in concentration per year.
Logistic Regression
The ReconLogistic dataset presents atrazine concentrations at streams
across the midwestern United States. There were multiple detection
limits, and a health advisory of 1 ug/L. Several characteristics of the
basin at the time of sampling, including streamflow, are also
recorded.
data(ReconLogistic)
head(ReconLogistic,3)
#> CORNpct SOILGP PRECIP DYPLANT FPCTL APPLIC TEMP GT_1
#> 1 15 2.9 94.5 10 97 19.0 10.2 1
#> 2 15 2.9 94.5 130 19 19.0 10.2 0
#> 3 19 2.9 94.0 284 86 25.4 10.8 0
We will model the above/below 1 ug/L pattern using the
GT_1 variable. The primary assumption is that there is a
linear relationship between the X variables and the log(odds). Start by
checking VIFs for all 6 candidate variables. The glm
command using the family=binomial(logit) link function
produces the equation:
glm.1 <- glm(GT_1~
APPLIC+
CORNpct+
SOILGP+
PRECIP+
DYPLANT+
FPCTL,
ReconLogistic,family=binomial(logit))
vif(glm.1)
#> APPLIC CORNpct SOILGP PRECIP DYPLANT FPCTL
#> 1.802862 1.738165 1.467392 1.550877 1.172150 1.119903
There is no multicollinearity between the variables, as indicated by
VIF values of less than 2. Therefore the reported p-values should be
trustworthy.
summary(glm.1)
#>
#> Call:
#> glm(formula = GT_1 ~ APPLIC + CORNpct + SOILGP + PRECIP + DYPLANT +
#> FPCTL, family = binomial(logit), data = ReconLogistic)
#>
#> Deviance Residuals:
#> Min 1Q Median 3Q Max
#> -2.1267 -0.4117 -0.1715 0.3839 3.4336
#>
#> Coefficients:
#> Estimate Std. Error z value Pr(>|z|)
#> (Intercept) -6.323566 1.693157 -3.735 0.000188 ***
#> APPLIC 0.017530 0.020220 0.867 0.385984
#> CORNpct 0.034572 0.023263 1.486 0.137250
#> SOILGP 0.439216 0.503541 0.872 0.383070
#> PRECIP 0.039064 0.015949 2.449 0.014315 *
#> DYPLANT -0.016791 0.001919 -8.749 < 2e-16 ***
#> FPCTL 0.036820 0.006130 6.006 1.9e-09 ***
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#>
#> (Dispersion parameter for binomial family taken to be 1)
#>
#> Null deviance: 558.0 on 422 degrees of freedom
#> Residual deviance: 245.9 on 416 degrees of freedom
#> AIC: 259.9
#>
#> Number of Fisher Scoring iterations: 6
To compute the overall test of whether this model is better than no
model at all, the test that all slopes are zero, first compute the null
model by using a 1 instead of any X variables in the right-hand side of
the equation:
glm.null <- glm(GT_1~1,ReconLogistic,family=binomial(logit))
anova(glm.null,glm.1,test="Chisq")
#> Analysis of Deviance Table
#>
#> Model 1: GT_1 ~ 1
#> Model 2: GT_1 ~ APPLIC + CORNpct + SOILGP + PRECIP + DYPLANT + FPCTL
#> Resid. Df Resid. Dev Df Deviance Pr(>Chi)
#> 1 422 558.0
#> 2 416 245.9 6 312.1 < 2.2e-16 ***
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
The test statistic often named G equals 312.1. Compared to a
chi-squared distribution with 6 degrees of freedom, the difference in
the number of X variables between the two models, this statistic has a
p-value of 2 x 10^-16, and so is very significant. We conclude that
there is information in this model for predicting atrazine occurrence
above 1 and proceed to try and find the best model. The simplest way to
have the computer tell you the best logistic regression model is to use
the bestglm command. But let’s try manually first and see how we do. We
need to see if any of the X variables need to be transformed.
residualPlots(glm.1,type="deviance")
#> Test stat Pr(>|Test stat|)
#> APPLIC 2.2736 0.1315939
#> CORNpct 0.9137 0.3391386
#> SOILGP 8.8183 0.0029822 **
#> PRECIP 4.2759 0.0386572 *
#> DYPLANT 12.0510 0.0005177 ***
#> FPCTL 3.4084 0.0648658 .
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
It is hard to judge residuals plots with logistic regression because
there are separate groups of residuals for the Y = 0 and 1 data. The
smooth for DYPLANT appears curved, but is not an improvement. If
log(Dyplant) is substituted for Dyplant the
AIC increases, so is not a better model.
Using glm.1, Applic had the highest
p-value, so we’ll drop it and see if AIC goes below 259.9.
glm.3 <- glm(GT_1~
CORNpct+
SOILGP+
PRECIP+
DYPLANT+
FPCTL,
ReconLogistic,family=binomial(logit))
summary(glm.3)
#>
#> Call:
#> glm(formula = GT_1 ~ CORNpct + SOILGP + PRECIP + DYPLANT + FPCTL,
#> family = binomial(logit), data = ReconLogistic)
#>
#> Deviance Residuals:
#> Min 1Q Median 3Q Max
#> -2.2041 -0.4140 -0.1716 0.3764 3.4801
#>
#> Coefficients:
#> Estimate Std. Error z value Pr(>|z|)
#> (Intercept) -6.818507 1.599850 -4.262 2.03e-05 ***
#> CORNpct 0.045965 0.019277 2.384 0.01711 *
#> SOILGP 0.509550 0.495964 1.027 0.30424
#> PRECIP 0.043378 0.015074 2.878 0.00401 **
#> DYPLANT -0.016598 0.001876 -8.848 < 2e-16 ***
#> FPCTL 0.037835 0.005993 6.313 2.73e-10 ***
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#>
#> (Dispersion parameter for binomial family taken to be 1)
#>
#> Null deviance: 558.00 on 422 degrees of freedom
#> Residual deviance: 246.66 on 417 degrees of freedom
#> AIC: 258.66
#>
#> Number of Fisher Scoring iterations: 6
This 5-variable model is improved over glm.1.
SoilGp remains insignificant so let’s drop that and see the
effect on AIC.
glm.4 <- glm(GT_1~
CORNpct+
PRECIP+
DYPLANT+
FPCTL,
ReconLogistic,family=binomial(logit))
summary(glm.4)
#>
#> Call:
#> glm(formula = GT_1 ~ CORNpct + PRECIP + DYPLANT + FPCTL, family = binomial(logit),
#> data = ReconLogistic)
#>
#> Deviance Residuals:
#> Min 1Q Median 3Q Max
#> -2.1097 -0.4123 -0.1875 0.3781 3.3046
#>
#> Coefficients:
#> Estimate Std. Error z value Pr(>|z|)
#> (Intercept) -5.856860 1.267119 -4.622 3.80e-06 ***
#> CORNpct 0.039045 0.017943 2.176 0.029554 *
#> PRECIP 0.050493 0.013412 3.765 0.000167 ***
#> DYPLANT -0.016507 0.001866 -8.846 < 2e-16 ***
#> FPCTL 0.037770 0.005970 6.327 2.51e-10 ***
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#>
#> (Dispersion parameter for binomial family taken to be 1)
#>
#> Null deviance: 558.00 on 422 degrees of freedom
#> Residual deviance: 247.72 on 418 degrees of freedom
#> AIC: 257.72
#>
#> Number of Fisher Scoring iterations: 6
AIC is lower (improved). All terms are significant. We’ll settle on
this as our final model. This is the model found ‘best’ using
bestglm (note that bestglm requires the
data.frame to have the 0/1 Y variable in the last column,
which Recon does).
bestglm(ReconLogistic,family = binomial(logit), IC = "AIC")
#> Morgan-Tatar search since family is non-gaussian.
#> AIC
#> BICq equivalent for q in (0.388676345462891, 0.914863500091258)
#> Best Model:
#> Estimate Std. Error z value Pr(>|z|)
#> (Intercept) -6.53329819 1.181212560 -5.531010 3.183926e-08
#> CORNpct 0.04814197 0.018616040 2.586048 9.708340e-03
#> DYPLANT -0.01743824 0.002029210 -8.593611 8.427813e-18
#> FPCTL 0.03585721 0.006189977 5.792787 6.922797e-09
#> TEMP 0.50653496 0.102252080 4.953786 7.278310e-07
glm.4 has slopes with algebraic signs that make scientific sense, and
all explanatory variables are significant at alpha = 0.05. The model can
be compared to the original 6 variable model using either a partial
test, or with the AIC. The partial test determines whether the two
variables that were dropped add significantly to the explanatory power
of the model, just as in multiple linear regression. The null hypothesis
is that they do not; not rejecting the null hypothesis says to keep the
simpler model.
anova(glm.4, glm.1, test="Chisq")
#> Analysis of Deviance Table
#>
#> Model 1: GT_1 ~ CORNpct + PRECIP + DYPLANT + FPCTL
#> Model 2: GT_1 ~ APPLIC + CORNpct + SOILGP + PRECIP + DYPLANT + FPCTL
#> Resid. Df Resid. Dev Df Deviance Pr(>Chi)
#> 1 418 247.72
#> 2 416 245.90 2 1.8257 0.4014
We do not reject that the smaller model is just as good. Go with
glm.4. Finally, if we want to get an r-squared or a Brier
score, and get some plots of the relation between the log-odds and each
X variable, compute the same 4-variable model using the lrm
command from the rms package:
Recon.frame = with(ReconLogistic, datadist(CORNpct, DYPLANT, FPCTL, TEMP, GT_1))
options(datadist = "Recon.frame")
lrm4 <- lrm(GT_1 ~ CORNpct + DYPLANT + FPCTL + TEMP, data = ReconLogistic)
lrm4
#> Logistic Regression Model
#>
#> lrm(formula = GT_1 ~ CORNpct + DYPLANT + FPCTL + TEMP, data = ReconLogistic)
#>
#> Model Likelihood Discrimination Rank Discrim.
#> Ratio Test Indexes Indexes
#> Obs 423 LR chi2 323.35 R2 0.729 C 0.939
#> 0 266 d.f. 4 R2(4,423)0.530 Dxy 0.877
#> 1 157 Pr(> chi2) <0.0001 R2(4,296.2)0.660 gamma 0.877
#> max |deriv| 8e-07 Brier 0.076 tau-a 0.410
#>
#> Coef S.E. Wald Z Pr(>|Z|)
#> Intercept -6.5333 1.1812 -5.53 <0.0001
#> CORNpct 0.0481 0.0186 2.59 0.0097
#> DYPLANT -0.0174 0.0020 -8.59 <0.0001
#> FPCTL 0.0359 0.0062 5.79 <0.0001
#> TEMP 0.5065 0.1023 4.95 <0.0001
Multivariate Methods for Censored Data
Symonds et al (2016) used microbial source tracking (MST) markers to
detect fecal pollution in waters along the coast of Florida. Six MST
markers are in the dataset Markers in interval-censored
format, where (0 to MDL) indicate values below a limit of detection.
Nonzero lower ends of the interval indicate either (MDL to QL) data or
detected values above the QL. Also included is the US EPA total
enterococci marker Entero1A, a general fecal pollution
indicator.
Test whether the pattern of the six MST markers plus the
Entero1A indicator differs among the five sites using
ANOSIM.
Test whether there is a ‘trend’ (correlation) between the six MST
markers versus the general fecal pollution indicator using the Mantel
test for trend analysis using Kendall’s tau correlation between the
matrices (multivariate nonparametric correlation).
Solution
- Load the data and compute the ranks of the uscores; then compute the
anosim test for group differences and illustrate the results with an
MDS:
data(Markers)
head(Markers,3)
Mdat <- Markers[, -15] # remove the Site Name column
M.usc <- uscoresi(Mdat) # uscoresi drops rows with NAs (row 13 here)
M.euclid <- dist(M.usc)
Site <- Markers$Site_Name[-13] # delete the site entry for row 13 with NAs
M.anosim <- anosim(M.euclid, Site)
M.anosim
anosimPlot(M.anosim)
There is a difference between the five sites. To illustrate which
sites appear different than others, draw an NMDS plot:
uMDS(M.usc, group = Site,
legend.pos = "topright",
title = "NMDS of rank(uscores) for markers + entero")
From the left-right axis (NMDS1) we see that all three inlets (Port
of Miami, Port Everglades and Baker’s Haulover) are on the left side,
while Miami Central, an ocean outfall site, is on the right side. That
is the main contrast between sites. From the second axis (NMDS2) there
are one or two samples within Sites that are ‘outliers’ towards the top
as compared to others in that site. Some characteristic differs in those
samples. Sample 11 compared to the rest of Miami_North and site 14
compared to the rest of Baker’s Haulover, for example. Second axes for
NMDS aren’t always interpretable, but the first axis should be expected
to show the main difference detected by the ANOSIM result.
If you’d like to draw the NMDS plot manually, perhaps to change some
options from what is in the function, here is the R code that will draw
the same plot as the uMDS function:
M.nmds <- metaMDS(M.euclid)
Site <- as.factor(Site)
gp.color <- as.integer(Site)
Mplot <- ordiplot(M.nmds, type="none", display = "sites",
main="NMDS of rank(uscores) for markers + entero")
points(Mplot, "sites", pch=19, col=gp.color)
text(Mplot, "sites", pos=4, cex = 0.8)
leg.col <- c(1: length(levels(Site)))
legend("topright", legend=levels(Site), bty="n", col = leg.col, text.col = leg.col, pch = 19)
- Create two triangular distance matrices, one for the 6 MST markers
and the second for the general fecal pollution indicator (entero1A)
data. Then correlate the two matrices using the Mantel command. This is
to see if there is a ‘trend’ in the MST marker pattern with increasing
entero1A.
Mdata <- Markers[-13,] # delete row with NAs
M6 <- Mdata[, -(13:15)] # only keep columns for the 6 markers
M6.usc <- uscoresi(M6)
M6.euclid <- dist(M6.usc) # matrix for the 6 MST markers
ent <- Mdata[, 13:14] # the entero1A data
ent.usc<-uscoresi(ent)
#> [1] 29
ent.euclid<-dist(ent.usc) # matrix for the entero1A data
M6.Ktau <- mantel(ent.euclid, M6.euclid, method="kendall", permutations = 9999)
M6.Ktau
#>
#> Mantel statistic based on Kendall's rank correlation tau
#>
#> Call:
#> mantel(xdis = ent.euclid, ydis = M6.euclid, method = "kendall", permutations = 9999)
#>
#> Mantel statistic r: 0.3627
#> Significance: 1e-04
#>
#> Upper quantiles of permutations (null model):
#> 90% 95% 97.5% 99%
#> 0.0488 0.0673 0.0863 0.1136
#> Permutation: free
#> Number of permutations: 9999
There is a significant positive correlation between the MST marker
values and the entero1A values. One way to picture the correlation is to
plot their distance matrix entries against one another. The x,y pairs
are those in the triangular matrices that are being correlated using
Kendall’s tau in the mantel test. If there was a positive correlation,
for example, larger distances (relating small to large entero1A data)
would also have larger M6 distances between MST marker concentrations.
This is the pattern seen in the plot.
Site <- as.factor(Markers$Site_Name)
gp.color <- as.numeric(Site) # assigns numbers to group names in Site_Name
plot(ent.euclid, M6.euclid, pch = 19,
col = gp.color,
main = "Correlation of distance matrix of rank(uscores)")
lws <- lowess(ent.euclid, M6.euclid)
lines(lws)
legend("bottomright", legend=levels(Site), bty="n", col = 1:nlevels(Site), pch = 19)
The plot, in addition to being colorful, shows the general increase
of M6 distances as a function of entero1A distances, as shown by the
lowess smooth. As entero1A values increase, the pattern of 6 marker
concentrations also increases.
Extra Credit: Which MST marker(s) have the highest
correlation with the Entero1A values? This can be determined with the
bioenv command in the vegan package. This command performs iterative
mantel tests with subsets of the marker data. The highest correlation
coefficient is the set of best predictors, which may be 1 up to all 6 of
the markers.
bioenv(ent.euclid, M6.usc, method = "kendall")
#>
#> Call:
#> bioenv(comm = ent.euclid, env = M6.usc, method = "kendall")
#>
#> Subset of environmental variables with best correlation to community data.
#>
#> Correlations: kendall
#> Dissimilarities: euclidean
#> Metric: euclidean
#>
#> Best model has 1 parameters (max. 6 allowed):
#> usc.HF183_lo
#> with correlation 0.4591559
The highest correlation model is with one marker, HF183. This can be
illustrated by plotting the entero1A uscores against the HF183 uscores.
HF183 was the 4th of six columns within M6.
plot(ent.usc, M6.usc[,4], pch = 19)
You can plot ent.usc against the other five MST markers if you wish
to see which markers appear to be related to the Entero1A
concentrations. My look at the plots – BacHum is also correlated, and
HPyV has a binary style correlation – as Entero1A increases the
probability of a high rather than low PHyV increases – rather than a
linear relationship. This is because HPyV values occur in just two
categories, (0 to 249) and (250 to 499). Had all data been censored to
<500 this relationship would not have been visible. See the plot
below. HpyV was the 5th of the 6 marker columns. The other three markers
(all animal MST markers) show little correlation with Entero1A.
plot(ent.usc, M6.usc[,5], ylab = "HPyv rank of uscores")
