Extensions of Penalized Spline Regression for Small Area Estimation of Soil Profiles, Study Guides, Projects, Research of School management&administration

Download Extensions of Penalized Spline Regression for Small Area Estimation of Soil Profiles and more Study Guides, Projects, Research School management&administration in PDF only on Docsity! Extensions of Penalized Spline Regression for Natural Resource Monitoring Applications: Small Area Estimation of Soil Profiles F. Jay Breidt Colorado State University The work reported here was developed under STAR Research Assistance Agreements CR-829095 awarded by the U.S. Environmental Protection Agency (EPA) to Colorado State University. This presentation has not been formally reviewed by EPA. EPA does not endorse any products or commercial services mentioned in this report. 1 Environmental and Ecological Importance of Soils • Soils support life – viruses, bacteria, fungi, algae, protozoa, mites, nematodes, worms, insects and insect larvae, larger animals, and plant roots • Soils buffer environmental change – provide sink for greenhouse gases – regulate and partition water flow (infiltration vs. runoff) – decompose organic wastes (pesticides, sewage, etc.) – store and degrade nitrates, phosphorus, etc. 2 Soil Horizons and Soil Core Samples • Horizon thickness, number, and order vary from site to site • Use probe to select core sample; separate by horizon • Obtain lab analyses on cylindrical volume from each horizon 5 Horizon-Averaged Data: pH • pH vs. depth by map unit symbol 0 10 20 30 40 Depth (inches) 5 6 7 8 9 pH 9D2 0 10 20 30 40 Depth (inches) 5 6 7 8 9 pH 268F 0 10 20 30 40 Depth (inches) 4 5 6 7 8 9 pH 9B 0 10 20 30 40 Depth (inches) 4 5 6 7 8 9 pH 99D3 6 Sampled Cores from Soil Survey Update •Major Land Resource Area 107 Pilot Project: – two counties in western Iowa – multiphase sampling design – focus on third phase: horizon-specific lab data •Many study variables of interest – pH, carbon, sand, coarse silt, fine silt, clay, nitrogen, calcite, dolomite, CaCo3 •Want area-specific profiles – properties of soil as a function of depth – standard is a text description 7 Horizon-Averaged Data: Percent Clay • Clay vs. depth by map unit symbol 0 10 20 30 40 Depth (inches) 0 20 40 60 cl ay 99D 0 10 20 30 40 Depth (inches) 0 20 40 60 cl ay 9D 0 10 20 30 40 Depth (inches) 0 20 40 60 cl ay 99D2 0 10 20 30 40 Depth (inches) 0 20 40 60 cl ay 9E3 10 Key Data Features • Horizon averaging – irregular, wide; measurements not depth-specific – difficult to specify parametric model • Possible non-linear features (asymptotes) •Within-core dependence – horizons within same core may be correlated • Small sample sizes for many small areas – e.g., 97 map unit symbols, 193 soil cores •Many study variables 11 Small Area Estimation of Soil Profiles • Goal: borrow strength across small areas to obtain profiles for each map unit symbol 0 10 20 30 40 Depth (inches) 5 6 7 8 9 pH 1D3 0 10 20 30 40 Depth (inches) 5 6 7 8 9 pH 112D3 0 10 20 30 40 Depth (inches) 5 6 7 8 9 pH 10C3 0 10 20 30 40 Depth (inches) 5 6 7 8 9 pH 1E3 12 Extension: Continuous Domain • From finite population to infinite population θg = 1 |Dg| ∫ Dg θgs ds = 1 |Dg| ∫ Dg xTgsβ + zTgsu + 0  ds •Measure Ygs = θgs + gs • Need auxiliary information ∫ Dg xgs ds and ∫ Dg zgs ds 15 Extension: Continuous Parameter • Soil profile indexed by depth θg(t) = 1 |Dg| ∫ Dg θgs(t) ds = 1 |Dg| ∫ Dg xTgs(t)β + zTgs(t)u  ds • Need ∫ Dg xgs(t) ds and ∫ Dg zgs(t) ds •Measure horizon averages Ygsj = 1 dgsj − dgs,j−1 ∫ dgsj dgs,j−1 θgs(t) dt + gsj 16 Semiparametric Mixed Model: Longitudinal • Observation on subject i at jth time point • Subject-specific deviation is stochastic process Yij = x T ijβ +m(tij;wi) + z T ijui + Ui(tij) + ij e.g., Zhang, Lin, Raz, Sowers (1998) • Subject-specific deviation is smooth Yij = x T ijβ +m(tij;wi) + z T ijui +mi(tij) + ij e.g., Brumback and Rice (1998), Durbán, Harezlak, Wand, and Carroll 2005 • Need to adapt for small areas and horizon averages 17 Semiparametric Small Area Mixed Model: Horizons • jth horizon average in gth map unit symbol, sth core Ygsj = 1 dgsj − dgs,j−1 ∫ dgsj dgs,j−1 θgs(t) dt + gsj = β0 + β1 dgs,j−1 + dgsj 2 + K∑ k=1 ak 2 { (dgsj − κk)2+ − (dgs,j−1 − κk)2+ } +b0g + b1g dgs,j−1 + dgsj 2 + K∗∑ k=1 bk+1,g 2 { (dgsj − κ∗k)2+ − (dgs,j−1 − κ∗k)2+ } + 1 dgsj − dgs,j−1 ∫ dgsj dgs,j−1 Ugs(t) dt + gsj 20 Mixed Model Formulation of Penalized Splines • Regard semiparametric profile model as linear mixed model (Ruppert, Wand and Carroll 2003; Durbán et al. 2005) {ak} iid N(0, σ2a),   b0g b1g   iid N(0,Σ), {bk+1,g} iid N(0, σ2b), {gsj}iid N(0, σ2) • For now, assume simple within-core dependence model: {Ugs(t)} ≡ {Ugs} iid N(0, σ2U) 21 Model Fitting: Global Spline + Area Line • Straightforward estimation with standard software: S-Plus # # Set up global linear splines, then integrate over horizons. # Z1 <- outer(d1, knots, "-") Z1 <- Z1 * (Z1 > 0) Z2 <- outer(d2, knots, "-") Z2 <- Z2 * (Z2 > 0) Z <- 0.5 * diag(1/(d2 - d1)) %*% (Z2^2 - Z1^2) # # Fit global spline plus area-specific line using lme. Midd is horizon midpoint. # Adapted from Durban, Harezlak, Wand and Carroll 2005. # SplineCoeff <- factor(rep(1, length(y))) # fit <- lme(y ~ Midd, random = list(SplineCoeff = pdIdent( ~ Z - 1), MUS = pdSymm( ~ Midd), Core = pdIdent( ~ 1))) 22 Small Area Profile Estimates: Global Spline, Local Line • Nitrogen vs. depth profiles by map unit symbol 0 10 20 30 40 Depth (inches) 0. 05 0. 10 0. 15 0. 20 0. 25 ni tr og en 9D2 0 10 20 30 40 Depth (inches) 0. 05 0. 10 0. 15 0. 20 0. 25 ni tr og en 9B 0 10 20 30 40 Depth (inches) 0. 05 0. 10 0. 15 0. 20 0. 25 ni tr og en 99D3 0 10 20 30 40 Depth (inches) 0. 05 0. 10 0. 15 0. 20 0. 25 ni tr og en 10D2 25 Small Area Profile Estimates: Global Line, Local Line • Clay vs. depth profiles by map unit symbol 0 10 20 30 40 Depth (inches) 0 20 40 60 cl ay 527B 0 10 20 30 40 Depth (inches) 0 20 40 60 cl ay 212+ 0 10 20 30 40 Depth (inches) 0 20 40 60 cl ay 237 0 10 20 30 40 Depth (inches) 0 20 40 60 cl ay 237C 26 Types of Small Area Estimates for Soil Profiles • Point: classified on site as a given map unit symbol θ̂g(t)± 2σ̂ √√√√√Cgt(CTC + Λ̂)−1CTgt • Field: small unit with detailed soil map – field proportions by map unit symbol: ψ = (ψ1, . . . , ψG) T ψT (θ̂1(t), . . . , θ̂G(t)) T±2σ̂ √√√√ψTC∗t(CTC + Λ̂)−1CT∗tψ • Soil map unit: all delineations with same map unit symbol – replace ψ by ψ̂; adjust uncertainty 27