RBFY Pen-Level DMI Spatial Analysis

Description

This report characterizes spatial patterns of dry matter intake (DMI) across pens at Royal Beef Feedyard using operational data from June 2022 to current. The analytical unit is the pen-week. To be included, a pen-week was required to represent a designated feeding pen with a weekly average of at least 30 animals, have at least five days of feeding records during the week, and maintain consistent lot occupancy throughout the week. Hospital, receiving, processing, sort, and shipping pens were excluded.

Descriptive DMI heat map

June 2022 to current

MBW = metabolic body weight (BW^0.75)

summaries   = transpose(pen_summaries_ojs)
coords      = transpose(pen_coords_ojs)
metric_info = transpose(metric_meta_ojs)

// Dropdowns — metric first, defaults to first level (DMI, lb/hd/day)
viewof metric_sel = Inputs.select(metric_levels_ojs, {label: "Metric"})
viewof sex_sel    = Inputs.select(sex_levels_ojs,    {label: "Sex"})
viewof weight_sel = Inputs.select(weight_levels_ojs, {label: "Weight"})
viewof season_sel = Inputs.select(season_levels_ojs, {label: "Season"})

selected = summaries.filter(d =>
  d.metric === metric_sel &&
  d.sex    === sex_sel    &&
  d.weight === weight_sel &&
  d.season === season_sel
)

L = require("leaflet@1.9.4/dist/leaflet-src.js")
html`<link rel="stylesheet" href="https://unpkg.com/leaflet@1.9.4/dist/leaflet.css">`

{
  const geojson = JSON.parse(pen_geojson_ojs)

  function makeLookup(sel) {
    return Object.fromEntries(sel.map(d => [d.PEN, d]))
  }

  function getColor(val, minV, maxV) {
    if (val === null || val === undefined || isNaN(val)) return "#808080"
    const t = Math.max(0, Math.min(1, (val - minV) / (maxV - minV)))
    const stops = [
      [215,  48,  39],
      [252, 141,  89],
      [254, 224, 139],
      [255, 255, 191],
      [217, 239, 139],
      [145, 207,  96],
      [ 26, 152,  80]
    ]
    const scaled = t * (stops.length - 1)
    const lo     = Math.floor(scaled)
    const hi     = Math.min(lo + 1, stops.length - 1)
    const frac   = scaled - lo
    const r = Math.round(stops[lo][0] + frac * (stops[hi][0] - stops[lo][0]))
    const g = Math.round(stops[lo][1] + frac * (stops[hi][1] - stops[lo][1]))
    const b = Math.round(stops[lo][2] + frac * (stops[hi][2] - stops[lo][2]))
    return `rgb(${r},${g},${b})`
  }

function quantileSorted(values, p) {
    if (values.length === 0) return null
    const sorted = [...values].sort((a, b) => a - b)
    const idx    = (sorted.length - 1) * p
    const lo     = Math.floor(idx)
    const hi     = Math.ceil(idx)
    if (lo === hi) return sorted[lo]
    return sorted[lo] + (idx - lo) * (sorted[hi] - sorted[lo])
  }

  const div = document.createElement("div")
  div.style.height = "600px"
  div.style.width  = "100%"
  yield div

  const lats   = coords.map(d => d.lat)
  const lons   = coords.map(d => d.lon)
  const minLat = Math.min(...lats), maxLat = Math.max(...lats)
  const minLon = Math.min(...lons), maxLon = Math.max(...lons)

  const map = L.map(div, { minZoom: 15.25 })
    .fitBounds([[minLat - 0.01, minLon - 0.01],
                [maxLat + 0.01, maxLon + 0.01]])

  map.setMaxBounds([
    [minLat - 0.003, minLon - 0.003],
    [maxLat + 0.003, maxLon + 0.003]
  ])

  L.tileLayer(
    "https://server.arcgisonline.com/ArcGIS/rest/services/World_Imagery/MapServer/tile/{z}/{y}/{x}",
    { attribution: "Tiles © Esri" }
  ).addTo(map)

  let polyLayer = L.geoJSON(null).addTo(map)

  function updateMap(sel) {
    // Metadata for the selected metric (units label + decimals)
    const meta     = metric_info.find(m => m.label === metric_sel)
    const decimals = meta ? meta.decimals : 2
    // Units string for legend/popups: strip "DMI, " prefix from the label
    const units    = meta ? meta.label.replace(/^DMI,\s*/, "") : ""

    const lkp  = makeLookup(sel)
    const vals = Object.values(lkp)
      .map(d => d.mean_val)
      .filter(v => v !== null && v !== undefined && !isNaN(v))
    // Percentile-capped scale bounds (1st/99th) so outlier pens don't
    // compress the color range; values beyond the caps take the end colors.
    const minV = vals.length > 0 ? quantileSorted(vals, 0.005) : 0
    const maxV = vals.length > 0 ? quantileSorted(vals, 0.995) : 1

    polyLayer.clearLayers()
    polyLayer = L.geoJSON(geojson, {
      style: feature => {
        const d   = lkp[feature.properties.pen_id]
        const val = (d && d.mean_val !== null && !isNaN(d.mean_val))
                      ? d.mean_val : null
        return {
          fillColor:   getColor(val, minV, maxV),
          fillOpacity: 0.6,
          color:       "white",
          weight:      1
        }
      },
      onEachFeature: (feature, layer) => {
        const penId   = feature.properties.pen_id
        const d       = lkp[penId]
        const hasData = d && d.mean_val !== null && !isNaN(d.mean_val)
        const fmt     = v => (v !== null && v !== undefined && !isNaN(v))
                               ? `${Number(v).toFixed(decimals)} ${units}` : "N/A"
        layer.bindTooltip(
          hasData ? `Pen ${penId}: ${fmt(d.mean_val)}`
                  : `Pen ${penId}: no data`,
          { direction: "auto" }
        )
        layer.bindPopup(
          hasData
            ? `<b>Pen ${penId}</b><br>
               Mean: ${fmt(d.mean_val)}<br>
               Median: ${fmt(d.median_val)}<br>
               SD: ${fmt(d.sd_val)}<br>
               Pen-weeks: ${d.n_pen_weeks}`
            : `<b>Pen ${penId}</b><br>No data for this selection`
        )
      }
    }).addTo(map)

    const legend = L.control({ position: "bottomright" })
    legend.onAdd = function() {
      const el = L.DomUtil.create("div")
      el.style.cssText = `
        background: white; padding: 8px 10px; border-radius: 4px;
        box-shadow: 0 1px 4px rgba(0,0,0,0.3); font-size: 14px; line-height: 1.6;
      `
      const steps  = 100
      const height = 150
      let bars = ""
      for (let i = steps; i >= 0; i--) {
        const t   = i / steps
        const val = minV + t * (maxV - minV)
        bars += `<div style="width:16px; height:${height/steps}px;
                             background:${getColor(val, minV, maxV)};"></div>`
      }
      const nTicks = 5
      let ticks = ""
      for (let i = nTicks; i >= 0; i--) {
        const t   = i / nTicks
        const val = minV + t * (maxV - minV)
        const top = ((nTicks - i) / nTicks) * height
        const tickLabel = (i === nTicks)      ? `≤ ${val.toFixed(decimals)}`
                        : (i === 0)           ? `≥ ${val.toFixed(decimals)}`
                        :                        val.toFixed(decimals)
        ticks += `<div style="position:absolute; top:${top}px; left:22px;
                              font-size:14px; white-space:nowrap;
                              transform:translateY(-50%);">
                    ${tickLabel}
                  </div>`
      }
      el.innerHTML = `
        <b style="display:block; margin-bottom:4px;">${metric_sel}</b>
        <div style="position:relative; display:flex; align-items:stretch;">
          <div style="display:flex; flex-direction:column;">${bars}</div>
          <div style="position:relative; width:64px;">${ticks}</div>
        </div>
      `
      return el
    }
    if (map._legend) map.removeControl(map._legend)
    map._legend = legend
    legend.addTo(map)
  }

  invalidation.then(() => map.remove())
  updateMap(selected)
}

The color scale anchors to the minimum and maximum values within the currently selected combination (i.e., the legend scale changes with each change to filterable parameters).

Elevation and slope

Elevation heat map

Slope heat map

Elevation and slope values were derived from a 1-meter resolution bare-earth digital elevation model (DEM) produced by the USGS 3D Elevation Program, from statewide lidar collected for Kansas in 2018. For each pen, elevation statistics were summarized from all 1-meter DEM cells falling within the pen’s footprint, yielding a mean elevation, within-pen elevation range, and mean and maximum slope. Slope was calculated as the degree of elevation change between adjacent DEM cells, expressed as percent grade and averaged across each pen.

A few limitations are worth noting. The lidar predates the cattle intake data by several years, so any regrading or construction after 2018 is not reflected. The DEM’s vertical accuracy is on the order of 10 cm, which is small relative to the roughly 3-meter elevation spread across the yard but large enough that fine within-pen detail should be read cautiously. The DEM captures the bare-earth surface, so engineered drainage grade, manure mounds, and concrete bunk aprons all contribute to the elevation and slope values; the per-pen slope figures therefore reflect built drainage design more than natural terrain. Finally, elevation is a proxy for water-supply pressure and slope is a proxy for drainage and mud risk, but neither measures those conditions directly.

Baseline mixed model

A multivariable linear mixed model was fit to estimate each pen’s intake after adjusting for the major non-spatial sources of variation in DMI. The fixed effects below were included as adjustment terms so that differences in cattle type, ration, health, and season are not mistaken for pen-driven differences in intake:

• Sex, breed, background, risk classification, and origin
• Current estimated BW
• Days on feed
• Ration NEg
• Natural-program status
• Cumulative mortality rate
• Weekly morbidity rate
• Days since last implant
• Month of year

Lot and year were included as random effects. Lot accounts for intake being heavily shaped by the specific cohort of cattle occupying a pen, and year accounts for broad temporal shifts from year-to-year variability; including them keeps cohort and time-period effects from leaking into the pen estimates. Pen was also included as a random effect, and its estimated values are mapped - representing each pen’s tendency to feed above or below expectation once all of the above has been accounted for.

Model-adjusted heat map

The color scale shows each pen’s adjusted intake as a deviation from the yard average - these values are adjusted for cattle characteristics, management, and season, but not for physical pen attributes such as size, location, or bunk space. A pen that remains consistently above or below average therefore indicates a location effect worth further examination against physical pen features.

These adjusted values control for measured differences in cattle, ration, health, and season, but they cannot account for unmeasured cattle characteristics that influence intake. Where cattle are placed non-randomly - for example, if lighter cattle are concentrated in certain pens - any intake differences not captured by the measured variables (e.g., genetics of the lighter cattle) will be attributed to the pen rather than to the cattle. A consistently above- or below-average pen should therefore be read as a flag for further investigation rather than as a confirmed location effect.