diff --git a/.github/workflows/build_and_upload_wheels.yaml b/.github/workflows/build_and_upload_wheels.yaml
index 2f2a7dfe..9272f300 100644
--- a/.github/workflows/build_and_upload_wheels.yaml
+++ b/.github/workflows/build_and_upload_wheels.yaml
@@ -21,7 +21,7 @@ jobs:
     runs-on: ubuntu-22.04
     strategy:
       matrix:
-        python-version: [3.7, 3.8, 3.9, '3.10', 3.11]
+        python-version: ['3.11', '3.12', '3.13']
     steps:
       - uses: actions/checkout@v3
 
diff --git a/.github/workflows/coverage.yaml b/.github/workflows/coverage.yaml
index 55b17d92..096c13ef 100644
--- a/.github/workflows/coverage.yaml
+++ b/.github/workflows/coverage.yaml
@@ -14,7 +14,7 @@ jobs:
     strategy:
       matrix:
         os: [ubuntu-22.04]
-        python-version: [3.9, "3.10", "3.11", "3.12", "3.13"]
+        python-version: ["3.11", "3.12", "3.13"]
 
     steps:
     - uses: actions/checkout@v4
@@ -24,20 +24,12 @@ jobs:
       with:
         python-version: ${{ matrix.python-version }}
 
-    - name: Install Python Dependencies for 3.12 and later
-      if: ${{ matrix.python-version >= '3.12' }}
+    - name: Install Python3 Dependencies
       run: |
         python -m pip install --upgrade pip pytest setuptools wheel
-        pip install -r requirements.txt
-        # Optional Dependency for HDF Checkpointing
-        pip install tables
-        python setup.py install
-        python setup.py build_ext --inplace
-        python -m pip list
-
-    - name: Standard Install Python3 Dependencies
-      run: |
-        python -m pip install --upgrade pip pytest
+        # Install hatchet first: pandas 3 support is on its develop branch, and
+        # requirements.txt cannot be resolved against released hatchet versions
+        python -m pip install git+https://github.com/LLNL/hatchet.git@develop
         pip install -r requirements.txt
         # Optional Dependency for HDF Checkpointing
         pip install tables
diff --git a/.github/workflows/test-for-pandas-3.yml b/.github/workflows/test-for-pandas-3.yml
deleted file mode 100644
index 8a968d90..00000000
--- a/.github/workflows/test-for-pandas-3.yml
+++ /dev/null
@@ -1,103 +0,0 @@
-name: Run unit tests with Pandas 3.0
-
-on: 
-  push:
-    branches: [ develop ]
-  pull_request:
-    branches: [ develop, releases/** ]
-
-jobs:
-  run_pandas_3_unit_tests:
-
-    runs-on: ubuntu-latest
-    
-    steps:
-    - uses: actions/checkout@v4
-    
-    - name: Set up Python 3.11 
-      uses: actions/setup-python@v5
-      with:
-        python-version: 3.11 
-    
-    - name: Set up Node
-      uses: actions/setup-node@v4
-      with:
-        node-version: node
-    
-    - name: Install Pandas 3.0 (i.e., Nightly)
-      run: |
-        python3 -m pip install --upgrade pip pytest
-        python3 -m pip install --upgrade --force-reinstall --pre --extra-index \
-            https://pypi.anaconda.org/scientific-python-nightly-wheels/simple pandas[pyarrow]
-    
-    - name: Install Python 3 dependencies
-      run: | 
-        python3 -m pip install git+https://github.com/LLNL/hatchet.git@develop
-        python3 -m pip install -r requirements.txt
-        python3 -m pip install .
-        python3 -m pip list
-    
-    # From here on, I'm making use of the trick from the following
-    # GitHub repo to not fail the full CI/CD pipeline if this fails:
-    # https://github.com/burningmantech/ranger-ims-server/pull/1347/files
-    - name: Basic test with PyTest
-      id: test
-      run: |
-        set +e  # Prevent immediate exit if pytest reports errors
-        PYTHONPATH=. $(which pytest)
-        status=$?
-        if [ ${status} -ne 0 ]; then
-            echo "=================================="
-            echo "Unit tests with Pandas 3.0 FAILED!"
-            echo "=================================="
-            echo "::warning::Optional CI test with Pandas 3 failed"
-            echo "optional_fail=true" >> "${GITHUB_OUTPUT}"
-            echo "optional_fail_status=${status}" >> "${GITHUB_OUTPUT}"
-        else
-            echo "=================================="
-            echo "Unit tests with Pandas 3.0 PASSED!"
-            echo "=================================="
-        fi
-        exit 0
-    
-    - name: Find issue comment
-      if: github.event_name == 'push'
-      uses: peter-evans/find-comment@v3
-      id: find_comment
-      with:
-        issue-number: 221
-        comment-author: 'github-actions[bot]'
-        body-includes: Pandas 3.0 Unit Tests
-
-    - name: Add comment to PR if test failed
-      if: github.event_name == 'push' && steps.test.outputs.optional_fail == 'true'
-      uses: peter-evans/create-or-update-comment@v4
-      with:
-        comment-id: ${{ steps.find_comment.outputs.comment-id }}
-        issue-number: 221
-        body: |
-          ### Pandas 3.0 Unit Tests Failed!
-          Due to breaking changes in Pandas 3.0 (namely, copy-on-write), we are performing optional
-          tests of Thicket against the nightly release of Pandas 3.0.
-
-          This is not a full testing failure, and it will not prevent your PR from being merged
-          at this time. However, as we prepare for the release of Pandas 3.0, we encourage
-          all developers to design their code to work with Pandas 3.0, if possible.
-          
-          Pytest status code: ${{ steps.test.outputs.optional_fail_status }}
-          Action log: ${{ github.server_url }}/${{ github.repository }}/actions/runs/${{ github.run_id }}
-        edit-mode: 'replace'
-
-    - name: Add comment to PR if test passed
-      if: github.event_name == 'push' && steps.test.outputs.optional_fail == 'false'
-      uses: peter-evans/create-or-update-comment@v4
-      with:
-        comment-id: ${{ steps.find_comment.outputs.comment-id }}
-        issue-number: 221
-        body: |
-          ### Pandas 3.0 Unit Tests Passed!
-          Due to breaking changes in Pandas 3.0 (namely, copy-on-write), we are performing optional
-          tests of Thicket against the nightly release of Pandas 3.0.
-
-          This PR passed unit tests when run with Pandas 3.0!
-        edit-mode: 'replace'
\ No newline at end of file
diff --git a/.github/workflows/unit-tests.yaml b/.github/workflows/unit-tests.yaml
index cc4c30f4..205bb839 100644
--- a/.github/workflows/unit-tests.yaml
+++ b/.github/workflows/unit-tests.yaml
@@ -14,7 +14,7 @@ jobs:
       matrix:
         # TODO: add macos-latest
         os: [ubuntu-latest]
-        python-version: [3.9, "3.10", "3.11", "3.12", "3.13"]
+        python-version: ["3.11", "3.12", "3.13"]
 
     steps:
     - uses: actions/checkout@v2
@@ -30,35 +30,26 @@ jobs:
         # Obtain the latest cached version of npm (either local cache or action cache)
         node-version: node
 
-    - name: Install Python3 dependencies for 3.12 and later
-      if: ${{ matrix.python-version >= '3.12' }}
+    - name: Install Python3 Dependencies
       run: |
         python -m pip install --upgrade pip pytest setuptools wheel
+        # Install hatchet first: pandas 3 support is on its develop branch, and
+        # requirements.txt cannot be resolved against released hatchet versions
+        python -m pip install git+https://github.com/LLNL/hatchet.git@develop
         pip install -r requirements.txt
-        python -m pip install --upgrade --force-reinstall git+https://github.com/LLNL/hatchet.git@develop
-        python setup.py install
-        python setup.py build_ext --inplace
-        python -m pip list
-
-    - name: Standard Install Python3 Dependencies
-      run: |
-        python -m pip install --upgrade pip pytest
-        pip install -r requirements.txt
-        python -m pip install --upgrade --force-reinstall git+https://github.com/LLNL/hatchet.git@develop
         python setup.py install
         python setup.py build_ext --inplace
         python -m pip list
 
     - name: Update Black
-      if: ${{ matrix.python-version == 3.9 }}
+      if: ${{ matrix.python-version == '3.11' }}
       run: |
-        pip install flake8-pytest-importorskip
         pip install --upgrade click==8.0.4
         pip install black==21.12b0
-        pip install flake8==4.0.1
+        pip install flake8==7.3.0
 
     - name: Lint and Format Check with Flake8 and Black
-      if: ${{ matrix.python-version == 3.9 }}
+      if: ${{ matrix.python-version == '3.11' }}
       run: |
         black --diff --check .
         flake8
@@ -67,10 +58,8 @@ jobs:
       run: |
         python license.py verify
 
-    - name: Install ExtraP
-      if: ${{ matrix.python-version >= 3.7 }}
-      run: |
-        pip install extrap
+    # Extra-P is not installed: all Extra-P releases require NumPy < 2, which
+    # pandas 3 does not support, so thicket/tests/test_model_extrap.py is skipped.
 
     - name: Basic Test with pytest
       run: |
diff --git a/docs/getting_started.rst b/docs/getting_started.rst
index 758479d0..7454d37d 100644
--- a/docs/getting_started.rst
+++ b/docs/getting_started.rst
@@ -15,9 +15,9 @@
 Thicket has the following minimum requirements, which must be installed before Thicket
 is run:
 
-#. Python 3 (3.6 - 3.11)
+#. Python 3 (3.11 or newer)
 #. hatchet
-#. pandas >= 1.1
+#. pandas >= 3.0
 #. numpy
 #. matplotlib, and
 #. scipy
diff --git a/requirements.txt b/requirements.txt
index 3215a087..37e60ac5 100644
--- a/requirements.txt
+++ b/requirements.txt
@@ -1,7 +1,7 @@
 scipy
 numpy
-pandas>=1.1,<3
-llnl-hatchet
+pandas>=3
+llnl-hatchet>=2026.2.0
 matplotlib
 seaborn
 beautifulsoup4
diff --git a/setup.py b/setup.py
index 095c0f90..6b4a4efe 100644
--- a/setup.py
+++ b/setup.py
@@ -33,7 +33,7 @@ def readme():
         "Source Code": "https://github.com/LLNL/thicket",
         "Documentation": "https://thicket.readthedocs.io/",
     },
-    python_requires=">=3.6.1",
+    python_requires=">=3.11",
     packages=[
         "thicket",
         "thicket.stats",
@@ -44,8 +44,11 @@ def readme():
     install_requires=[
         "scipy",
         "numpy",
-        "pandas >= 1.1,<3",
-        "llnl-hatchet",
+        "pandas >= 3",
+        # pandas 3 support is on hatchet develop (LLNL/hatchet#185). The 2026.2.0
+        # release itself requires pandas<3, so this floor admits hatchet develop and
+        # later releases while excluding older releases that lack a pandas pin.
+        "llnl-hatchet >= 2026.2.0",
         "tqdm",
         "more-itertools",
     ],
diff --git a/thicket/ensemble.py b/thicket/ensemble.py
index 08dee352..f79caabe 100644
--- a/thicket/ensemble.py
+++ b/thicket/ensemble.py
@@ -154,7 +154,9 @@ def _handle_metadata():
             # Update index to reflect performance data table index
             if metadata_key != inner_idx:
                 for i in range(len(thickets_cp)):
-                    thickets_cp[i].metadata.reset_index(drop=True, inplace=True)
+                    thickets_cp[i].metadata = thickets_cp[i].metadata.reset_index(
+                        drop=True
+                    )
             if metadata_key is None:
                 for i in range(len(thickets_cp)):
                     thickets_cp[i].metadata.index.set_names(
@@ -163,8 +165,10 @@ def _handle_metadata():
             else:
                 for i in range(len(thickets_cp)):
                     if metadata_key != inner_idx:
-                        thickets_cp[i].metadata.set_index(metadata_key, inplace=True)
-                    thickets_cp[i].metadata.sort_index(inplace=True)
+                        thickets_cp[i].metadata = thickets_cp[i].metadata.set_index(
+                            metadata_key
+                        )
+                    thickets_cp[i].metadata = thickets_cp[i].metadata.sort_index()
 
             # Create multi-index columns
             for i in range(len(thickets_cp)):
@@ -215,7 +219,9 @@ def _handle_perfdata():
                     thickets_cp[i].metadata_columns_to_perfdata(
                         "new_profiles", drop=True
                     )
-                    thickets_cp[i].dataframe.reset_index(level=inner_idx, inplace=True)
+                    thickets_cp[i].dataframe = thickets_cp[i].dataframe.reset_index(
+                        level=inner_idx
+                    )
                     new_mappings.update(
                         pd.Series(
                             thickets_cp[i]
@@ -225,9 +231,11 @@ def _handle_perfdata():
                             index=thickets_cp[i].dataframe[inner_idx],
                         ).to_dict()
                     )
-                    thickets_cp[i].dataframe.drop(inner_idx, axis=1, inplace=True)
-                    thickets_cp[i].dataframe.set_index(
-                        "new_profiles", append=True, inplace=True
+                    thickets_cp[i].dataframe = thickets_cp[i].dataframe.drop(
+                        inner_idx, axis=1
+                    )
+                    thickets_cp[i].dataframe = thickets_cp[i].dataframe.set_index(
+                        "new_profiles", append=True
                     )
                     thickets_cp[i].dataframe.index.rename(
                         thickets_cp[i].profile_idx_name,
@@ -238,7 +246,9 @@ def _handle_perfdata():
                 for i in range(len(thickets_cp)):
                     if metadata_key not in thickets_cp[i].dataframe.index.names:
                         thickets_cp[i].metadata_columns_to_perfdata(metadata_key)
-                    thickets_cp[i].dataframe.reset_index(level=inner_idx, inplace=True)
+                    thickets_cp[i].dataframe = thickets_cp[i].dataframe.reset_index(
+                        level=inner_idx
+                    )
                     new_mappings.update(
                         pd.Series(
                             thickets_cp[i]
@@ -249,11 +259,13 @@ def _handle_perfdata():
                         ).to_dict()
                     )
                     if inner_idx != metadata_key:
-                        thickets_cp[i].dataframe.drop(inner_idx, axis=1, inplace=True)
-                    thickets_cp[i].dataframe.set_index(
-                        metadata_key, append=True, inplace=True
+                        thickets_cp[i].dataframe = thickets_cp[i].dataframe.drop(
+                            inner_idx, axis=1
+                        )
+                    thickets_cp[i].dataframe = thickets_cp[i].dataframe.set_index(
+                        metadata_key, append=True
                     )
-                    thickets_cp[i].dataframe.sort_index(inplace=True)
+                    thickets_cp[i].dataframe = thickets_cp[i].dataframe.sort_index()
 
             # Create tuple columns
             new_columns = [
@@ -283,16 +295,16 @@ def _handle_perfdata():
             )
 
             # Extract "name" columns to upper level
-            nodes = list(set(combined_th.dataframe.reset_index()["node"]))
-            for node in nodes:
-                combined_th.dataframe.loc[node, "name"] = node.frame["name"]
-            combined_th.dataframe.drop(
+            combined_th.dataframe[("name", "")] = [
+                node.frame["name"]
+                for node in combined_th.dataframe.index.get_level_values("node")
+            ]
+            combined_th.dataframe = combined_th.dataframe.drop(
                 columns=[(headers[i], "name") for i in range(len(headers))],
-                inplace=True,
             )
 
             # Sort DataFrame
-            combined_th.dataframe.sort_index(inplace=True)
+            combined_th.dataframe = combined_th.dataframe.sort_index()
 
             return new_mappings
 
@@ -397,9 +409,9 @@ def _index(
             unify_df = _fill_perfdata(unify_df)
 
         # Sort PerfData
-        unify_df.sort_index(inplace=True)
+        unify_df = unify_df.sort_index()
         # Sort Metadata
-        unify_metadata.sort_index(inplace=True)
+        unify_metadata = unify_metadata.sort_index()
 
         # Remove duplicates in metrics
         unify_inc_metrics = list(set(unify_inc_metrics))
diff --git a/thicket/external/console.py b/thicket/external/console.py
index 02433a60..336331bb 100644
--- a/thicket/external/console.py
+++ b/thicket/external/console.py
@@ -128,9 +128,9 @@ def render(self, roots, dataframe, **kwargs):
             metric_series = (dataframe.xs(self.rank, level=1))[self.primary_metric]
         else:
             metric_series = dataframe[self.primary_metric]
-        isfinite_mask = np.isfinite(metric_series.values)
+        isfinite_mask = np.isfinite(metric_series.to_numpy())
         filtered_series = pd.Series(
-            metric_series.values[isfinite_mask], metric_series.index[isfinite_mask]
+            metric_series.to_numpy()[isfinite_mask], metric_series.index[isfinite_mask]
         )
 
         self.max_metric = self.max_value if self.max_value else filtered_series.max()
diff --git a/thicket/groupby.py b/thicket/groupby.py
index 7b0fcfcf..2dd20f19 100644
--- a/thicket/groupby.py
+++ b/thicket/groupby.py
@@ -17,7 +17,11 @@ def agg(self, func, disable_tqdm=False):
         """Aggregate the Thickets' PerfData numerical columns in a GroupBy object.
 
         Arguments:
-            func (function, list or dict): See pandas.DataFrame.aggregate API
+            func (function, str, list or dict): See pandas.DataFrame.aggregate API.
+                Functions are passed to pandas as-is, so NumPy functions keep NumPy
+                semantics (e.g., np.var computes the population variance). Use a
+                pandas aggregation name (e.g., "var" or "std") for pandas semantics.
+                Output columns are named "<column>_<function name or string>".
 
         Returns:
             (Thicket): Aggregated Thicket object.
@@ -50,10 +54,14 @@ def _rename_col(total_cols, tname):
                 tcols[col] = col + "_" + tname
             return tcols
 
+        def _func_name(f):
+            """Name used in output columns for an aggregation function or name."""
+            return f if isinstance(f, str) else f.__name__
+
         def _agg_rows(col_series):
             """Aggregate values in 'col_series' into a set to remove duplicates."""
-            if len(col_series) <= 1:
-                return col_series
+            if len(col_series) == 1:
+                return col_series.iloc[0]
             else:
                 if isinstance(col_series.iloc[0], list) or isinstance(
                     col_series.iloc[0], set
@@ -154,7 +162,7 @@ def _compute_agg_df(col_names, functions, _tk, _agg_cols, _perf_indices):
             cols = list(func.keys())
             # force list structure
             funcs = [[f] if type(f) is not list else f for f in func.values()]
-            snames = [[f.__name__ for f in func] for func in funcs]
+            snames = [[_func_name(f) for f in func] for func in funcs]
             agg_df = _compute_agg_df(snames[0], funcs[0], tk_c, [cols[0]], perf_indices)
             for i in range(1, len(funcs)):
                 t_agg_df = _compute_agg_df(
@@ -166,10 +174,10 @@ def _compute_agg_df(col_names, functions, _tk, _agg_cols, _perf_indices):
                 )
         else:
             if isinstance(func, list):
-                snames = [f.__name__ for f in func]
+                snames = [_func_name(f) for f in func]
                 sfuncs = func
             else:  # Assume single value function
-                snames = [func.__name__]
+                snames = [_func_name(func)]
                 sfuncs = [func]
             agg_df = _compute_agg_df(snames, sfuncs, tk_c, agg_cols, perf_indices)
 
diff --git a/thicket/helpers.py b/thicket/helpers.py
index c6e7314c..55f889b4 100644
--- a/thicket/helpers.py
+++ b/thicket/helpers.py
@@ -121,7 +121,7 @@ def _resolve_missing_indicies(th_list):
             if idx not in th.dataframe.index.names:
                 print(f"Resolving '{idx}' in thicket: ({id(th)})")
                 th.dataframe[idx] = 0
-                th.dataframe.set_index(idx, append=True, inplace=True)
+                th.dataframe = th.dataframe.set_index(idx, append=True)
 
 
 def _set_node_ordering(thickets):
diff --git a/thicket/model_extrap.py b/thicket/model_extrap.py
index 9f779fb3..e7818a69 100644
--- a/thicket/model_extrap.py
+++ b/thicket/model_extrap.py
@@ -320,7 +320,7 @@ def componentize_statsframe(self, columns=None):
             columns = [
                 col
                 for col in self.tht.statsframe.dataframe
-                if isinstance(self.tht.statsframe.dataframe[col][0], ModelWrapper)
+                if isinstance(self.tht.statsframe.dataframe[col].iloc[0], ModelWrapper)
             ]
 
         # Error checking
@@ -329,7 +329,7 @@ def componentize_statsframe(self, columns=None):
                 raise ValueError(
                     "column " + c + " is not in the aggregated statistics table."
                 )
-            elif not isinstance(self.tht.statsframe.dataframe[c][0], ModelWrapper):
+            elif not isinstance(self.tht.statsframe.dataframe[c].iloc[0], ModelWrapper):
                 raise TypeError(
                     "column "
                     + c
diff --git a/thicket/ncu.py b/thicket/ncu.py
index e0692ee4..581ff9e4 100644
--- a/thicket/ncu.py
+++ b/thicket/ncu.py
@@ -138,7 +138,7 @@ def _predicate_builder(kernel, is_regex=False):
         if is_regex:
             return (
                 lambda row: row["name"]
-                .apply(lambda x: kernel in x if x is not None else False)
+                .apply(lambda x: kernel in x if isinstance(x, str) else False)
                 .all()
             )
         else:
diff --git a/thicket/stats/calc_boxplot_statistics.py b/thicket/stats/calc_boxplot_statistics.py
index f2e02059..2b8ee385 100644
--- a/thicket/stats/calc_boxplot_statistics.py
+++ b/thicket/stats/calc_boxplot_statistics.py
@@ -102,9 +102,12 @@ def calc_boxplot_statistics(thicket, columns=[], quartiles=[0.25, 0.5, 0.75], **
             else:
                 thicket.statsframe.inc_metrics.extend(list(boxplot_dict.keys()))
 
-            df_box = pd.DataFrame(boxplot_dict)
-            df_box["node"] = thicket.statsframe.dataframe.index.tolist()
-            df_box.set_index("node", inplace=True)
+            df_box = pd.DataFrame(
+                boxplot_dict,
+                index=pd.Index(
+                    thicket.statsframe.dataframe.index.tolist(), name="node"
+                ),
+            )
             thicket.statsframe.dataframe = thicket.statsframe.dataframe.join(df_box)
     # columnar joined thicket object
     else:
@@ -171,9 +174,12 @@ def calc_boxplot_statistics(thicket, columns=[], quartiles=[0.25, 0.5, 0.75], **
             else:
                 thicket.statsframe.inc_metrics.extend(list(create_multi_index.keys()))
 
-            df_box = pd.DataFrame(create_multi_index)
-            df_box["node"] = thicket.statsframe.dataframe.reset_index()["node"].tolist()
-            df_box = df_box.set_index("node")
+            df_box = pd.DataFrame(
+                create_multi_index,
+                index=pd.Index(
+                    thicket.statsframe.dataframe.index.tolist(), name="node"
+                ),
+            )
             thicket.statsframe.dataframe = thicket.statsframe.dataframe.join(df_box)
 
         # sort columns in index
diff --git a/thicket/stats/check_normality.py b/thicket/stats/check_normality.py
index a5166109..74c1ceaa 100644
--- a/thicket/stats/check_normality.py
+++ b/thicket/stats/check_normality.py
@@ -51,7 +51,7 @@ def check_normality(thicket, columns=None):
         for column in columns:
             output_column_names.append(column + "_normality")
             for i in range(0, len(df[column])):
-                pvalue = df[column][i].pvalue
+                pvalue = df[column].iloc[i].pvalue
 
                 if pvalue < 0.05:
                     thicket.statsframe.dataframe.loc[
@@ -62,7 +62,7 @@ def check_normality(thicket, columns=None):
                         df.index[i], column + "_normality"
                     ] = "True"
                 else:
-                    thicket.stataframe.dataframe.loc[
+                    thicket.statsframe.dataframe.loc[
                         df.index[i], column + "_normality"
                     ] = pd.NA
                 # check to see if exclusive metric
@@ -82,7 +82,7 @@ def check_normality(thicket, columns=None):
         for idx, column in columns:
             output_column_names.append((idx, column + "_normality"))
             for i in range(0, len(df[(idx, column)])):
-                pvalue = df[(idx, column)][i].pvalue
+                pvalue = df[(idx, column)].iloc[i].pvalue
 
                 if pvalue < 0.05:
                     thicket.statsframe.dataframe.loc[
diff --git a/thicket/stats/display_violinplot.py b/thicket/stats/display_violinplot.py
index 542e2cf8..90e4ef90 100644
--- a/thicket/stats/display_violinplot.py
+++ b/thicket/stats/display_violinplot.py
@@ -482,7 +482,7 @@ def display_violinplot_thicket(
             filtered_dfs[thicket_idx]["node"] = thicket_name
 
     master_df = pd.concat(filtered_dfs, ignore_index=True)
-    master_df.rename(columns={"node": "thicket"}, inplace=True)
+    master_df = master_df.rename(columns={"node": "thicket"})
 
     graph = sns.violinplot(
         data=master_df,
diff --git a/thicket/stats/distance.py b/thicket/stats/distance.py
index a19dee67..2e733bf7 100644
--- a/thicket/stats/distance.py
+++ b/thicket/stats/distance.py
@@ -12,6 +12,14 @@
 from .stats_utils import cache_stats_op
 
 
+def _column_key(df, name):
+    """Key for a new column 'name' in 'df': for MultiIndex columns, a full-length
+    tuple ('name', '', ...), as pandas no longer pads shorter keys."""
+    if df.columns.nlevels > 1 and not isinstance(name, tuple):
+        return (name,) + ("",) * (df.columns.nlevels - 1)
+    return name
+
+
 def _calc_bhattacharyya(means_1, means_2, stds_1, stds_2, num_nodes):
     results = []
 
@@ -102,11 +110,11 @@ def bhattacharyya_distance(thicket, columns=None, output_column_name=None):
     mean_columns = th.stats.mean(thicket, columns)
     std_columns = th.stats.std(thicket, columns)
 
-    means_target1 = thicket.statsframe.dataframe[mean_columns[0]]
-    means_target2 = thicket.statsframe.dataframe[mean_columns[1]]
+    means_target1 = thicket.statsframe.dataframe[mean_columns[0]].to_numpy()
+    means_target2 = thicket.statsframe.dataframe[mean_columns[1]].to_numpy()
 
-    stds_target1 = thicket.statsframe.dataframe[std_columns[0]]
-    stds_target2 = thicket.statsframe.dataframe[std_columns[1]]
+    stds_target1 = thicket.statsframe.dataframe[std_columns[0]].to_numpy()
+    stds_target2 = thicket.statsframe.dataframe[std_columns[1]].to_numpy()
 
     result = _calc_bhattacharyya(
         means_target1, means_target2, stds_target1, stds_target2, num_nodes
@@ -130,7 +138,9 @@ def bhattacharyya_distance(thicket, columns=None, output_column_name=None):
     else:
         stats_frame_column_name = output_column_name
 
-    thicket.statsframe.dataframe[stats_frame_column_name] = result
+    thicket.statsframe.dataframe[
+        _column_key(thicket.statsframe.dataframe, stats_frame_column_name)
+    ] = result
     thicket.statsframe.dataframe = thicket.statsframe.dataframe.sort_index(axis=1)
 
     return [stats_frame_column_name]
@@ -186,11 +196,11 @@ def hellinger_distance(thicket, columns=None, output_column_name=None):
     mean_columns = th.stats.mean(thicket, columns)
     std_columns = th.stats.std(thicket, columns)
 
-    means_target1 = thicket.statsframe.dataframe[mean_columns[0]]
-    means_target2 = thicket.statsframe.dataframe[mean_columns[1]]
+    means_target1 = thicket.statsframe.dataframe[mean_columns[0]].to_numpy()
+    means_target2 = thicket.statsframe.dataframe[mean_columns[1]].to_numpy()
 
-    stds_target1 = thicket.statsframe.dataframe[std_columns[0]]
-    stds_target2 = thicket.statsframe.dataframe[std_columns[1]]
+    stds_target1 = thicket.statsframe.dataframe[std_columns[0]].to_numpy()
+    stds_target2 = thicket.statsframe.dataframe[std_columns[1]].to_numpy()
 
     result = _calc_hellinger(
         means_target1, means_target2, stds_target1, stds_target2, num_nodes
@@ -214,7 +224,9 @@ def hellinger_distance(thicket, columns=None, output_column_name=None):
     else:
         stats_frame_column_name = output_column_name
 
-    thicket.statsframe.dataframe[stats_frame_column_name] = result
+    thicket.statsframe.dataframe[
+        _column_key(thicket.statsframe.dataframe, stats_frame_column_name)
+    ] = result
     thicket.statsframe.dataframe = thicket.statsframe.dataframe.sort_index(axis=1)
 
     return [stats_frame_column_name]
diff --git a/thicket/stats/maximum.py b/thicket/stats/maximum.py
index 7c30e408..424ed772 100644
--- a/thicket/stats/maximum.py
+++ b/thicket/stats/maximum.py
@@ -36,7 +36,7 @@ def maximum(thicket, columns=None):
 
     # thicket object without columnar index
     if thicket.dataframe.columns.nlevels == 1:
-        df = thicket.dataframe[columns].reset_index().groupby("node").agg(max)
+        df = thicket.dataframe[columns].reset_index().groupby("node").agg("max")
         for column in columns:
             output_column_names.append(column + "_max")
             thicket.statsframe.dataframe[column + "_max"] = df[column]
@@ -49,7 +49,7 @@ def maximum(thicket, columns=None):
 
     # columnar joined thicket object
     else:
-        df = thicket.dataframe[columns].reset_index(level=1).groupby("node").agg(max)
+        df = thicket.dataframe[columns].reset_index(level=1).groupby("node").agg("max")
         for idx, column in columns:
             output_column_names.append((idx, column + "_max"))
             thicket.statsframe.dataframe[(idx, column + "_max")] = df[(idx, column)]
diff --git a/thicket/stats/mean.py b/thicket/stats/mean.py
index 6f43b155..48ef915e 100644
--- a/thicket/stats/mean.py
+++ b/thicket/stats/mean.py
@@ -3,8 +3,6 @@
 #
 # SPDX-License-Identifier: MIT
 
-import numpy as np
-
 from ..utils import verify_thicket_structures
 from .stats_utils import cache_stats_op
 
@@ -36,7 +34,7 @@ def mean(thicket, columns=None):
 
     # thicket object without columnar index
     if thicket.dataframe.columns.nlevels == 1:
-        df = thicket.dataframe[columns].reset_index().groupby("node").agg(np.mean)
+        df = thicket.dataframe[columns].reset_index().groupby("node").agg("mean")
         for column in columns:
             output_column_names.append(column + "_mean")
             thicket.statsframe.dataframe[column + "_mean"] = df[column]
@@ -48,9 +46,7 @@ def mean(thicket, columns=None):
                 thicket.statsframe.inc_metrics.append(column + "_mean")
     # columnar joined thicket object
     else:
-        df = (
-            thicket.dataframe[columns].reset_index(level=1).groupby("node").agg(np.mean)
-        )
+        df = thicket.dataframe[columns].reset_index(level=1).groupby("node").agg("mean")
         for idx, column in columns:
             output_column_names.append((idx, column + "_mean"))
             thicket.statsframe.dataframe[(idx, column + "_mean")] = df[(idx, column)]
diff --git a/thicket/stats/median.py b/thicket/stats/median.py
index ce9fbd7c..6d8c08bb 100644
--- a/thicket/stats/median.py
+++ b/thicket/stats/median.py
@@ -3,8 +3,6 @@
 #
 # SPDX-License-Identifier: MIT
 
-import numpy as np
-
 from ..utils import verify_thicket_structures
 from .stats_utils import cache_stats_op
 
@@ -36,7 +34,7 @@ def median(thicket, columns=None):
 
     # thicket object without columnar index
     if thicket.dataframe.columns.nlevels == 1:
-        df = thicket.dataframe[columns].reset_index().groupby("node").agg(np.median)
+        df = thicket.dataframe[columns].reset_index().groupby("node").agg("median")
         for column in columns:
             output_column_names.append(column + "_median")
             thicket.statsframe.dataframe[column + "_median"] = df[column]
@@ -53,7 +51,7 @@ def median(thicket, columns=None):
             thicket.dataframe[columns]
             .reset_index(level=1)
             .groupby("node")
-            .agg(np.median)
+            .agg("median")
         )
         for idx, column in columns:
             output_column_names.append(str((idx, column + "_median")))
diff --git a/thicket/stats/minimum.py b/thicket/stats/minimum.py
index 8e4a27bd..c12f6528 100644
--- a/thicket/stats/minimum.py
+++ b/thicket/stats/minimum.py
@@ -36,7 +36,7 @@ def minimum(thicket, columns=None):
 
     # thicket object without columnar index
     if thicket.dataframe.columns.nlevels == 1:
-        df = thicket.dataframe[columns].reset_index().groupby("node").agg(min)
+        df = thicket.dataframe[columns].reset_index().groupby("node").agg("min")
         for column in columns:
             output_column_names.append(column + "_min")
             thicket.statsframe.dataframe[column + "_min"] = df[column]
@@ -48,7 +48,7 @@ def minimum(thicket, columns=None):
                 thicket.statsframe.inc_metrics.append(column + "_min")
     # columnar joined thicket object
     else:
-        df = thicket.dataframe[columns].reset_index(level=1).groupby("node").agg(min)
+        df = thicket.dataframe[columns].reset_index(level=1).groupby("node").agg("min")
         for idx, column in columns:
             output_column_names.append((idx, column + "_min"))
             thicket.statsframe.dataframe[(idx, column + "_min")] = df[(idx, column)]
diff --git a/thicket/stats/percentiles.py b/thicket/stats/percentiles.py
index 22b037f4..8635ecda 100644
--- a/thicket/stats/percentiles.py
+++ b/thicket/stats/percentiles.py
@@ -67,7 +67,7 @@ def percentiles(thicket, columns=None, percentiles=[0.25, 0.50, 0.75]):
         df = df_num.reset_index().groupby("node").quantile(percentiles)
         for column in columns:
             calculated_percentiles = []
-            for node in pd.unique(df.reset_index()["node"].tolist()):
+            for node in pd.unique(df.reset_index()["node"]):
                 calculated_percentiles.append(list(df.loc[node][column]))
 
             for index, percentile in enumerate(percentiles):
@@ -99,7 +99,7 @@ def percentiles(thicket, columns=None, percentiles=[0.25, 0.50, 0.75]):
             calculated_percentiles = []
 
             # Get all the calculated values into a list for each node
-            for node in pd.unique(df.reset_index()["node"].tolist()):
+            for node in pd.unique(df.reset_index()["node"]):
                 calculated_percentiles.append(list(df.loc[node][(idx_level, column)]))
 
             # Go through each of the percentiles, and make them it's own column
diff --git a/thicket/stats/preference.py b/thicket/stats/preference.py
index 91b47d7c..11835ab9 100644
--- a/thicket/stats/preference.py
+++ b/thicket/stats/preference.py
@@ -62,14 +62,14 @@ def preference(thicket, columns, comparison_func, *args, test="ttest", **kwargs)
             if t_statistic < -1 * tvalue or t_statistic > tvalue:
                 pref_mean.append(
                     comparison_func(
-                        thicket.statsframe.dataframe[columns[0] + "_mean"][i],
-                        thicket.statsframe.dataframe[columns[1] + "_mean"][i],
+                        thicket.statsframe.dataframe[columns[0] + "_mean"].iloc[i],
+                        thicket.statsframe.dataframe[columns[1] + "_mean"].iloc[i],
                     )
                 )
                 pref_std.append(
                     comparison_func(
-                        thicket.statsframe.dataframe[columns[0] + "_std"][i],
-                        thicket.statsframe.dataframe[columns[1] + "_std"][i],
+                        thicket.statsframe.dataframe[columns[0] + "_std"].iloc[i],
+                        thicket.statsframe.dataframe[columns[1] + "_std"].iloc[i],
                     )
                 )
             else:
@@ -89,14 +89,14 @@ def preference(thicket, columns, comparison_func, *args, test="ttest", **kwargs)
             if t_statistic < -1 * tvalue or t_statistic > tvalue:
                 pref_mean.append(
                     comparison_func(
-                        thicket.statsframe.dataframe[idx_mean[0]][i],
-                        thicket.statsframe.dataframe[idx_mean[1]][i],
+                        thicket.statsframe.dataframe[idx_mean[0]].iloc[i],
+                        thicket.statsframe.dataframe[idx_mean[1]].iloc[i],
                     )
                 )
                 pref_std.append(
                     comparison_func(
-                        thicket.statsframe.dataframe[idx_std[0]][i],
-                        thicket.statsframe.dataframe[idx_std[1]][i],
+                        thicket.statsframe.dataframe[idx_std[0]].iloc[i],
+                        thicket.statsframe.dataframe[idx_std[1]].iloc[i],
                     )
                 )
             else:
diff --git a/thicket/stats/scoring.py b/thicket/stats/scoring.py
index ce8d0cbd..e0461541 100644
--- a/thicket/stats/scoring.py
+++ b/thicket/stats/scoring.py
@@ -116,11 +116,11 @@ def _score_delta(thicket, columns, output_column_name, scoring_function):
     mean_columns = th.stats.mean(thicket, columns)
     std_columns = th.stats.std(thicket, columns)
 
-    means_target1 = thicket.statsframe.dataframe[mean_columns[0]]
-    means_target2 = thicket.statsframe.dataframe[mean_columns[1]]
+    means_target1 = thicket.statsframe.dataframe[mean_columns[0]].to_numpy()
+    means_target2 = thicket.statsframe.dataframe[mean_columns[1]].to_numpy()
 
-    stds_target1 = thicket.statsframe.dataframe[std_columns[0]]
-    stds_target2 = thicket.statsframe.dataframe[std_columns[1]]
+    stds_target1 = thicket.statsframe.dataframe[std_columns[0]].to_numpy()
+    stds_target2 = thicket.statsframe.dataframe[std_columns[1]].to_numpy()
 
     # Call the scoring function that the user specified
     resulting_scores = scoring_function(
diff --git a/thicket/stats/std.py b/thicket/stats/std.py
index f432c90e..2d890d14 100644
--- a/thicket/stats/std.py
+++ b/thicket/stats/std.py
@@ -3,8 +3,6 @@
 #
 # SPDX-License-Identifier: MIT
 
-import numpy as np
-
 from ..utils import verify_thicket_structures
 from .stats_utils import cache_stats_op
 
@@ -38,7 +36,7 @@ def std(thicket, columns=None):
 
     # thicket object without columnar index
     if thicket.dataframe.columns.nlevels == 1:
-        df = thicket.dataframe[columns].reset_index().groupby("node").agg(np.std)
+        df = thicket.dataframe[columns].reset_index().groupby("node").agg("std")
         for column in columns:
             output_column_names.append(column + "_std")
             thicket.statsframe.dataframe[column + "_std"] = df[column]
@@ -50,7 +48,7 @@ def std(thicket, columns=None):
                 thicket.statsframe.inc_metrics.append(column + "_std")
     # columnar joined thicket object
     else:
-        df = thicket.dataframe[columns].reset_index(level=1).groupby("node").agg(np.std)
+        df = thicket.dataframe[columns].reset_index(level=1).groupby("node").agg("std")
         for idx, column in columns:
             output_column_names.append((idx, column + "_std"))
             thicket.statsframe.dataframe[(idx, column + "_std")] = df[(idx, column)]
diff --git a/thicket/stats/sum.py b/thicket/stats/sum.py
index c6cf5720..14f83d1d 100644
--- a/thicket/stats/sum.py
+++ b/thicket/stats/sum.py
@@ -3,8 +3,6 @@
 #
 # SPDX-License-Identifier: MIT
 
-import numpy as np
-
 from ..utils import verify_thicket_structures
 from .stats_utils import cache_stats_op
 
@@ -36,7 +34,7 @@ def sum(thicket, columns=None):
 
     # thicket object without columnar index
     if thicket.dataframe.columns.nlevels == 1:
-        df = thicket.dataframe[columns].reset_index().groupby("node").agg(np.sum)
+        df = thicket.dataframe[columns].reset_index().groupby("node").agg("sum")
         for column in columns:
             output_column_names.append(column + "_sum")
             thicket.statsframe.dataframe[column + "_sum"] = df[column]
@@ -48,7 +46,7 @@ def sum(thicket, columns=None):
                 thicket.statsframe.inc_metrics.append(column + "_sum")
     # columnar joined thicket object
     else:
-        df = thicket.dataframe[columns].reset_index(level=1).groupby("node").agg(np.sum)
+        df = thicket.dataframe[columns].reset_index(level=1).groupby("node").agg("sum")
         for idx, column in columns:
             output_column_names.append((idx, column + "_sum"))
             thicket.statsframe.dataframe[(idx, column + "_sum")] = df[(idx, column)]
diff --git a/thicket/stats/ttest.py b/thicket/stats/ttest.py
index aadbd760..87f87b70 100644
--- a/thicket/stats/ttest.py
+++ b/thicket/stats/ttest.py
@@ -58,11 +58,11 @@ def __ttest(thicket, columns, alpha=0.05, *args, **kwargs):
         t_statistics = []
         for i in range(0, len(thicket.statsframe.dataframe)):
             tStatistic = ttest_ind_from_stats(
-                mean1=thicket.statsframe.dataframe[mean_columns[0]][i],
-                std1=thicket.statsframe.dataframe[std_columns[0]][i],
+                mean1=thicket.statsframe.dataframe[mean_columns[0]].iloc[i],
+                std1=thicket.statsframe.dataframe[std_columns[0]].iloc[i],
                 nobs1=nobs_column1,
-                mean2=thicket.statsframe.dataframe[mean_columns[1]][i],
-                std2=thicket.statsframe.dataframe[std_columns[1]][i],
+                mean2=thicket.statsframe.dataframe[mean_columns[1]].iloc[i],
+                std2=thicket.statsframe.dataframe[std_columns[1]].iloc[i],
                 nobs2=nobs_column2,
                 equal_var=False,
             )
@@ -96,11 +96,11 @@ def __ttest(thicket, columns, alpha=0.05, *args, **kwargs):
         t_statistics = []
         for i in range(0, len(thicket.statsframe.dataframe)):
             tStatistic = ttest_ind_from_stats(
-                mean1=thicket.statsframe.dataframe[mean_columns[0]][i],
-                std1=thicket.statsframe.dataframe[std_columns[0]][i],
+                mean1=thicket.statsframe.dataframe[mean_columns[0]].iloc[i],
+                std1=thicket.statsframe.dataframe[std_columns[0]].iloc[i],
                 nobs1=nobs_column1,
-                mean2=thicket.statsframe.dataframe[mean_columns[1]][i],
-                std2=thicket.statsframe.dataframe[std_columns[1]][i],
+                mean2=thicket.statsframe.dataframe[mean_columns[1]].iloc[i],
+                std2=thicket.statsframe.dataframe[std_columns[1]].iloc[i],
                 nobs2=nobs_column2,
                 equal_var=False,
             )
diff --git a/thicket/stats/variance.py b/thicket/stats/variance.py
index 73b0c5bf..c1285463 100644
--- a/thicket/stats/variance.py
+++ b/thicket/stats/variance.py
@@ -3,8 +3,6 @@
 #
 # SPDX-License-Identifier: MIT
 
-import numpy as np
-
 from ..utils import verify_thicket_structures
 from .stats_utils import cache_stats_op
 
@@ -39,7 +37,7 @@ def variance(thicket, columns=None):
 
     # thicket object without columnar index
     if thicket.dataframe.columns.nlevels == 1:
-        df = thicket.dataframe[columns].reset_index().groupby("node").agg(np.var)
+        df = thicket.dataframe[columns].reset_index().groupby("node").agg("var")
         for column in columns:
             output_column_names.append(column + "_var")
             thicket.statsframe.dataframe[column + "_var"] = df[column]
@@ -51,7 +49,7 @@ def variance(thicket, columns=None):
                 thicket.statsframe.inc_metrics.append(column + "_var")
     # columnar joined thicket object
     else:
-        df = thicket.dataframe[columns].reset_index(level=1).groupby("node").agg(np.var)
+        df = thicket.dataframe[columns].reset_index(level=1).groupby("node").agg("var")
         for idx, column in columns:
             output_column_names.append((idx, column + "_var"))
             thicket.statsframe.dataframe[(idx, column + "_var")] = df[(idx, column)]
diff --git a/thicket/tests/test_concat_thickets.py b/thicket/tests/test_concat_thickets.py
index a156142f..2432d7c1 100644
--- a/thicket/tests/test_concat_thickets.py
+++ b/thicket/tests/test_concat_thickets.py
@@ -97,8 +97,8 @@ def test_concat_thickets_columns(thicket_axis_columns):
 
     # PerfData and StatsFrame nodes should be in the same order.
     assert (
-        pd.unique(combined_th.dataframe.reset_index()["node"].tolist())
-        == pd.unique(combined_th.statsframe.dataframe.reset_index()["node"].tolist())
+        pd.unique(combined_th.dataframe.reset_index()["node"])
+        == pd.unique(combined_th.statsframe.dataframe.reset_index()["node"])
     ).all()
 
 
diff --git a/thicket/tests/test_copy.py b/thicket/tests/test_copy.py
index 12a4a838..0e974caa 100644
--- a/thicket/tests/test_copy.py
+++ b/thicket/tests/test_copy.py
@@ -57,15 +57,16 @@ def test_copy(rajaperf_seq_O3_1M_cali, intersection, fill_perfdata):
     other.graph.roots[0]._hatchet_nid += 1
     assert self.graph.roots[0]._hatchet_nid == other.graph.roots[0]._hatchet_nid
     assert self.graph.roots[0] is other.graph.roots[0]
-
-    # Shallow copy of data
-    node = other.dataframe.index.get_level_values("node")[0]
-    profile = other.dataframe.index.get_level_values(other.profile_idx_name)[0]
-    other.dataframe.loc[(node, profile), "nid"] = -1
-    assert (
-        other.dataframe.loc[(node, profile), "nid"]
-        == self.dataframe.loc[(node, profile), "nid"]
-    )
+    # Restore the nid: nodes are hashed by nid, and the DataFrame indices hold them
+    other.graph.roots[0]._hatchet_nid -= 1
+
+    # Data is copied lazily (pandas Copy-on-Write): modifying the copy's data does
+    # not change the original's data
+    nid_col = other.dataframe.columns.get_loc("nid")
+    original_value = self.dataframe.iloc[0, nid_col]
+    other.dataframe.iloc[0, nid_col] = -1
+    assert other.dataframe.iloc[0, nid_col] == -1
+    assert self.dataframe.iloc[0, nid_col] == original_value
     # Deep copy of structure
     assert len(self.dataframe.columns) + len(self.dataframe.index[0]) == len(
         other.dataframe.reset_index().columns
diff --git a/thicket/tests/test_groupby.py b/thicket/tests/test_groupby.py
index 20717393..d1925f82 100644
--- a/thicket/tests/test_groupby.py
+++ b/thicket/tests/test_groupby.py
@@ -92,7 +92,11 @@ def test_aggregate(rajaperf_cuda_block128_1M_cali, intersection, fill_perfdata):
 
     epsilon = 0.000001
 
-    def _check_values(_tk_agg):
+    # Variance of "Min time/rank" for the "Algorithm" node over the 10 profiles
+    population_var = 1.825684e-08  # NumPy semantics (ddof=0), e.g., np.var
+    sample_var = 2.0285377777777793e-08  # pandas semantics (ddof=1), e.g., "var"
+
+    def _check_values(_tk_agg, expected_var):
         node = [
             node
             for node in _tk_agg.dataframe.index.get_level_values("node")
@@ -108,21 +112,27 @@ def _check_values(_tk_agg):
             for node in _tk_agg.dataframe.index.get_level_values("node")
             if node.frame["name"] == "Algorithm"
         ][0]
-        assert (
-            abs(
-                _tk_agg.dataframe.loc[algorithm_node, 2]["Min time/rank_var"]
-                - 2.0285377777777793e-08
-            )
-            < epsilon
-        )
+        assert _tk_agg.dataframe.loc[algorithm_node, 2][
+            "Min time/rank_var"
+        ] == pytest.approx(expected_var, rel=1e-6)
 
+    # Functions are passed to pandas as-is, so np.var computes the population variance
     tk_agg = gb.agg(
         func={"Min time/rank": [np.mean, np.var], "Total time": np.mean},
         disable_tqdm=True,
     )
-    _check_values(tk_agg)
+    _check_values(tk_agg, population_var)
     tk_agg = gb.agg(func=[np.mean, np.var], disable_tqdm=True)
-    _check_values(tk_agg)
+    _check_values(tk_agg, population_var)
+
+    # pandas aggregation names use pandas semantics and name the columns the same way
+    tk_agg = gb.agg(
+        func={"Min time/rank": ["mean", "var"], "Total time": "mean"},
+        disable_tqdm=True,
+    )
+    _check_values(tk_agg, sample_var)
+    tk_agg = gb.agg(func=["mean", "var"], disable_tqdm=True)
+    _check_values(tk_agg, sample_var)
 
 
 def test_groupby(rajaperf_seq_O3_1M_cali, intersection, fill_perfdata):
diff --git a/thicket/tests/test_query_stats.py b/thicket/tests/test_query_stats.py
index 833b5974..b15d0dba 100644
--- a/thicket/tests/test_query_stats.py
+++ b/thicket/tests/test_query_stats.py
@@ -47,7 +47,7 @@ def check_query(th_x, hnids, query):
     assert len(filt_nodes) == len(match)
     assert all([n.frame in match_frames for n in filt_nodes])
     assert all([n.frame["name"] in match_names for n in filt_nodes])
-    assert sorted(pd.unique(filt_th_df_nodes)) == sorted(pd.unique(filt_nodes))
+    assert sorted(set(filt_th_df_nodes)) == sorted(set(filt_nodes))
     assert sorted(filt_th_df_profiles.unique().to_list()) == sorted(
         th_df_profiles.unique().to_list()
     )
@@ -55,8 +55,8 @@ def check_query(th_x, hnids, query):
     assert len(sframe_nodes) == len(match)
     assert all([n.frame in match_frames for n in sframe_nodes])
     assert all([n.frame["name"] in match_names for n in sframe_nodes])
-    assert all([n in pd.unique(filt_th_df_nodes) for n in sframe_nodes])
-    assert sorted(pd.unique(filt_th_df_nodes)) == sorted(pd.unique(sframe_nodes))
+    assert all([n in set(filt_th_df_nodes) for n in sframe_nodes])
+    assert sorted(set(filt_th_df_nodes)) == sorted(set(sframe_nodes))
 
     check_identity(th_x, filt_th, ["default_metric", "profile_idx_name"])
 
diff --git a/thicket/tests/test_stats.py b/thicket/tests/test_stats.py
index 442d2f2a..5c39940f 100644
--- a/thicket/tests/test_stats.py
+++ b/thicket/tests/test_stats.py
@@ -4,6 +4,7 @@
 # SPDX-License-Identifier: MIT
 
 import math
+import warnings
 
 import numpy as np
 import pytest
@@ -303,6 +304,45 @@ def test_std_columnar_join(thicket_axis_columns):
     assert (idx, "Min time/rank_std") in combined_th.statsframe.show_metric_columns()
 
 
+def test_aggregation_values(
+    rajaperf_cuda_block128_1M_cali, rajaperf_seq_O3_1M_cali, intersection, fill_perfdata
+):
+    """Check the values of the aggregation functions against NumPy.
+
+    Missing values (e.g., rows added by fill_perfdata) are skipped, and std and
+    variance are sample statistics (ddof=1).
+    """
+    th_ens = th.Thicket.from_caliperreader(
+        rajaperf_cuda_block128_1M_cali[:3] + rajaperf_seq_O3_1M_cali[:3],
+        intersection=intersection,
+        fill_perfdata=fill_perfdata,
+        disable_tqdm=True,
+    )
+    column = "Min time/rank"
+    stats_funcs = {
+        "mean": (th.stats.mean, np.nanmean),
+        "median": (th.stats.median, np.nanmedian),
+        "sum": (th.stats.sum, np.nansum),
+        "max": (th.stats.maximum, np.nanmax),
+        "min": (th.stats.minimum, np.nanmin),
+        "std": (th.stats.std, lambda x: np.nanstd(x, ddof=1)),
+        "var": (th.stats.variance, lambda x: np.nanvar(x, ddof=1)),
+    }
+    for stats_func, _ in stats_funcs.values():
+        stats_func(th_ens, columns=[column])
+
+    for node in th_ens.statsframe.dataframe.index:
+        values = th_ens.dataframe.loc[node, column].to_numpy()
+        for suffix, (_, numpy_func) in stats_funcs.items():
+            with warnings.catch_warnings():
+                # e.g., "Degrees of freedom <= 0" for single values
+                warnings.simplefilter("ignore", RuntimeWarning)
+                expected = numpy_func(values)
+            assert th_ens.statsframe.dataframe.loc[
+                node, f"{column}_{suffix}"
+            ] == pytest.approx(expected, nan_ok=True)
+
+
 def test_percentiles(rajaperf_seq_O3_1M_cali, intersection, fill_perfdata):
     th_ens = th.Thicket.from_caliperreader(
         rajaperf_seq_O3_1M_cali,
@@ -547,7 +587,7 @@ def test_normality(rajaperf_cuda_block128_1M_cali, intersection, fill_perfdata):
 
     th.stats.check_normality(th_ens, columns=["Min time/rank"])
 
-    assert th_ens.statsframe.dataframe["Min time/rank_normality"][0] in {
+    assert th_ens.statsframe.dataframe["Min time/rank_normality"].iloc[0] in {
         "True",
         "False",
     }
@@ -897,11 +937,11 @@ def test_bhattacharyya_distance(thicket_axis_columns):
     mean_columns = th.stats.mean(combined_th_cpy, columns)
     std_columns = th.stats.std(combined_th_cpy, columns)
 
-    means_target1 = combined_th_cpy.statsframe.dataframe[mean_columns[0]]
-    means_target2 = combined_th_cpy.statsframe.dataframe[mean_columns[1]]
+    means_target1 = combined_th_cpy.statsframe.dataframe[mean_columns[0]].to_numpy()
+    means_target2 = combined_th_cpy.statsframe.dataframe[mean_columns[1]].to_numpy()
 
-    stds_target1 = combined_th_cpy.statsframe.dataframe[std_columns[0]]
-    stds_target2 = combined_th_cpy.statsframe.dataframe[std_columns[1]]
+    stds_target1 = combined_th_cpy.statsframe.dataframe[std_columns[0]].to_numpy()
+    stds_target2 = combined_th_cpy.statsframe.dataframe[std_columns[1]].to_numpy()
 
     results = []
 
@@ -997,11 +1037,11 @@ def test_hellinger_distance(thicket_axis_columns):
     mean_columns = th.stats.mean(combined_th_cpy, columns)
     std_columns = th.stats.std(combined_th_cpy, columns)
 
-    means_target1 = combined_th_cpy.statsframe.dataframe[mean_columns[0]]
-    means_target2 = combined_th_cpy.statsframe.dataframe[mean_columns[1]]
+    means_target1 = combined_th_cpy.statsframe.dataframe[mean_columns[0]].to_numpy()
+    means_target2 = combined_th_cpy.statsframe.dataframe[mean_columns[1]].to_numpy()
 
-    stds_target1 = combined_th_cpy.statsframe.dataframe[std_columns[0]]
-    stds_target2 = combined_th_cpy.statsframe.dataframe[std_columns[1]]
+    stds_target1 = combined_th_cpy.statsframe.dataframe[std_columns[0]].to_numpy()
+    stds_target2 = combined_th_cpy.statsframe.dataframe[std_columns[1]].to_numpy()
 
     results = []
 
@@ -1268,10 +1308,10 @@ def test_score_bhattacharyya(thicket_axis_columns):
     )
 
     # NaN, and Inf are not comparable, removed them
-    combined_th_cpy.statsframe.dataframe.replace(
-        [np.inf, -np.inf, np.nan], -1, inplace=True
+    combined_th_cpy.statsframe.dataframe = combined_th_cpy.statsframe.dataframe.replace(
+        [np.inf, -np.inf, np.nan], -1
     )
-    distance_data.replace([np.inf, -np.inf, np.nan], -1, inplace=True)
+    distance_data = distance_data.replace([np.inf, -np.inf, np.nan], -1)
 
     assert (
         combined_th_cpy.statsframe.dataframe[stats_out[0]].tolist()
@@ -1319,10 +1359,10 @@ def test_score_hellinger(thicket_axis_columns):
     )
 
     # NaN, and Inf are not comparable, removed them
-    combined_th_cpy.statsframe.dataframe.replace(
-        [np.inf, -np.inf, np.nan], -1, inplace=True
+    combined_th_cpy.statsframe.dataframe = combined_th_cpy.statsframe.dataframe.replace(
+        [np.inf, -np.inf, np.nan], -1
     )
-    distance_data.replace([np.inf, -np.inf, np.nan], -1, inplace=True)
+    distance_data = distance_data.replace([np.inf, -np.inf, np.nan], -1)
 
     assert (
         combined_th_cpy.statsframe.dataframe[stats_out[0]].tolist()
@@ -1417,4 +1457,4 @@ def test_confidence_interval(thicket_axis_columns):
 
     for col in correct_data.keys():
         for idx, val in enumerate(correct_data[col]):
-            assert combined_th.statsframe.dataframe[col][idx] == val
+            assert combined_th.statsframe.dataframe[col].iloc[idx] == val
diff --git a/thicket/tests/test_thicket.py b/thicket/tests/test_thicket.py
index aaff1bda..ab52524e 100644
--- a/thicket/tests/test_thicket.py
+++ b/thicket/tests/test_thicket.py
@@ -198,9 +198,10 @@ def test_thicketize_graphframe(rajaperf_seq_O3_1M_cali):
     assert ht1.graph == th1.graph
 
     # Check dataframes are equivalent when profile level is dropped
-    th1.dataframe.reset_index(level=th1.profile_idx_name, inplace=True)
-    th1.dataframe.drop(th1.profile_idx_name, axis=1, inplace=True)
-    assert ht1.dataframe.equals(th1.dataframe)
+    th1_df = th1.dataframe.reset_index(level=th1.profile_idx_name).drop(
+        th1.profile_idx_name, axis=1
+    )
+    assert ht1.dataframe.equals(th1_df)
 
 
 def test_unique_metadata_base_cuda(
diff --git a/thicket/tests/test_timeseries.py b/thicket/tests/test_timeseries.py
index 9339e146..fc5c0349 100644
--- a/thicket/tests/test_timeseries.py
+++ b/thicket/tests/test_timeseries.py
@@ -3,9 +3,13 @@
 #
 # SPDX-License-Identifier: MIT
 
+import copy
 import re
 
+import caliperreader as cr
 import hatchet as ht
+import numpy as np
+import pytest
 
 import thicket as tt
 
@@ -128,3 +132,39 @@ def test_timeseries_temporal_pattern(mem_power_timeseries):
     # Check if tree output is correct.
     assert bool(re.search("0.000", mem_tree))
     assert bool(re.search("0.013", pow_tree))
+
+
+@pytest.mark.parametrize("level", ["loop.start_iteration", "iteration#lulesh.cycle"])
+def test_timeseries_missing_observations(mem_power_timeseries, level):
+    """Missing observations stay NaN, including before a node's first sample."""
+    reader = cr.CaliperReader()
+    reader.read(mem_power_timeseries)
+    records = []
+    for record in reader.records:
+        record = record.copy()
+        iteration = record.get("loop.start_iteration")
+        if iteration is not None:
+            record[level] = record.pop("loop.start_iteration")
+        if iteration == "40":
+            record["memstat.vmrss"] = "0"
+        records.append(record)
+        if iteration in ("20", "80"):
+            optional = record.copy()
+            optional["path"] = ["main", "optional"]
+            optional["memstat.vmrss"] = "0" if iteration == "20" else "743"
+            records.append(optional)
+    reader.records = records
+    original_records = copy.deepcopy(records)
+
+    th = tt.Thicket.from_timeseries(reader, level=level)
+    nodes = {node.frame["name"]: node for node in th.graph.traverse()}
+    expected_values = {
+        "main": [679, np.nan, np.nan, np.nan, np.nan],
+        "lulesh.cycle": [743, 743, 0, 743, 743],
+        "optional": [np.nan, 0, np.nan, np.nan, 743],
+    }
+    for name, expected in expected_values.items():
+        values = th.dataframe.xs(nodes[name], level="node")["memstat.vmrss"]
+        assert list(values.index) == [0, 20, 40, 60, 80]
+        np.testing.assert_equal(values.to_numpy(), expected)
+    assert reader.records == original_records
diff --git a/thicket/thicket.py b/thicket/thicket.py
index 95e77f64..9da374d3 100644
--- a/thicket/thicket.py
+++ b/thicket/thicket.py
@@ -234,12 +234,15 @@ def thicketize_graphframe(gf, prf):
             th.metadata = pd.DataFrame.from_dict(temp_meta, orient="index")
             th.metadata.index.set_names(th.profile_idx_name, inplace=True)
 
-            # Add profile to dataframe index
-            th.dataframe[th.profile_idx_name] = hash_arg
+            # Add profile to dataframe index. Build a new DataFrame rather than
+            # modifying gf.dataframe, which th.dataframe refers to.
             index_names = list(th.dataframe.index.names)
             index_names.insert(1, th.profile_idx_name)
-            th.dataframe.reset_index(inplace=True)
-            th.dataframe.set_index(index_names, inplace=True)
+            th.dataframe = (
+                th.dataframe.assign(**{th.profile_idx_name: hash_arg})
+                .reset_index()
+                .set_index(index_names)
+            )
 
         return th
 
@@ -271,12 +274,14 @@ def thicketize_timeseries_graphframe(gf, iter_column="loop.start_iteration"):
 
         # Use iter_column to set iteration dataframe index then delete it
         iteration = int(th.dataframe[iter_column].dropna().max())
-        th.dataframe.drop(columns=[iter_column], inplace=True)
-        th.dataframe["iteration"] = iteration
         index_names = list(th.dataframe.index.names)
         index_names.insert(1, "iteration")
-        th.dataframe.reset_index(inplace=True)
-        th.dataframe.set_index(index_names, inplace=True)
+        th.dataframe = (
+            th.dataframe.drop(columns=[iter_column])
+            .assign(iteration=iteration)
+            .reset_index()
+            .set_index(index_names)
+        )
 
         return th
 
@@ -585,15 +590,17 @@ def from_timeseries(
         """Read in a single Caliper .cali file that has timeseries records
         gets split into a list of graphframes by hatchet.
 
+        Nodes without an observation in a time step have NaN metric values.
+
         Arguments:
-            filename_or_stream (str or file-like): name of a Caliper timeseries output file in
-                `.cali` format, or an open file object to read one
+            filename_or_caliperreader (str or CaliperReader): name of a Caliper
+                timeseries output file in `.cali` format, or a CaliperReader object
             level (str): level to split the timeseries, default "loop.start_iteration"
             intersection (bool): whether to perform intersection or union (default)
         """
         # if we are reading a timeseries file we will expect a list back from the reader
         ens_list = []
-        gf_list = GraphFrame.from_timeseries(filename_or_caliperreader)
+        gf_list = GraphFrame.from_timeseries(filename_or_caliperreader, level=level)
         for gf in gf_list:
             ens_list.append(Thicket.thicketize_timeseries_graphframe(gf, level))
 
@@ -912,7 +919,7 @@ def from_json(json_thicket):
 
         if "metadata" in thicket_dict:
             mf = pd.DataFrame(thicket_dict["metadata"])
-            mf.set_index(mf[th.profile_idx_name], inplace=True)
+            mf = mf.set_index(mf[th.profile_idx_name])
             if th.profile_idx_name in mf.columns:
                 mf = mf.drop(columns=[th.profile_idx_name])
             th.metadata = mf
@@ -1062,7 +1069,7 @@ def metadata_columns_to_perfdata(
             if mkey in self.dataframe.columns:
                 # Drop column to overwrite, otherwise warn and return
                 if overwrite:
-                    self.dataframe.drop(mkey, axis=1, inplace=True)
+                    self.dataframe = self.dataframe.drop(mkey, axis=1)
                 else:
                     warnings.warn(
                         "Column "
@@ -1078,7 +1085,7 @@ def metadata_columns_to_perfdata(
 
         # Drop column
         if drop:
-            self.metadata.drop(metadata_columns, axis=1, inplace=True)
+            self.metadata = self.metadata.drop(metadata_columns, axis=1)
 
     def squash(self, update_inc_cols=True, new_statsframe=True):
         """Rewrite the Graph to include only nodes present in the performance
@@ -1099,10 +1106,10 @@ def squash(self, update_inc_cols=True, new_statsframe=True):
         # Hatchet's squash code
         #####
         index_names = self.dataframe.index.names
-        self.dataframe.reset_index(inplace=True)
+        df = self.dataframe.reset_index()
 
         # create new nodes for each unique node in the old dataframe
-        old_to_new = {n: n.copy() for n in set(self.dataframe["node"])}
+        old_to_new = {n: n.copy() for n in set(df["node"])}
         for i in old_to_new:
             old_to_new[i]._hatchet_nid = i._hatchet_nid
 
@@ -1155,7 +1162,6 @@ def rewire(node, new_parent, visited):
         graph.enumerate_traverse()
 
         # reindex new dataframe with new nodes
-        df = self.dataframe.copy()
         df["node"] = df["node"].apply(lambda x: old_to_new[x])
 
         # at this point, the graph is potentially invalid, as some nodes
@@ -1163,8 +1169,7 @@ def rewire(node, new_parent, visited):
         merges = graph.normalize()
         df["node"] = df["node"].apply(lambda n: merges.get(n, n))
 
-        self.dataframe.set_index(index_names, inplace=True)
-        df.set_index(index_names, inplace=True)
+        df = df.set_index(index_names)
         # create dict that stores aggregation function for each column
         agg_dict = {}
         for col in df.columns.tolist():
@@ -1178,8 +1183,7 @@ def rewire(node, new_parent, visited):
                 agg_dict[col] = lambda x: x.iloc[0]
 
         # perform a groupby to merge nodes with the same callpath
-        agg_df = df.groupby(index_names).agg(agg_dict)
-        agg_df.sort_index(inplace=True)
+        agg_df = df.groupby(index_names).agg(agg_dict).sort_index()
 
         #####
         # Thicket's squash code
@@ -1230,7 +1234,10 @@ def rewire(node, new_parent, visited):
     def copy(self):
         """Return a partially shallow copy of the Thicket.
 
-        See GraphFrame.copy() for more details
+        See GraphFrame.copy() for more details. The graph, and the Node objects in
+        the DataFrame indices, are shared with self. The DataFrames are shallow
+        copies, which pandas Copy-on-Write copies lazily: modifying the data of
+        one Thicket never changes the data of the other.
 
         Arguments:
             self (Thicket): object to make a copy of
@@ -1503,10 +1510,9 @@ def from_statsframes(tk_list, metadata_key=None, disable_tqdm=False):
             tk_copy.graph = tk_copy.statsframe.graph
 
             # Modify the performance data table
-            df = tk_copy.statsframe.dataframe
+            df = tk_copy.statsframe.dataframe.copy()
             df[idx_name] = tk_id
-            df.set_index(idx_name, inplace=True, append=True)
-            tk_copy.dataframe = df
+            tk_copy.dataframe = df.set_index(idx_name, append=True)
 
             # Adjust profile and profile_mapping
             tk_copy.profile = [tk_id]
@@ -1515,12 +1521,12 @@ def from_statsframes(tk_list, metadata_key=None, disable_tqdm=False):
 
             # Modify metadata dataframe
             tk_copy.metadata[idx_name] = tk_id
-            tk_copy.metadata.set_index(idx_name, inplace=True)
+            tk_copy.metadata = tk_copy.metadata.set_index(idx_name)
 
             def _agg_to_set(obj):
                 """Aggregate values in 'obj' into a set to remove duplicates."""
-                if len(obj) <= 1:
-                    return obj
+                if len(obj) == 1:
+                    return obj.iloc[0]
                 else:
                     if isinstance(obj.iloc[0], list) or isinstance(obj.iloc[0], set):
                         _set = set(tuple(i) for i in obj)
@@ -1617,11 +1623,9 @@ def intersection(self):
             return tkc
         else:
             # Check for padded perfdata
-            if self.dataframe["name"].isnull().any():
-                # Row that didn't exist will contain "None" in the name column.
-                query = Query().match(
-                    ".", lambda df: df["name"].apply(lambda n: n is not None).all()
-                )
+            if self.dataframe["name"].isna().any():
+                # Rows that didn't exist contain a missing value in the name column.
+                query = Query().match(".", lambda df: df["name"].notna().all())
             else:
                 # If perfdata not padded
                 query = Query().match(".", lambda df: len(df) == len(self.profile))
@@ -1753,9 +1757,8 @@ def query(
                     type(query_obj)
                 )
             )
-        dframe_copy = self.dataframe.copy()
         index_names = self.dataframe.index.names
-        dframe_copy.reset_index(inplace=True)
+        dframe_copy = self.dataframe.reset_index()
         query = (
             local_query_obj
             if not is_old_style_query(query_obj)
@@ -1765,7 +1768,7 @@ def query(
         filtered_df = dframe_copy.loc[dframe_copy["node"].isin(query_matches)]
         if filtered_df.shape[0] == 0:
             raise EmptyQuery("The provided query would have produced an empty Thicket.")
-        filtered_df.set_index(index_names, inplace=True)
+        filtered_df = filtered_df.set_index(index_names)
 
         filtered_th = self.deepcopy()
         filtered_th.dataframe = filtered_df
@@ -1799,9 +1802,8 @@ def query_stats(self, query_obj, squash=True, update_inc_cols=True):
                     type(query_obj)
                 )
             )
-        sframe_copy = self.statsframe.dataframe.copy()
         sf_index_names = self.statsframe.dataframe.index.names
-        sframe_copy.reset_index(inplace=True)
+        sframe_copy = self.statsframe.dataframe.reset_index()
 
         query = (
             local_query_obj
@@ -1817,13 +1819,12 @@ def query_stats(self, query_obj, squash=True, update_inc_cols=True):
             raise EmptyQuery("The provided query would have produced an empty Thicket.")
 
         df_index_names = self.dataframe.index.names
-        dframe_copy = self.dataframe.copy()
-        dframe_copy.reset_index(inplace=True)
+        dframe_copy = self.dataframe.reset_index()
 
         filtered_df = dframe_copy[dframe_copy["node"].isin(query_matches)]
-        filtered_df.set_index(df_index_names, inplace=True)
+        filtered_df = filtered_df.set_index(df_index_names)
 
-        filtered_sf_df.set_index(sf_index_names, inplace=True)
+        filtered_sf_df = filtered_sf_df.set_index(sf_index_names)
 
         filtered_th = Thicket(
             graph=self.graph,
diff --git a/thicket/utils.py b/thicket/utils.py
index 16cef3ac..eeeb2937 100644
--- a/thicket/utils.py
+++ b/thicket/utils.py
@@ -91,12 +91,12 @@ def _check_missing_hnid(df):
             i += 1
 
     def _validate_name_column(df):
-        """Check if all of the values in a node's name column are either its name or None."""
+        """Check if all of the values in a node's name column are either its name or missing."""
         for node in set(df.index.get_level_values("node")):
             names = set(df.loc[node]["name"].tolist())
             node_name = node.frame["name"]
             for name in names:
-                if name != node_name and name is not None:
+                if name != node_name and not pd.isna(name):
                     raise InvalidNameError(
                         f"Value in the Thicket.dataframe's 'name' column is not valid. {name} != {node_name}"
                     )
@@ -262,7 +262,7 @@ def validate_nodes(tk):
 
 
 def _fill_perfdata(df):
-    """Create full index for DataFrame and fill created rows with NaN's or Nones where applicable.
+    """Create full index for DataFrame and fill created rows with missing values (NaN).
 
     Arguments:
         df (DataFrame): DataFrame to fill missing rows in
@@ -274,15 +274,12 @@ def _fill_perfdata(df):
     try:
         # Value used to fill new rows
         fill_value = np.nan
-        # Fill missing rows in dataframe with NaN's
+        # Fill missing rows in dataframe with NaN's. String columns keep their
+        # dtype, which also uses NaN as its missing value.
         new_df = new_df.reindex(
             pd.MultiIndex.from_product(new_df.index.levels),
             fill_value=fill_value,
         )
-        # Replace "NaN" with "None" in columns of string type
-        for col in new_df.columns:
-            if pd.api.types.is_string_dtype(new_df[col].dtype):
-                new_df[col] = new_df[col].replace({fill_value: None})
     except ValueError as e:
         estr = str(e)
         if estr == "cannot handle a non-unique multi-index!":