diff --git a/.Rbuildignore b/.Rbuildignore
new file mode 100644
index 0000000..89b0d98
--- /dev/null
+++ b/.Rbuildignore
@@ -0,0 +1,4 @@
+^.*\.Rproj$
+^\.Rproj\.user$
+^\.github$
+^inst$
diff --git a/.github/.gitignore b/.github/.gitignore
new file mode 100644
index 0000000..2d19fc7
--- /dev/null
+++ b/.github/.gitignore
@@ -0,0 +1 @@
+*.html
diff --git a/.github/workflows/R-CMD-check.yaml b/.github/workflows/R-CMD-check.yaml
new file mode 100644
index 0000000..f4b17a4
--- /dev/null
+++ b/.github/workflows/R-CMD-check.yaml
@@ -0,0 +1,29 @@
+# Workflow derived from https://github.com/r-lib/actions/tree/v2/examples
+# Need help debugging build failures? Start at https://github.com/r-lib/actions#where-to-find-help
+on:
+  push:
+    branches: [main, master]
+  pull_request:
+    branches: [main, master]
+
+name: R-CMD-check
+
+jobs:
+  R-CMD-check:
+    runs-on: ubuntu-latest
+    env:
+      GITHUB_PAT: ${{ secrets.GITHUB_TOKEN }}
+      R_KEEP_PKG_SOURCE: yes
+    steps:
+      - uses: actions/checkout@v3
+
+      - uses: r-lib/actions/setup-r@v2
+        with:
+          use-public-rspm: true
+
+      - uses: r-lib/actions/setup-r-dependencies@v2
+        with:
+          extra-packages: any::rcmdcheck
+          needs: check
+
+      - uses: r-lib/actions/check-r-package@v2
diff --git a/.gitignore b/.gitignore
new file mode 100644
index 0000000..38f60e1
--- /dev/null
+++ b/.gitignore
@@ -0,0 +1,15 @@
+.Rproj.user
+.Rhistory
+.RData
+.Ruserdata
+*.Rproj
+*.csv
+*.pdf
+inst/data/
+inst/ignore/
+mimic
+*.Rout
+*.bib
+*.rds
+*.sh
+
diff --git a/DESCRIPTION b/DESCRIPTION
new file mode 100644
index 0000000..f6d5f87
--- /dev/null
+++ b/DESCRIPTION
@@ -0,0 +1,18 @@
+Package: comet
+Type: Package
+Title: Covariance Measure Tests for Conditional Independence 
+Version: 0.0-1
+Authors@R: person("Lucas", "Kook", email = "lucasheinrich.kook@gmail.com",
+    role = c("aut", "cre"))
+Description: Covariance measure tests for conditional independence testing 
+    against conditional covariance and nonlinear conditional mean alternatives.
+    Contains versions of the generalised covariance measure test (Shah and Peters,
+    2020, <doi:TODO>) and projected covariance measure test (Lundborg et al., 2023,
+    <doi:TODO>). Applications can be found in Kook and Lundborg (2024, <doi:>).
+Imports: mlt, sandwich, ranger, glmnet
+License: GPL-3
+Encoding: UTF-8
+RoxygenNote: 7.2.3
+Suggests: 
+    testthat (>= 3.0.0), tram, ggplot2, tidyr, ggpubr, dplyr
+Config/testthat/edition: 3
diff --git a/LICENSE b/LICENSE
deleted file mode 100644
index f288702..0000000
--- a/LICENSE
+++ /dev/null
@@ -1,674 +0,0 @@
-                    GNU GENERAL PUBLIC LICENSE
-                       Version 3, 29 June 2007
-
- Copyright (C) 2007 Free Software Foundation, Inc. <https://fsf.org/>
- Everyone is permitted to copy and distribute verbatim copies
- of this license document, but changing it is not allowed.
-
-                            Preamble
-
-  The GNU General Public License is a free, copyleft license for
-software and other kinds of works.
-
-  The licenses for most software and other practical works are designed
-to take away your freedom to share and change the works.  By contrast,
-the GNU General Public License is intended to guarantee your freedom to
-share and change all versions of a program--to make sure it remains free
-software for all its users.  We, the Free Software Foundation, use the
-GNU General Public License for most of our software; it applies also to
-any other work released this way by its authors.  You can apply it to
-your programs, too.
-
-  When we speak of free software, we are referring to freedom, not
-price.  Our General Public Licenses are designed to make sure that you
-have the freedom to distribute copies of free software (and charge for
-them if you wish), that you receive source code or can get it if you
-want it, that you can change the software or use pieces of it in new
-free programs, and that you know you can do these things.
-
-  To protect your rights, we need to prevent others from denying you
-these rights or asking you to surrender the rights.  Therefore, you have
-certain responsibilities if you distribute copies of the software, or if
-you modify it: responsibilities to respect the freedom of others.
-
-  For example, if you distribute copies of such a program, whether
-gratis or for a fee, you must pass on to the recipients the same
-freedoms that you received.  You must make sure that they, too, receive
-or can get the source code.  And you must show them these terms so they
-know their rights.
-
-  Developers that use the GNU GPL protect your rights with two steps:
-(1) assert copyright on the software, and (2) offer you this License
-giving you legal permission to copy, distribute and/or modify it.
-
-  For the developers' and authors' protection, the GPL clearly explains
-that there is no warranty for this free software.  For both users' and
-authors' sake, the GPL requires that modified versions be marked as
-changed, so that their problems will not be attributed erroneously to
-authors of previous versions.
-
-  Some devices are designed to deny users access to install or run
-modified versions of the software inside them, although the manufacturer
-can do so.  This is fundamentally incompatible with the aim of
-protecting users' freedom to change the software.  The systematic
-pattern of such abuse occurs in the area of products for individuals to
-use, which is precisely where it is most unacceptable.  Therefore, we
-have designed this version of the GPL to prohibit the practice for those
-products.  If such problems arise substantially in other domains, we
-stand ready to extend this provision to those domains in future versions
-of the GPL, as needed to protect the freedom of users.
-
-  Finally, every program is threatened constantly by software patents.
-States should not allow patents to restrict development and use of
-software on general-purpose computers, but in those that do, we wish to
-avoid the special danger that patents applied to a free program could
-make it effectively proprietary.  To prevent this, the GPL assures that
-patents cannot be used to render the program non-free.
-
-  The precise terms and conditions for copying, distribution and
-modification follow.
-
-                       TERMS AND CONDITIONS
-
-  0. Definitions.
-
-  "This License" refers to version 3 of the GNU General Public License.
-
-  "Copyright" also means copyright-like laws that apply to other kinds of
-works, such as semiconductor masks.
-
-  "The Program" refers to any copyrightable work licensed under this
-License.  Each licensee is addressed as "you".  "Licensees" and
-"recipients" may be individuals or organizations.
-
-  To "modify" a work means to copy from or adapt all or part of the work
-in a fashion requiring copyright permission, other than the making of an
-exact copy.  The resulting work is called a "modified version" of the
-earlier work or a work "based on" the earlier work.
-
-  A "covered work" means either the unmodified Program or a work based
-on the Program.
-
-  To "propagate" a work means to do anything with it that, without
-permission, would make you directly or secondarily liable for
-infringement under applicable copyright law, except executing it on a
-computer or modifying a private copy.  Propagation includes copying,
-distribution (with or without modification), making available to the
-public, and in some countries other activities as well.
-
-  To "convey" a work means any kind of propagation that enables other
-parties to make or receive copies.  Mere interaction with a user through
-a computer network, with no transfer of a copy, is not conveying.
-
-  An interactive user interface displays "Appropriate Legal Notices"
-to the extent that it includes a convenient and prominently visible
-feature that (1) displays an appropriate copyright notice, and (2)
-tells the user that there is no warranty for the work (except to the
-extent that warranties are provided), that licensees may convey the
-work under this License, and how to view a copy of this License.  If
-the interface presents a list of user commands or options, such as a
-menu, a prominent item in the list meets this criterion.
-
-  1. Source Code.
-
-  The "source code" for a work means the preferred form of the work
-for making modifications to it.  "Object code" means any non-source
-form of a work.
-
-  A "Standard Interface" means an interface that either is an official
-standard defined by a recognized standards body, or, in the case of
-interfaces specified for a particular programming language, one that
-is widely used among developers working in that language.
-
-  The "System Libraries" of an executable work include anything, other
-than the work as a whole, that (a) is included in the normal form of
-packaging a Major Component, but which is not part of that Major
-Component, and (b) serves only to enable use of the work with that
-Major Component, or to implement a Standard Interface for which an
-implementation is available to the public in source code form.  A
-"Major Component", in this context, means a major essential component
-(kernel, window system, and so on) of the specific operating system
-(if any) on which the executable work runs, or a compiler used to
-produce the work, or an object code interpreter used to run it.
-
-  The "Corresponding Source" for a work in object code form means all
-the source code needed to generate, install, and (for an executable
-work) run the object code and to modify the work, including scripts to
-control those activities.  However, it does not include the work's
-System Libraries, or general-purpose tools or generally available free
-programs which are used unmodified in performing those activities but
-which are not part of the work.  For example, Corresponding Source
-includes interface definition files associated with source files for
-the work, and the source code for shared libraries and dynamically
-linked subprograms that the work is specifically designed to require,
-such as by intimate data communication or control flow between those
-subprograms and other parts of the work.
-
-  The Corresponding Source need not include anything that users
-can regenerate automatically from other parts of the Corresponding
-Source.
-
-  The Corresponding Source for a work in source code form is that
-same work.
-
-  2. Basic Permissions.
-
-  All rights granted under this License are granted for the term of
-copyright on the Program, and are irrevocable provided the stated
-conditions are met.  This License explicitly affirms your unlimited
-permission to run the unmodified Program.  The output from running a
-covered work is covered by this License only if the output, given its
-content, constitutes a covered work.  This License acknowledges your
-rights of fair use or other equivalent, as provided by copyright law.
-
-  You may make, run and propagate covered works that you do not
-convey, without conditions so long as your license otherwise remains
-in force.  You may convey covered works to others for the sole purpose
-of having them make modifications exclusively for you, or provide you
-with facilities for running those works, provided that you comply with
-the terms of this License in conveying all material for which you do
-not control copyright.  Those thus making or running the covered works
-for you must do so exclusively on your behalf, under your direction
-and control, on terms that prohibit them from making any copies of
-your copyrighted material outside their relationship with you.
-
-  Conveying under any other circumstances is permitted solely under
-the conditions stated below.  Sublicensing is not allowed; section 10
-makes it unnecessary.
-
-  3. Protecting Users' Legal Rights From Anti-Circumvention Law.
-
-  No covered work shall be deemed part of an effective technological
-measure under any applicable law fulfilling obligations under article
-11 of the WIPO copyright treaty adopted on 20 December 1996, or
-similar laws prohibiting or restricting circumvention of such
-measures.
-
-  When you convey a covered work, you waive any legal power to forbid
-circumvention of technological measures to the extent such circumvention
-is effected by exercising rights under this License with respect to
-the covered work, and you disclaim any intention to limit operation or
-modification of the work as a means of enforcing, against the work's
-users, your or third parties' legal rights to forbid circumvention of
-technological measures.
-
-  4. Conveying Verbatim Copies.
-
-  You may convey verbatim copies of the Program's source code as you
-receive it, in any medium, provided that you conspicuously and
-appropriately publish on each copy an appropriate copyright notice;
-keep intact all notices stating that this License and any
-non-permissive terms added in accord with section 7 apply to the code;
-keep intact all notices of the absence of any warranty; and give all
-recipients a copy of this License along with the Program.
-
-  You may charge any price or no price for each copy that you convey,
-and you may offer support or warranty protection for a fee.
-
-  5. Conveying Modified Source Versions.
-
-  You may convey a work based on the Program, or the modifications to
-produce it from the Program, in the form of source code under the
-terms of section 4, provided that you also meet all of these conditions:
-
-    a) The work must carry prominent notices stating that you modified
-    it, and giving a relevant date.
-
-    b) The work must carry prominent notices stating that it is
-    released under this License and any conditions added under section
-    7.  This requirement modifies the requirement in section 4 to
-    "keep intact all notices".
-
-    c) You must license the entire work, as a whole, under this
-    License to anyone who comes into possession of a copy.  This
-    License will therefore apply, along with any applicable section 7
-    additional terms, to the whole of the work, and all its parts,
-    regardless of how they are packaged.  This License gives no
-    permission to license the work in any other way, but it does not
-    invalidate such permission if you have separately received it.
-
-    d) If the work has interactive user interfaces, each must display
-    Appropriate Legal Notices; however, if the Program has interactive
-    interfaces that do not display Appropriate Legal Notices, your
-    work need not make them do so.
-
-  A compilation of a covered work with other separate and independent
-works, which are not by their nature extensions of the covered work,
-and which are not combined with it such as to form a larger program,
-in or on a volume of a storage or distribution medium, is called an
-"aggregate" if the compilation and its resulting copyright are not
-used to limit the access or legal rights of the compilation's users
-beyond what the individual works permit.  Inclusion of a covered work
-in an aggregate does not cause this License to apply to the other
-parts of the aggregate.
-
-  6. Conveying Non-Source Forms.
-
-  You may convey a covered work in object code form under the terms
-of sections 4 and 5, provided that you also convey the
-machine-readable Corresponding Source under the terms of this License,
-in one of these ways:
-
-    a) Convey the object code in, or embodied in, a physical product
-    (including a physical distribution medium), accompanied by the
-    Corresponding Source fixed on a durable physical medium
-    customarily used for software interchange.
-
-    b) Convey the object code in, or embodied in, a physical product
-    (including a physical distribution medium), accompanied by a
-    written offer, valid for at least three years and valid for as
-    long as you offer spare parts or customer support for that product
-    model, to give anyone who possesses the object code either (1) a
-    copy of the Corresponding Source for all the software in the
-    product that is covered by this License, on a durable physical
-    medium customarily used for software interchange, for a price no
-    more than your reasonable cost of physically performing this
-    conveying of source, or (2) access to copy the
-    Corresponding Source from a network server at no charge.
-
-    c) Convey individual copies of the object code with a copy of the
-    written offer to provide the Corresponding Source.  This
-    alternative is allowed only occasionally and noncommercially, and
-    only if you received the object code with such an offer, in accord
-    with subsection 6b.
-
-    d) Convey the object code by offering access from a designated
-    place (gratis or for a charge), and offer equivalent access to the
-    Corresponding Source in the same way through the same place at no
-    further charge.  You need not require recipients to copy the
-    Corresponding Source along with the object code.  If the place to
-    copy the object code is a network server, the Corresponding Source
-    may be on a different server (operated by you or a third party)
-    that supports equivalent copying facilities, provided you maintain
-    clear directions next to the object code saying where to find the
-    Corresponding Source.  Regardless of what server hosts the
-    Corresponding Source, you remain obligated to ensure that it is
-    available for as long as needed to satisfy these requirements.
-
-    e) Convey the object code using peer-to-peer transmission, provided
-    you inform other peers where the object code and Corresponding
-    Source of the work are being offered to the general public at no
-    charge under subsection 6d.
-
-  A separable portion of the object code, whose source code is excluded
-from the Corresponding Source as a System Library, need not be
-included in conveying the object code work.
-
-  A "User Product" is either (1) a "consumer product", which means any
-tangible personal property which is normally used for personal, family,
-or household purposes, or (2) anything designed or sold for incorporation
-into a dwelling.  In determining whether a product is a consumer product,
-doubtful cases shall be resolved in favor of coverage.  For a particular
-product received by a particular user, "normally used" refers to a
-typical or common use of that class of product, regardless of the status
-of the particular user or of the way in which the particular user
-actually uses, or expects or is expected to use, the product.  A product
-is a consumer product regardless of whether the product has substantial
-commercial, industrial or non-consumer uses, unless such uses represent
-the only significant mode of use of the product.
-
-  "Installation Information" for a User Product means any methods,
-procedures, authorization keys, or other information required to install
-and execute modified versions of a covered work in that User Product from
-a modified version of its Corresponding Source.  The information must
-suffice to ensure that the continued functioning of the modified object
-code is in no case prevented or interfered with solely because
-modification has been made.
-
-  If you convey an object code work under this section in, or with, or
-specifically for use in, a User Product, and the conveying occurs as
-part of a transaction in which the right of possession and use of the
-User Product is transferred to the recipient in perpetuity or for a
-fixed term (regardless of how the transaction is characterized), the
-Corresponding Source conveyed under this section must be accompanied
-by the Installation Information.  But this requirement does not apply
-if neither you nor any third party retains the ability to install
-modified object code on the User Product (for example, the work has
-been installed in ROM).
-
-  The requirement to provide Installation Information does not include a
-requirement to continue to provide support service, warranty, or updates
-for a work that has been modified or installed by the recipient, or for
-the User Product in which it has been modified or installed.  Access to a
-network may be denied when the modification itself materially and
-adversely affects the operation of the network or violates the rules and
-protocols for communication across the network.
-
-  Corresponding Source conveyed, and Installation Information provided,
-in accord with this section must be in a format that is publicly
-documented (and with an implementation available to the public in
-source code form), and must require no special password or key for
-unpacking, reading or copying.
-
-  7. Additional Terms.
-
-  "Additional permissions" are terms that supplement the terms of this
-License by making exceptions from one or more of its conditions.
-Additional permissions that are applicable to the entire Program shall
-be treated as though they were included in this License, to the extent
-that they are valid under applicable law.  If additional permissions
-apply only to part of the Program, that part may be used separately
-under those permissions, but the entire Program remains governed by
-this License without regard to the additional permissions.
-
-  When you convey a copy of a covered work, you may at your option
-remove any additional permissions from that copy, or from any part of
-it.  (Additional permissions may be written to require their own
-removal in certain cases when you modify the work.)  You may place
-additional permissions on material, added by you to a covered work,
-for which you have or can give appropriate copyright permission.
-
-  Notwithstanding any other provision of this License, for material you
-add to a covered work, you may (if authorized by the copyright holders of
-that material) supplement the terms of this License with terms:
-
-    a) Disclaiming warranty or limiting liability differently from the
-    terms of sections 15 and 16 of this License; or
-
-    b) Requiring preservation of specified reasonable legal notices or
-    author attributions in that material or in the Appropriate Legal
-    Notices displayed by works containing it; or
-
-    c) Prohibiting misrepresentation of the origin of that material, or
-    requiring that modified versions of such material be marked in
-    reasonable ways as different from the original version; or
-
-    d) Limiting the use for publicity purposes of names of licensors or
-    authors of the material; or
-
-    e) Declining to grant rights under trademark law for use of some
-    trade names, trademarks, or service marks; or
-
-    f) Requiring indemnification of licensors and authors of that
-    material by anyone who conveys the material (or modified versions of
-    it) with contractual assumptions of liability to the recipient, for
-    any liability that these contractual assumptions directly impose on
-    those licensors and authors.
-
-  All other non-permissive additional terms are considered "further
-restrictions" within the meaning of section 10.  If the Program as you
-received it, or any part of it, contains a notice stating that it is
-governed by this License along with a term that is a further
-restriction, you may remove that term.  If a license document contains
-a further restriction but permits relicensing or conveying under this
-License, you may add to a covered work material governed by the terms
-of that license document, provided that the further restriction does
-not survive such relicensing or conveying.
-
-  If you add terms to a covered work in accord with this section, you
-must place, in the relevant source files, a statement of the
-additional terms that apply to those files, or a notice indicating
-where to find the applicable terms.
-
-  Additional terms, permissive or non-permissive, may be stated in the
-form of a separately written license, or stated as exceptions;
-the above requirements apply either way.
-
-  8. Termination.
-
-  You may not propagate or modify a covered work except as expressly
-provided under this License.  Any attempt otherwise to propagate or
-modify it is void, and will automatically terminate your rights under
-this License (including any patent licenses granted under the third
-paragraph of section 11).
-
-  However, if you cease all violation of this License, then your
-license from a particular copyright holder is reinstated (a)
-provisionally, unless and until the copyright holder explicitly and
-finally terminates your license, and (b) permanently, if the copyright
-holder fails to notify you of the violation by some reasonable means
-prior to 60 days after the cessation.
-
-  Moreover, your license from a particular copyright holder is
-reinstated permanently if the copyright holder notifies you of the
-violation by some reasonable means, this is the first time you have
-received notice of violation of this License (for any work) from that
-copyright holder, and you cure the violation prior to 30 days after
-your receipt of the notice.
-
-  Termination of your rights under this section does not terminate the
-licenses of parties who have received copies or rights from you under
-this License.  If your rights have been terminated and not permanently
-reinstated, you do not qualify to receive new licenses for the same
-material under section 10.
-
-  9. Acceptance Not Required for Having Copies.
-
-  You are not required to accept this License in order to receive or
-run a copy of the Program.  Ancillary propagation of a covered work
-occurring solely as a consequence of using peer-to-peer transmission
-to receive a copy likewise does not require acceptance.  However,
-nothing other than this License grants you permission to propagate or
-modify any covered work.  These actions infringe copyright if you do
-not accept this License.  Therefore, by modifying or propagating a
-covered work, you indicate your acceptance of this License to do so.
-
-  10. Automatic Licensing of Downstream Recipients.
-
-  Each time you convey a covered work, the recipient automatically
-receives a license from the original licensors, to run, modify and
-propagate that work, subject to this License.  You are not responsible
-for enforcing compliance by third parties with this License.
-
-  An "entity transaction" is a transaction transferring control of an
-organization, or substantially all assets of one, or subdividing an
-organization, or merging organizations.  If propagation of a covered
-work results from an entity transaction, each party to that
-transaction who receives a copy of the work also receives whatever
-licenses to the work the party's predecessor in interest had or could
-give under the previous paragraph, plus a right to possession of the
-Corresponding Source of the work from the predecessor in interest, if
-the predecessor has it or can get it with reasonable efforts.
-
-  You may not impose any further restrictions on the exercise of the
-rights granted or affirmed under this License.  For example, you may
-not impose a license fee, royalty, or other charge for exercise of
-rights granted under this License, and you may not initiate litigation
-(including a cross-claim or counterclaim in a lawsuit) alleging that
-any patent claim is infringed by making, using, selling, offering for
-sale, or importing the Program or any portion of it.
-
-  11. Patents.
-
-  A "contributor" is a copyright holder who authorizes use under this
-License of the Program or a work on which the Program is based.  The
-work thus licensed is called the contributor's "contributor version".
-
-  A contributor's "essential patent claims" are all patent claims
-owned or controlled by the contributor, whether already acquired or
-hereafter acquired, that would be infringed by some manner, permitted
-by this License, of making, using, or selling its contributor version,
-but do not include claims that would be infringed only as a
-consequence of further modification of the contributor version.  For
-purposes of this definition, "control" includes the right to grant
-patent sublicenses in a manner consistent with the requirements of
-this License.
-
-  Each contributor grants you a non-exclusive, worldwide, royalty-free
-patent license under the contributor's essential patent claims, to
-make, use, sell, offer for sale, import and otherwise run, modify and
-propagate the contents of its contributor version.
-
-  In the following three paragraphs, a "patent license" is any express
-agreement or commitment, however denominated, not to enforce a patent
-(such as an express permission to practice a patent or covenant not to
-sue for patent infringement).  To "grant" such a patent license to a
-party means to make such an agreement or commitment not to enforce a
-patent against the party.
-
-  If you convey a covered work, knowingly relying on a patent license,
-and the Corresponding Source of the work is not available for anyone
-to copy, free of charge and under the terms of this License, through a
-publicly available network server or other readily accessible means,
-then you must either (1) cause the Corresponding Source to be so
-available, or (2) arrange to deprive yourself of the benefit of the
-patent license for this particular work, or (3) arrange, in a manner
-consistent with the requirements of this License, to extend the patent
-license to downstream recipients.  "Knowingly relying" means you have
-actual knowledge that, but for the patent license, your conveying the
-covered work in a country, or your recipient's use of the covered work
-in a country, would infringe one or more identifiable patents in that
-country that you have reason to believe are valid.
-
-  If, pursuant to or in connection with a single transaction or
-arrangement, you convey, or propagate by procuring conveyance of, a
-covered work, and grant a patent license to some of the parties
-receiving the covered work authorizing them to use, propagate, modify
-or convey a specific copy of the covered work, then the patent license
-you grant is automatically extended to all recipients of the covered
-work and works based on it.
-
-  A patent license is "discriminatory" if it does not include within
-the scope of its coverage, prohibits the exercise of, or is
-conditioned on the non-exercise of one or more of the rights that are
-specifically granted under this License.  You may not convey a covered
-work if you are a party to an arrangement with a third party that is
-in the business of distributing software, under which you make payment
-to the third party based on the extent of your activity of conveying
-the work, and under which the third party grants, to any of the
-parties who would receive the covered work from you, a discriminatory
-patent license (a) in connection with copies of the covered work
-conveyed by you (or copies made from those copies), or (b) primarily
-for and in connection with specific products or compilations that
-contain the covered work, unless you entered into that arrangement,
-or that patent license was granted, prior to 28 March 2007.
-
-  Nothing in this License shall be construed as excluding or limiting
-any implied license or other defenses to infringement that may
-otherwise be available to you under applicable patent law.
-
-  12. No Surrender of Others' Freedom.
-
-  If conditions are imposed on you (whether by court order, agreement or
-otherwise) that contradict the conditions of this License, they do not
-excuse you from the conditions of this License.  If you cannot convey a
-covered work so as to satisfy simultaneously your obligations under this
-License and any other pertinent obligations, then as a consequence you may
-not convey it at all.  For example, if you agree to terms that obligate you
-to collect a royalty for further conveying from those to whom you convey
-the Program, the only way you could satisfy both those terms and this
-License would be to refrain entirely from conveying the Program.
-
-  13. Use with the GNU Affero General Public License.
-
-  Notwithstanding any other provision of this License, you have
-permission to link or combine any covered work with a work licensed
-under version 3 of the GNU Affero General Public License into a single
-combined work, and to convey the resulting work.  The terms of this
-License will continue to apply to the part which is the covered work,
-but the special requirements of the GNU Affero General Public License,
-section 13, concerning interaction through a network will apply to the
-combination as such.
-
-  14. Revised Versions of this License.
-
-  The Free Software Foundation may publish revised and/or new versions of
-the GNU General Public License from time to time.  Such new versions will
-be similar in spirit to the present version, but may differ in detail to
-address new problems or concerns.
-
-  Each version is given a distinguishing version number.  If the
-Program specifies that a certain numbered version of the GNU General
-Public License "or any later version" applies to it, you have the
-option of following the terms and conditions either of that numbered
-version or of any later version published by the Free Software
-Foundation.  If the Program does not specify a version number of the
-GNU General Public License, you may choose any version ever published
-by the Free Software Foundation.
-
-  If the Program specifies that a proxy can decide which future
-versions of the GNU General Public License can be used, that proxy's
-public statement of acceptance of a version permanently authorizes you
-to choose that version for the Program.
-
-  Later license versions may give you additional or different
-permissions.  However, no additional obligations are imposed on any
-author or copyright holder as a result of your choosing to follow a
-later version.
-
-  15. Disclaimer of Warranty.
-
-  THERE IS NO WARRANTY FOR THE PROGRAM, TO THE EXTENT PERMITTED BY
-APPLICABLE LAW.  EXCEPT WHEN OTHERWISE STATED IN WRITING THE COPYRIGHT
-HOLDERS AND/OR OTHER PARTIES PROVIDE THE PROGRAM "AS IS" WITHOUT WARRANTY
-OF ANY KIND, EITHER EXPRESSED OR IMPLIED, INCLUDING, BUT NOT LIMITED TO,
-THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
-PURPOSE.  THE ENTIRE RISK AS TO THE QUALITY AND PERFORMANCE OF THE PROGRAM
-IS WITH YOU.  SHOULD THE PROGRAM PROVE DEFECTIVE, YOU ASSUME THE COST OF
-ALL NECESSARY SERVICING, REPAIR OR CORRECTION.
-
-  16. Limitation of Liability.
-
-  IN NO EVENT UNLESS REQUIRED BY APPLICABLE LAW OR AGREED TO IN WRITING
-WILL ANY COPYRIGHT HOLDER, OR ANY OTHER PARTY WHO MODIFIES AND/OR CONVEYS
-THE PROGRAM AS PERMITTED ABOVE, BE LIABLE TO YOU FOR DAMAGES, INCLUDING ANY
-GENERAL, SPECIAL, INCIDENTAL OR CONSEQUENTIAL DAMAGES ARISING OUT OF THE
-USE OR INABILITY TO USE THE PROGRAM (INCLUDING BUT NOT LIMITED TO LOSS OF
-DATA OR DATA BEING RENDERED INACCURATE OR LOSSES SUSTAINED BY YOU OR THIRD
-PARTIES OR A FAILURE OF THE PROGRAM TO OPERATE WITH ANY OTHER PROGRAMS),
-EVEN IF SUCH HOLDER OR OTHER PARTY HAS BEEN ADVISED OF THE POSSIBILITY OF
-SUCH DAMAGES.
-
-  17. Interpretation of Sections 15 and 16.
-
-  If the disclaimer of warranty and limitation of liability provided
-above cannot be given local legal effect according to their terms,
-reviewing courts shall apply local law that most closely approximates
-an absolute waiver of all civil liability in connection with the
-Program, unless a warranty or assumption of liability accompanies a
-copy of the Program in return for a fee.
-
-                     END OF TERMS AND CONDITIONS
-
-            How to Apply These Terms to Your New Programs
-
-  If you develop a new program, and you want it to be of the greatest
-possible use to the public, the best way to achieve this is to make it
-free software which everyone can redistribute and change under these terms.
-
-  To do so, attach the following notices to the program.  It is safest
-to attach them to the start of each source file to most effectively
-state the exclusion of warranty; and each file should have at least
-the "copyright" line and a pointer to where the full notice is found.
-
-    <one line to give the program's name and a brief idea of what it does.>
-    Copyright (C) <year>  <name of author>
-
-    This program is free software: you can redistribute it and/or modify
-    it under the terms of the GNU General Public License as published by
-    the Free Software Foundation, either version 3 of the License, or
-    (at your option) any later version.
-
-    This program is distributed in the hope that it will be useful,
-    but WITHOUT ANY WARRANTY; without even the implied warranty of
-    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
-    GNU General Public License for more details.
-
-    You should have received a copy of the GNU General Public License
-    along with this program.  If not, see <https://www.gnu.org/licenses/>.
-
-Also add information on how to contact you by electronic and paper mail.
-
-  If the program does terminal interaction, make it output a short
-notice like this when it starts in an interactive mode:
-
-    <program>  Copyright (C) <year>  <name of author>
-    This program comes with ABSOLUTELY NO WARRANTY; for details type `show w'.
-    This is free software, and you are welcome to redistribute it
-    under certain conditions; type `show c' for details.
-
-The hypothetical commands `show w' and `show c' should show the appropriate
-parts of the General Public License.  Of course, your program's commands
-might be different; for a GUI interface, you would use an "about box".
-
-  You should also get your employer (if you work as a programmer) or school,
-if any, to sign a "copyright disclaimer" for the program, if necessary.
-For more information on this, and how to apply and follow the GNU GPL, see
-<https://www.gnu.org/licenses/>.
-
-  The GNU General Public License does not permit incorporating your program
-into proprietary programs.  If your program is a subroutine library, you
-may consider it more useful to permit linking proprietary applications with
-the library.  If this is what you want to do, use the GNU Lesser General
-Public License instead of this License.  But first, please read
-<https://www.gnu.org/licenses/why-not-lgpl.html>.
diff --git a/NAMESPACE b/NAMESPACE
new file mode 100644
index 0000000..9e56258
--- /dev/null
+++ b/NAMESPACE
@@ -0,0 +1,13 @@
+# Generated by roxygen2: do not edit by hand
+
+S3method(plot,gcm)
+S3method(plot,pcm)
+export(gcm)
+export(pcm)
+importFrom(glmnet,cv.glmnet)
+importFrom(ranger,ranger)
+importFrom(stats,model.frame)
+importFrom(stats,model.response)
+importFrom(stats,pnorm)
+importFrom(stats,predict)
+importFrom(stats,terms)
diff --git a/R/gcm.R b/R/gcm.R
new file mode 100644
index 0000000..eeced97
--- /dev/null
+++ b/R/gcm.R
@@ -0,0 +1,166 @@
+#' GCM using random forests
+#'
+#' @param Y Response
+#' @param X Covariates
+#' @param Z Covariates
+#' @param alternative Alternative
+#' @param ... Additional arguments to ranger
+#'
+#' @return Object of class htest
+#' @export
+#'
+#' @examples
+#' X <- matrix(rnorm(2e3), ncol = 2)
+#' colnames(X) <- c("X1", "X2")
+#' Z <- matrix(rnorm(2e3), ncol = 2)
+#' colnames(Z) <- c("Z1", "Z2")
+#' Y <- rnorm(1e3) # X[, 2] + Z[, 2] + rnorm(1e3)
+#' (gcm1 <- gcm(Y, X, Z))
+#' plot(gcm1)
+#'
+gcm <- function(Y, X, Z, alternative = c("two.sided", "less", "greater"), ...) {
+  alternative <- match.arg(alternative)
+  args <- if (length(list(...)) > 0) list(...) else NULL
+  YZ <- do.call("pcm_ranger", c(list(y = Y, x = Z), args))
+  XZ <- apply(as.data.frame(X), 2, \(tX) {
+    do.call("pcm_ranger", c(list(y = tX, x = Z), args))
+  })
+  rY <- Y - predict(YZ, data = Z)
+  preds <- lapply(XZ, predict.pcm_ranger, data = Z)
+  rX <- X - do.call("cbind", preds)
+  stat <- .gcm(rY, rX)
+  pval <- .compute_normal_pval(stat, alternative)
+
+  structure(list(
+    statistic = c("Z" = stat), p.value = pval,
+    hypothesis = c("E[cov(Y, X | Z)]" = "0"),
+    null.value = c("E[cov(Y, X | Z)]" = "0"), alternative = alternative,
+    method = paste0("Generalized covariance measure test"),
+    data.name = deparse(match.call(), width.cutoff = 80),
+    rY = rY, rX = rX), class = c("gcm", "htest"))
+
+}
+
+# Helpers -----------------------------------------------------------------
+
+.gcm <- function (r, e) {
+  dR <- NCOL(r)
+  dE <- NCOL(e)
+  nn <- NROW(r)
+  if (dR > 1 || dE > 1) {
+    R_mat <- matrix(r, nrow = nn, ncol = dE) * e
+    sigma <- crossprod(R_mat)/nn - tcrossprod(colMeans(R_mat))
+    eig <- eigen(sigma)
+    if (min(eig$values) < .Machine$double.eps)
+      warning("`vcov` of test statistic is not invertible")
+    siginvhalf <- eig$vectors %*% diag(eig$values^(-1/2)) %*%
+      t(eig$vectors)
+    tstat <- siginvhalf %*% colSums(R_mat)/sqrt(nn)
+    stat <- structure(sum(tstat^2), df = dR * dE)
+  }
+  else {
+    R <- r * e
+    R.sq <- R^2
+    meanR <- mean(R)
+    stat <- sqrt(nn) * meanR/sqrt(mean(R.sq) - meanR^2)
+  }
+  stat
+}
+
+#' @importFrom stats pnorm
+.compute_normal_pval <- function(stat, alternative) {
+  if (!is.null(df <- attr(stat, "df")))
+    return(stats::pchisq(stat, df = df, lower.tail = FALSE))
+  switch(
+    alternative,
+    "two.sided" = 2 * stats::pnorm(-abs(stat)),
+    "greater" = stats::pnorm(-abs(stat)),
+    "less" = stats::pnorm(abs(stat))
+  )
+}
+
+.rm_int <- function(x) {
+  if (all(x[, 1] == 1))
+    return(x[, -1L, drop = FALSE])
+  x
+}
+
+#' @importFrom stats terms
+.get_terms <- function(formula) {
+  if (is.null(formula))
+    return(NULL)
+  atms <- stats::terms(formula)
+  tms <- attr(atms, "term.labels")
+  resp <- all.vars(formula)[1]
+  ridx <- grep("|", tms, fixed = TRUE)
+  tms[ridx] <- paste0("(", tms[ridx], ")")
+  ie <- grep(":", tms, value = TRUE)
+  me <- grep(":", tms, value = TRUE, invert = TRUE)
+  list(all = tms, me = me, ie = ie, response = resp, terms = atms,
+       fml = formula)
+}
+
+# Ranger ------------------------------------------------------------------
+
+#' @importFrom stats model.response model.frame
+.ranger <- function(formula, data, ...) {
+  response <- stats::model.response(stats::model.frame(formula, data))
+  is_factor <- is.factor(response)
+  tms <- .get_terms(formula)
+  resp <- if (is_factor)
+    .rm_int(stats::model.matrix(~ response, contrasts.arg = list(
+      "response" = "contr.treatment")))
+  else response
+  tmp <- list(data = data, response = resp, is_factor = is_factor)
+  if (identical(tms$me, character(0))) {
+    if (is_factor)
+      return(structure(c(list(mean = base::colMeans(resp)), tmp),
+                       class = "ranger"))
+    else return(structure(c(list(mean = mean(as.numeric(response))),
+                            tmp), class = "ranger"))
+  }
+  ret <- ranger::ranger(formula, data, probability = is_factor, ...)
+  structure(c(ret, tmp), class = "ranger")
+}
+
+#' @importFrom stats predict
+residuals.ranger <- function(object, newdata = NULL, newy = NULL, ...) {
+  if (is.null(newdata))
+    newdata <- object$data
+  if (!is.null(newy))
+    newy <- if (object$is_factor)
+      .rm_int(stats::model.matrix(~ newy, contrasts.arg = list(
+        "newy" = "contr.treatment")))
+  else newy
+  if (is.null(newy))
+    newy <- object$response
+  if (!is.null(object$mean))
+    return(newy - object$mean)
+  preds <- stats::predict(object, data = newdata)$predictions
+  if (object$is_factor)
+    preds <- preds[, -1]
+  unname(newy - preds)
+}
+
+# Diagnostics -------------------------------------------------------------
+
+#' @exportS3Method plot gcm
+plot.gcm <- function(x, ...) {
+  .data <- NULL
+  pd <- tidyr::pivot_longer(data.frame(rY = x$rY, rX = unname(x$rX)),
+                            dplyr::starts_with("rX"))
+  if (requireNamespace("ggplot2")) {
+    p1 <- ggplot2::ggplot(pd, ggplot2::aes(x = .data[["value"]] ,
+                                           y = .data[["rY"]],
+                                           color = .data[["name"]])) +
+      ggplot2::geom_point(alpha = 0.3, show.legend = FALSE) +
+      ggplot2::geom_smooth(method = "lm", se = FALSE, show.legend = FALSE) +
+      ggplot2::theme_bw() +
+      ggplot2::labs(x = "Residuals X | Z", y = "Residuals Y | Z")
+    print(p1)
+  }
+  return(invisible(p1))
+}
+
+.mm <- function(preds, data)
+  .rm_int(stats::model.matrix(stats::reformulate(preds), data = data))
diff --git a/R/pcm.R b/R/pcm.R
new file mode 100644
index 0000000..5422e01
--- /dev/null
+++ b/R/pcm.R
@@ -0,0 +1,213 @@
+#' Conditional mean independence test
+#'
+#' @param Y Numeric; response
+#' @param X Numeric; covariates
+#' @param Z Numceric; covariates
+#' @param rep Number of repetitions
+#' @param est_vhat Estimate variance functional
+#' @param ... Additional arguments passed to \code{reg}
+#' @param reg Character; regression method
+#' @param mtry Argument passed to \code{ranger}
+#' @param ghat_args Arguments passed to reg
+#' @param do.check Save check data
+#'
+#' @importFrom ranger ranger
+#'
+#' @return Object of class '\code{htest}'
+#' @export
+#'
+#' @examples
+#' X <- matrix(rnorm(2e3), ncol = 2)
+#' colnames(X) <- c("X1", "X2")
+#' Z <- matrix(rnorm(2e3), ncol = 2)
+#' colnames(Z) <- c("Z1", "Z2")
+#' Y <- rnorm(1e3) # X[, 2] + Z[, 2] + rnorm(1e3)
+#' (pcm1 <- pcm(Y, X, Z, ghat_args = list(mtry = NULL, max.depth = NULL),
+#'     est_vhat = TRUE, do.check = TRUE))
+#' plot(pcm1)
+#'
+pcm <- function(Y, X, Z, rep = 1, est_vhat = TRUE,
+                reg = c("pcm_ranger", "pcm_lasso"),
+                ghat_args = NULL, mtry = identity,
+                do.check = FALSE, ...) {
+  reg <- match.arg(reg)
+  if (rep != 1) {
+    pcms <- lapply(seq_len(rep), \(iter) {
+      pcm(Y = Y, X = X, Z = Z, rep = 1, est_vhat = est_vhat, reg = reg,
+          ghat_args = ghat_args, mtry = mtry, do.check = FALSE, ... = ...)
+    })
+    stat <- mean(unlist(lapply(pcms, \(tst) tst$statistic)))
+    pval <- pnorm(stat, lower.tail = FALSE)
+    return(structure(list(
+      statistic = c("Z" = stat), p.value = pval,
+      hypothesis = c("E[Y | X, Z]" = "E[Y | Z]"),
+      null.value = c("E[Y | X, Z]" = "E[Y | Z]"), alternative = "two.sided",
+      method = paste0("Projected covariance measure test (K = ", rep, " repetitions)"),
+      all_tests = pcms, data.name = deparse(match.call(), width.cutoff = 80)),
+      class = c("pcm", "htest")))
+  }
+  ### Sample splitting
+  idx <- sample.int(NROW(Y), ceiling(NROW(Y) / 2))
+  ### Split 1
+  Ytr <- Y[idx]
+  Xtr <- data.frame(X)[idx, , drop = FALSE]
+  Ztr <- data.frame(Z)[idx, , drop = FALSE]
+  ### Split 2
+  Yte <- Y[-idx]
+  Xte <- data.frame(X)[-idx, , drop = FALSE]
+  Zte <- data.frame(Z)[-idx, , drop = FALSE]
+
+  ### Obtain hat{h}
+  ghat <- do.call(reg, c(list(y = Ytr, x = cbind(Xtr, Ztr)), ghat_args))
+  mtilde <- ranger(x = Ztr, y = pghat <- predict(ghat, data = cbind(Xtr, Ztr)),
+                   mtry = mtry, ...)
+  htilde <- \(X, Z) {
+    predict(ghat, data = cbind(X, Z)) -
+      predict(mtilde, data = Z)$predictions
+  }
+  rho <- mean((Ytr - mtilde$predictions) * predict(ghat, data = cbind(Xtr, Ztr)))
+  hhat <- \(X, Z) sign(rho) * htilde(X, Z)
+
+  ### Obtain hat{v}
+  if (est_vhat) {
+    vtilde <- ranger(x = cbind(Xtr, Ztr), y = (sqr <- (Ytr - predict(
+      ghat, data = cbind(Xtr, Ztr)))^2), mtry = mtry, ...)
+    a <- function(c) mean(sqr / (pmax(vtilde$predictions, 0) + c))
+    chat <- if (a(0) < 1) 0 else stats::uniroot(\(c) a(c) - 1, c(0, 10), extendInt = "yes")$root
+    vhat <- \(X, Z) pmax(predict(vtilde, data = cbind(X, Z))$predictions, 0) + chat
+  }
+  else
+    vhat <- \(X, Z) 1
+
+  ### Obtain residuals for test
+  fhat <- \(X, Z) hhat(X, Z) / vhat(X, Z)
+  mhatfhat <- ranger(x = Zte, y = (fhats <- fhat(Xte, Zte)), mtry = mtry, ...)
+  mhat <- do.call(reg, c(list(y = Yte, x = Zte), ghat_args))
+
+  ### Test
+  L <- (Yte - predict(mhat, data = Zte)) * (fhats - mhatfhat$predictions)
+  stat <- sqrt(length(idx)) * mean(L) / sqrt(mean(L^2) - mean(L)^2)
+  if (is.nan(stat)) stat <- -Inf
+  pval <- pnorm(stat, lower.tail = FALSE)
+
+  dcheck <- if (do.check) {
+    pmtilde <- predict(mtilde, data = Ztr)$predictions
+    # pvtilde <- predict(vtilde, data = cbind(Xtr, Ztr))$predictions
+    pghat_test <- predict(ghat, data = cbind(Xte, Zte))
+    mhats <- predict(mhat, data = Zte)
+    list(
+      train = data.frame(
+        id = idx,
+        fitted_ghat = pghat,
+        resid_ghat = Ytr - pghat,
+        fitted_mtilde = pmtilde,
+        resid_mtilde = pghat - pmtilde,
+        # fitted_vtilde = pvtilde,
+        # resid_vtilde = sqr - pvtilde,
+        hhat = hhat(Xtr, Ztr),
+        vhat = vhat(Xtr, Ztr),
+        MSE_YZ = (Ytr - predict(mhat, data = Ztr))^2,
+        MSE_YXZ = (Ytr - pghat)^2
+      ),
+      test = data.frame(
+        id = setdiff(seq_len(NROW(Y)), idx),
+        fitted_ghat = pghat_test,
+        resid_ghat = Yte - pghat_test,
+        hhat = hhat(Xte, Zte),
+        vhat = vhat(Xte, Zte),
+        fhat = fhats,
+        resid_mhat = (Yte - mhats),
+        mhat = mhats,
+        resid_mhatfhat = (fhats - mhatfhat$predictions),
+        mhatfhat = mhats,
+        MSE_YXZ = (Yte - pghat_test)^2,
+        MSE_YZ = (Yte - mhats)^2
+      )
+    )
+  } else NULL
+
+  structure(list(
+    statistic = c("Z" = stat), p.value = pval,
+    hypothesis = c("E[Y | X, Z]" = "E[Y | Z]"),
+    null.value = c("E[Y | X, Z]" = "E[Y | Z]"), alternative = "two.sided",
+    method = paste0("Projected covariance measure test"),
+    data.name = deparse(match.call(), width.cutoff = 80),
+    check.data = dcheck), class = c("pcm", "htest"))
+}
+
+# Regressions -------------------------------------------------------------
+
+pcm_ranger <- function(y, x, ...) {
+  args <- list(...)
+  if (length(unique(y)) == 2) {
+    y <- as.factor(y)
+    args$probability <- TRUE
+  }
+  rf <- do.call("ranger", c(list(y = y, x = x), args))
+  class(rf) <- c("pcm_ranger", class(rf))
+  rf
+}
+
+predict.pcm_ranger <- function(object, data = NULL, ...) {
+  class(object) <- class(object)[-1]
+  preds <- predict(object, data = data)$predictions
+  if (object$treetype == "Probability estimation")
+    preds <- preds[, 2]
+  preds
+}
+
+#' @importFrom glmnet cv.glmnet
+pcm_lasso <- function(y, x, ...) {
+  obj <- cv.glmnet(y = y, x = as.matrix(x), ...)
+  class(obj) <- c("pcm_lasso", class(obj))
+  obj
+}
+
+predict.pcm_lasso <- function(object, data = NULL, ...) {
+  class(object) <- class(object)[-1]
+  predict(object, newx = as.matrix(data), s = object$lambda.min)[, 1]
+}
+
+# Diagnostics -------------------------------------------------------------
+
+#' @exportS3Method plot pcm
+plot.pcm <- function(x, ...) {
+  .data <- NULL
+  train <- x$check.data$train
+  test <- x$check.data$test
+  mse_train <- data.frame(MSE_YZ = train$MSE_YZ, MSE_YXZ = train$MSE_YXZ, set = "train")
+  mse_test <- data.frame(MSE_YZ = test$MSE_YZ, MSE_YXZ = test$MSE_YXZ, set = "test")
+  mses <- tidyr::pivot_longer(dplyr::bind_rows(mse_train, mse_test),
+                              dplyr::starts_with("MSE"))
+  if (is.null(train))
+    return("Nothing to plot. Consider running `pcm` with `do.check = TRUE`")
+  if (requireNamespace("ggplot2") && requireNamespace("tidyr") && requireNamespace("ggpubr")) {
+    mpl <- \(xx, yy, pdat, ...) {
+      ggplot2::ggplot(pdat, ggplot2::aes(y = .data[[yy]], x = .data[[xx]])) +
+        ggplot2::geom_point(alpha = 0.3) +
+        ggplot2::geom_smooth(se = FALSE, method = "lm") +
+        ggplot2::theme_bw()
+    }
+    # ptrain <- lapply(c("ghat", "mtilde", "vtilde"), \(tx) {
+    #   mpl(paste0("fitted_", tx), paste0("resid_", tx), train, "direction")
+    # })
+    # ptest <- apply(data.frame(y = c("resid_mhatfhat", "resid_mhat"),
+    #                           x = c("mhatfhat", "resid_mhatfhat")), 1,
+    #                \(x) mpl(x[2], x[1], test), simplify = FALSE)
+    # p1 <- mpl("mhatfhat", "resid_mhatfhat", test) +
+    #   ggplot2::labs(x = "f(X,Z)", y = "Residuals f(X,Z) | Z")
+    p2 <- mpl("resid_mhat", "resid_mhatfhat", test) +
+      ggplot2::labs(x = "Residuals f(X, Z) | Z", y = "Residuals Y | Z")
+    p3 <- ggplot2::ggplot(mses, ggplot2::aes(y = .data[["value"]],
+                                             x = .data[["name"]],
+                                             color = .data[["set"]])) +
+      ggplot2::geom_violin(width = 0.5, position = ggplot2::position_dodge(width = 0.7)) +
+      ggplot2::geom_boxplot(width = 0.3, outlier.shape = NA, position = ggplot2::position_dodge(width = 0.7)) +
+      ggplot2::geom_jitter(alpha = 0.1, position = ggplot2::position_dodge(width = 0.7)) +
+      ggplot2::theme_bw() +
+      ggplot2::scale_y_log10() +
+      ggplot2::labs(y = "MSE contributions", x = "Regression")
+    print(p2)
+  }
+  return(invisible(p2))
+}
diff --git a/README.md b/README.md
new file mode 100644
index 0000000..e486ad7
--- /dev/null
+++ b/README.md
@@ -0,0 +1,9 @@
+<!-- badges: start -->
+[![R-CMD-check](https://github.com/LucasKook/comet/actions/workflows/R-CMD-check.yaml/badge.svg)](https://github.com/LucasKook/comet/actions/workflows/R-CMD-check.yaml)
+<!-- badges: end -->
+
+# Covariance measure tests
+
+
+# Replication materials
+
diff --git a/inst/Makefile b/inst/Makefile
new file mode 100644
index 0000000..16fc0e4
--- /dev/null
+++ b/inst/Makefile
@@ -0,0 +1,23 @@
+dependencies:
+	(cd .. && Rscript --vanilla inst/dependencies.R > dependencies.out)
+
+mimic:
+	(cd .. && Rscript --vanilla inst/application/mimic.R > mimic.out)
+
+multiomics:
+	(cd .. && Rscript --vanilla inst/application/multiomics.R > multiomics.out)
+
+ccle:
+	(cd .. && Rscript --vanilla inst/application/ccle.R > ccle.out)
+
+dep:
+	bash dispatch.sh run-install.sh 1200
+
+mim:
+	bash dispatch.sh run-mimic.sh 21600
+
+mul:
+	bash dispatch.sh run-multiomics.sh 3600
+
+sig:
+  bash dispatch.sh run-ccle.sh 3600
diff --git a/inst/code/ccle.R b/inst/code/ccle.R
new file mode 100644
index 0000000..a54c8b8
--- /dev/null
+++ b/inst/code/ccle.R
@@ -0,0 +1,94 @@
+### Pre-process CCLE data
+### LK 2024
+
+bpath <- "inst/data/ccle"
+set.seed(2410)
+
+# DEPs --------------------------------------------------------------------
+
+library("tidyverse")
+library("survival")
+devtools::load_all()
+
+# Read --------------------------------------------------------------------
+
+cmp <- "PLX4720"
+nrep <- 10
+
+# expression <- read_table("inst/data/ccle/expression.txt")
+mutations <- read_table("inst/data/ccle/mutation.txt") |>
+  filter(Name %in% grep("_MUT$", Name, value = TRUE))
+response <- read_csv("inst/data/ccle/response.csv") |>
+  filter(Compound == cmp)
+
+idx <- match(response$`CCLE Cell Line Name`, colnames(mutations))
+
+resp <- response |> filter(!is.na(idx))
+mut <- mutations[, c(1, idx[!is.na(idx)])] |>
+  column_to_rownames("Name")
+
+cors <- apply(mut, 1, \(muts) {
+  tc <- cor(resp$Amax, muts)
+  ifelse(!is.na(tc), tc, 0)
+})
+
+screen <- cors[abs(cors) > 0.05 & !is.na(cors)]
+dat <- data.frame(Y = resp$Amax, t(mut[names(screen),]))
+
+# Run ---------------------------------------------------------------------
+
+response <- "Y"
+tgt <- paste0(c("BRAF.V600E", "BRAF.MC", "HIP1", "FLT3", "CDC42BPA",
+                "THBS3", "DNMT1", "PRKD1", "PIP5K1A", "MAP3K5"), "_MUT")
+muts <- setdiff(colnames(dat), c(response, tgt))
+
+tmtry <- NULL
+tmd <- NULL
+
+pb <- txtProgressBar(min = 0, max = length(tgt), style = 3)
+out <- lapply(seq_along(tgt), \(tmut) {
+  setTxtProgressBar(pb, tmut)
+  pcm <- pcm(dat[, response], dat[, tgt[tmut]], dat[, c(tgt[-tmut], muts)],
+             mtry = tmtry, rep = nrep, est_vhat = TRUE, ghat_args = list(
+               mtry = tmtry, max.depth = tmd))
+  gcm <- ranger_gcm(dat[, response], dat[, tgt[tmut]], dat[, c(tgt[-tmut], muts)],
+                    mtry = identity, max.depth = tmd)
+  naive <- cor.test(x = dat[, tgt[tmut]], y = dat[, response])
+  data.frame(gene = tgt[tmut],
+             pcm = pcm$p.value,
+             gcm = gcm$p.value,
+             naive = naive$p.value)
+})
+
+### Tabulate results
+(res <- bind_rows(out))
+knitr::kable(t(res[, c("gcm", "pcm")]), col.names = res$gene,
+             format = "latex", booktabs = TRUE, digits = 3)
+
+# # Subsampling -------------------------------------------------------------
+#
+# pb <- txtProgressBar(min = 0, max = length(tgt), style = 3)
+# out <- lapply(seq_along(tgt), \(tmut) {
+#   setTxtProgressBar(pb, tmut)
+#   n <- NROW(dat)
+#   folds <- sample(rep(1:2, ceiling(n/2)), n)
+#   lapply(1:2, \(iter) {
+#     tidx <- which(folds == iter)
+#     pcm <- pcm(dat[tidx, response], dat[tidx, tgt[tmut]], dat[tidx, c(tgt[-tmut], muts)],
+#                mtry = tmtry, rep = nrep, est_vhat = TRUE, ghat_args = list(
+#                  mtry = tmtry, max.depth = tmd))
+#     gcm <- ranger_gcm(dat[tidx, response], dat[tidx, tgt[tmut]], dat[tidx, c(tgt[-tmut], muts)],
+#                       mtry = identity, max.depth = tmd)
+#     naive <- cor.test(x = dat[tidx, tgt[tmut]], y = dat[tidx, response])
+#     data.frame(gene = tgt[tmut],
+#                pcm = pcm$p.value,
+#                gcm = gcm$p.value,
+#                naive = naive$p.value,
+#                fold = iter)
+#   }) |> bind_rows()
+# }) |> bind_rows()
+#
+# out |> pivot_wider(names_from = fold, values_from = pcm:naive) |>
+#   ggplot(aes(x = pcm_1, y = pcm_2)) +
+#   geom_point() +
+#   geom_abline(intercept = 0, slope = 1)
diff --git a/inst/code/mimic.R b/inst/code/mimic.R
new file mode 100644
index 0000000..3a6a285
--- /dev/null
+++ b/inst/code/mimic.R
@@ -0,0 +1,163 @@
+### MIMIC analysis GCM/PCM
+### LK 2023
+
+set.seed(2410)
+nrep <- 5
+run_full <- TRUE
+red <- 111 # 98% variance
+
+# Dependencies ------------------------------------------------------------
+
+library("readr")
+library("ranger")
+library("ggplot2")
+devtools::load_all()
+
+# Read data ---------------------------------------------------------------
+
+## Read meta data
+bpath <- "inst/data/mimic/MIMIC"
+
+if (file.exists(temb <- file.path(bpath, "emb.csv"))) {
+  emb <- as.data.frame(read_csv(temb))
+  emb$race <- as.factor(emb$race)
+  emb$sex <- as.factor(emb$sex)
+} else {
+  library("reticulate")
+  np <-import("numpy")
+  meta <- read_csv(file.path(bpath, "mimic_cfm_train_meta.csv"))
+  meta_test <- read_csv(file.path(bpath, "mimic_cfm_test_meta.csv"))
+  meta$sex <- as.factor(meta$sex)
+  meta_test$sex <- as.factor(meta_test$sex)
+  ## race: Merge Asian, White, Black
+  meta$race <- as.factor(meta$race)
+  meta_test$race <- as.factor(meta_test$race)
+
+  if(file.exists(file.path(bpath, "emb_svd.RDS"))) {
+
+    emb <- np$load(file.path(bpath, "mimic_cfm_train_emb.npy"))
+    emb_svd <- readRDS(file.path(bpath, "emb_svd.RDS"))
+
+  } else {
+
+    emb <- np$load(file.path(bpath, "mimic_cfm_train_emb.npy"))
+
+    ## SVD
+    emb_svd <- fast.svd(emb)
+    saveRDS(emb_svd, file = file.path(bpath, "emb_svd.RDS"))
+    plot(cumsum(emb_svd$d^2/sum(emb_svd$d^2)), type="b")
+
+  }
+
+  # create reduced embedding from SVD
+  emb <- emb %*% emb_svd$v[, 1:red]
+  emb_test <- tcrossprod(
+    np$load(file.path(bpath, "mimic_cfm_test_emb.npy")), emb_svd$v[1:red,])
+
+  ## Combine
+  emb <- cbind(as.data.frame(emb), meta)
+  emb_test <- cbind(as.data.frame(emb_test), meta_test)
+  rm(meta, meta_test); gc()
+
+  ## Create response
+  name_resp <- "Pleural Effusion"
+  emb$resp <- emb[[name_resp]]
+
+  write_csv(emb, temb)
+}
+
+# Test --------------------------------------------------------------------
+
+### Variables
+nY <- "resp"
+nX <- "race"
+nC <- c("sex", "age")
+nZ <- paste0("V", 1:red)
+
+### RF params
+tmtry <- NULL # \(n) ceiling(sqrt(n))
+rargs <- list(mtry = NULL, max.depth = NULL)
+
+### Test function
+run_tests <- function(splits = 20, max_size = 1e4, verbose = FALSE) {
+  pb <- txtProgressBar(min = 0, max = splits, style = 3)
+  n <- NROW(emb)
+  folds <- sample(rep(1:splits, ceiling(n/splits)), n)
+  lapply(seq_len(splits), \(iter) {
+    setTxtProgressBar(pb, iter)
+    idx <- which(folds == iter)
+    if (length(idx) > max_size)
+      idx <- idx[1:max_size]
+    nemb <- emb[idx, ]
+
+    ### Test resp _||_ race | emb, age, sex ### GCM/PCM
+    gcm1 <- ranger_gcm(nemb[, nY], .mm(nX, nemb), .mm(c(nZ, nC), nemb),
+                       mtry = rargs$mtry, max.depth = rargs$max.depth)
+    if (verbose)
+      cat("\nGCM1 done")
+    pcm1 <- pcm(nemb[, nY], .mm(nX, nemb), .mm(c(nZ, nC), nemb), rep = nrep,
+                mtry = tmtry, est_vhat = TRUE, ghat_args = c(rargs, list(
+                  probability = TRUE)), do.check = FALSE)
+    if (verbose)
+      cat("\nPCM1 done")
+    ### Test resp _||_ emb | race, age, sex ### GCM/PCM
+    gcm2 <- ranger_gcm(nemb[, nY], .mm(nZ, nemb), .mm(c(nX, nC), nemb),
+                       mtry = rargs$mtry, max.depth = rargs$max.depth)
+    if (verbose)
+      cat("\nGCM2 done")
+    pcm2 <- pcm(nemb[, nY], .mm(nZ, nemb), .mm(c(nX, nC), nemb), rep = nrep,
+                mtry = tmtry, est_vhat = TRUE, ghat_args = c(rargs, list(
+                  probability = TRUE)), do.check = FALSE)
+    if (verbose)
+      cat("\nPCM2 done")
+    ### Return
+    list(GCM1 = gcm1, PCM1 = pcm1, GCM2 = gcm2, PCM2 = pcm2)
+  })
+}
+
+### RUN either on full data or with splits
+if (run_full) {
+
+  print(full <- run_tests(1, nrow(emb), verbose = TRUE))
+  pvals <- c(
+    "GCM1" = -log10(full[[1]]$GCM1$p.value),
+    "PCM1" = -log10(full[[1]]$PCM1$p.value),
+    "GCM2" = -pchisq(unname(full[[1]]$GCM2$statistic), log.p = TRUE,
+                     df = 111, lower.tail = FALSE) / log(10),
+    "PCM2" = -pnorm(unname(full[[1]]$PCM2$statistic), log.p = TRUE,
+                    lower.tail = FALSE) / log(10)
+  )
+  cat("Negative log10 p-values\n")
+  print(pvals)
+  saveRDS(full, "inst/results/mimic-full.rds")
+
+} else {
+
+  ms <- 150 * 4^(0:2)
+  nsplit <- 75
+  out <- dplyr::bind_rows(lapply(ms, \(nmax) {
+    res <- run_tests(splits = nsplit, max_size = nmax)
+    dplyr::bind_rows(lapply(res, \(x) unlist(lapply(x, \(y) y$p.value)))) |>
+      dplyr::mutate(n = nmax)
+  }))
+
+  out |> tidyr::pivot_longer(-n) |>
+    dplyr::mutate(hypothesis = dplyr::case_when(
+      name %in% c("GCM1", "PCM1") ~ "PE~'_||_'~race~'|'~x*'-'*ray*','*~sex*','*~age",
+      name %in% c("GCM2", "PCM2") ~ "PE~'_||_'~x*'-'*ray~'|'~race*','*~sex*','*~age"
+    )) |>
+    ggplot(aes(x = ordered(n), color = stringr::str_remove(name, "[0-9]"), y = -log10(value))) +
+    geom_violin(width = 0.5, position = position_dodge(width = 0.5)) +
+    geom_boxplot(width = 0.3, position = position_dodge(width = 0.5), outlier.shape = NA) +
+    ggbeeswarm::geom_quasirandom(width = 0.1, dodge.width = 0.5, alpha = 0.3, size = 0.5) +
+    facet_wrap(~ hypothesis, labeller = label_parsed, scales = "free") +
+    labs(x = "Sample size", y = parse(text = "-log[10]~p*'-'*value"), color = "COMET") +
+    theme_bw() +
+    geom_hline(yintercept = -log10(0.05), color = "darkred", linetype = 3) +
+    scale_color_brewer(palette = "Dark2") +
+    theme(legend.position = "top")
+
+  readr::write_csv(out, "inst/results/mimic-pvals.csv")
+  ggsave("inst/figures/mimic-pval.pdf", height = 3.5, width = 6, scale = 0.85)
+
+}
diff --git a/inst/code/multiomics.R b/inst/code/multiomics.R
new file mode 100644
index 0000000..94909e3
--- /dev/null
+++ b/inst/code/multiomics.R
@@ -0,0 +1,113 @@
+### Pre-process multiomics data
+### LK 2024
+
+bpath <- "inst/data/multiomics"
+set.seed(2410)
+
+# DEPs --------------------------------------------------------------------
+
+library("tidyverse")
+library("survival")
+devtools::load_all()
+
+# Read --------------------------------------------------------------------
+
+rna <- read_tsv(file.path(bpath, "rna.tsv"))
+meth <- read_tsv(file.path(bpath, "meth.tsv"))
+mir <- read_tsv(file.path(bpath, "mir.tsv"))
+surv <- read_tsv(file.path(bpath, "survival.tsv")) %>%
+  mutate(surv = Surv(days, event))
+
+rna_idx <- match(rna$Samples, surv$Samples)
+mir_idx <- match(mir$Samples, surv$Samples)
+meth_idx <- match(meth$Samples, surv$Samples)
+
+stopifnot(all(rna_idx == mir_idx))
+stopifnot(all(meth_idx == mir_idx))
+
+dat <- data.frame(
+  surv = surv$surv[rna_idx],
+  rna = scale(as.matrix(rna[, -1])),
+  mir = scale(as.matrix(mir[, -1])),
+  meth = scale(as.matrix(meth[, -1]))
+)
+
+rm(list = c("mir", "rna", "surv", "meth"))
+
+# Run ---------------------------------------------------------------------
+
+response <- "surv"
+names_mir <- grep("^mir", colnames(dat), value = TRUE)
+names_rna <- grep("^rna", colnames(dat), value = TRUE)
+names_meth <- grep("^meth", colnames(dat), value = TRUE)
+
+rargs <- list(mtry = NULL, max.depth = NULL, probability = TRUE)
+tmtry <- NULL
+nrep <- 10
+
+(t1 <- pcm(dat[, response][, 2], dat[, names_rna], dat[, c(names_mir, names_meth)],
+           mtry = tmtry, rep = nrep, est_vhat = TRUE, ghat_args = rargs))
+(t2 <- pcm(dat[, response][, 2], dat[, names_mir], dat[, c(names_rna, names_meth)],
+           mtry = tmtry, rep = nrep, est_vhat = TRUE, ghat_args = rargs))
+(t3 <- pcm(dat[, response][, 2], dat[, names_meth], dat[, c(names_rna, names_mir)],
+           mtry = tmtry, rep = nrep, est_vhat = TRUE, ghat_args = rargs))
+
+(t4 <- pcm(dat[, response][, 2], dat[, names_rna], dat[, c(names_mir, names_meth)],
+           mtry = tmtry, rep = nrep, est_vhat = TRUE, reg = "pcm_lasso"))
+(t5 <- pcm(dat[, response][, 2], dat[, names_mir], dat[, c(names_rna, names_meth)],
+           mtry = tmtry, rep = nrep, est_vhat = TRUE, reg = "pcm_lasso"))
+(t6 <- pcm(dat[, response][, 2], dat[, names_meth], dat[, c(names_rna, names_mir)],
+           mtry = tmtry, rep = nrep, est_vhat = TRUE, reg = "pcm_lasso"))
+
+### Don't run
+# (g1 <- ranger_gcm(dat[, response][, 2], dat[, names_rna], dat[, c(names_mir, names_meth)],
+#                   mtry = rargs$mtry, max.depth = rargs$max.depth))
+# (g2 <- ranger_gcm(dat[, response][, 2], dat[, names_mir], dat[, c(names_rna, names_meth)],
+#                   mtry = rargs$mtry, max.depth = rargs$max.depth))
+# (g3 <- ranger_gcm(dat[, response][, 2], dat[, names_meth], dat[, c(names_rna, names_mir)],
+#                   mtry = rargs$mtry, max.depth = rargs$max.depth))
+
+# screen_rna <- names(which(sapply(names_rna, \(tx) cor(dat[, response][, 2], dat[, tx])) > 0.2))
+# screen_mir <- names(which(sapply(names_mir, \(tx) cor(dat[, response][, 2], dat[, tx])) > 0.2))
+# screen_meth <- names(which(sapply(names_meth, \(tx) cor(dat[, response][, 2], dat[, tx])) > 0.2))
+#
+# (t1 <- pcm(dat[, response][, 2], dat[, screen_rna], dat[, c(screen_mir, screen_meth)],
+#            mtry = tmtry, rep = nrep, est_vhat = TRUE, ghat_args = rargs))
+# (t2 <- pcm(dat[, response][, 2], dat[, screen_mir], dat[, c(screen_rna, screen_meth)],
+#            mtry = tmtry, rep = nrep, est_vhat = TRUE, ghat_args = rargs))
+# (t3 <- pcm(dat[, response][, 2], dat[, screen_meth], dat[, c(screen_rna, screen_mir)],
+#            mtry = tmtry, rep = nrep, est_vhat = TRUE, ghat_args = rargs))
+# (g1 <- ranger_gcm(dat[, response][, 2], as.matrix(dat[, screen_rna]),
+#                   as.matrix(dat[, c(screen_mir, screen_meth)]),
+#                   mtry = rargs$mtry, max.depth = rargs$max.depth))
+# (g2 <- ranger_gcm(dat[, response][, 2], as.matrix(dat[, screen_mir]),
+#                   as.matrix(dat[, c(screen_rna, screen_meth)]),
+#                   mtry = rargs$mtry, max.depth = rargs$max.depth))
+# (g3 <- ranger_gcm(dat[, response][, 2], as.matrix(dat[, screen_meth]),
+#                   as.matrix(dat[, c(screen_rna, screen_mir)]),
+#                   mtry = rargs$mtry, max.depth = rargs$max.depth))
+
+# SurvGCM -----------------------------------------------------------------
+
+# dmeth <- dat[, c(response, screen_meth)]
+# pb <- txtProgressBar(0, length(screen_meth), style = 3)
+# out <- lapply(seq_along(screen_meth), \(tmeth) {
+#   setTxtProgressBar(pb, tmeth)
+#   fm <- reformulate(screen_meth[-tmeth], response)
+#   mY <- survforest(fm, dmeth)
+#   tst <- gcm(mY, reformulate(screen_meth[-tmeth], screen_meth[tmeth]), data = dmeth)
+#   mY2 <- .ranger(reformulate(screen_meth[-tmeth], "surv[, 2]"), dmeth)
+#   tst2 <- gcm(mY2, reformulate(screen_meth[-tmeth], screen_meth[tmeth]), data = dmeth)
+#   # mY2 <- coxph(fm, dmeth)
+#   # tst2 <- gcm(mY2, reformulate(screen_meth[-tmeth], screen_meth[tmeth]), data = dmeth)
+#   data.frame(meth = screen_meth[tmeth], survgcm = tst$p.value, gcm = tst2$p.value)
+# })
+# bind_rows(out)
+#
+# ggplot(bind_rows(out), aes(x = gcm, y = coxph, color = gcm < 0.05 & coxph < 0.05)) +
+#   geom_point(show.legend = FALSE) +
+#   theme_bw() +
+#   scale_x_log10() +
+#   scale_y_log10() +
+#   geom_abline(slope = 1, intercept = 0) +
+#   labs(x = "Survival GCM", y = "Binary GCM")
diff --git a/inst/dependencies.R b/inst/dependencies.R
new file mode 100644
index 0000000..ba05f81
--- /dev/null
+++ b/inst/dependencies.R
@@ -0,0 +1,15 @@
+# Install all dependencies for reproducing the results in the paper
+# LK 2024
+
+.libPaths(c("~/tutu/lib", .libPaths()))
+
+### Install remotes for installing all other packages
+install.packages("remotes", repos = "https://cloud.r-project.org")
+
+### Install all other packages from CRAN
+pkgs <- c("tidyverse", "ranger", "mlt", "sandwich", "glmnet", "tram",
+          "reticulate", "mgcv", "survival")
+remotes::install_cran(pkgs)
+
+### Install comet
+remotes::install_local(force = TRUE)
diff --git a/inst/figures/components-overview.R b/inst/figures/components-overview.R
new file mode 100644
index 0000000..54839f3
--- /dev/null
+++ b/inst/figures/components-overview.R
@@ -0,0 +1,103 @@
+### Figure 1
+### LK 2024
+
+set.seed(2410)
+
+# DEPs --------------------------------------------------------------------
+
+library("tidyverse")
+library("mgcv")
+devtools::load_all()
+theme_set(theme_bw() + theme(text = element_text(size = 13.5)))
+
+pd1 <- data.frame(id = seq_len(1e3), pvalue = -log10(p.adjust(runif(1e3)^2)))
+
+ggplot(pd1, aes(x = id, y = pvalue, col = pvalue)) +
+  geom_col(show.legend = FALSE) +
+  geom_hline(yintercept = -log10(0.05), color = "darkred", linetype = 2) +
+  labs(x = element_blank(), y = parse(text = "-log[10]~p*'-'*value"),
+       color = parse(text = "-log[10]~p*'-'*value")) +
+  scale_color_viridis_c(begin = 0.1, end = 0.9, option = "C") +
+  theme(legend.position = "top", axis.text.x = element_blank(),
+        axis.ticks.x = element_blank())
+
+ggsave("inst/figures/manhatten.pdf", height = 2, width = 3)
+
+# GCM v PCM ---------------------------------------------------------------
+
+f <- \(x) 1 + sin(3 * X^2)
+g <- \(z) 1 + z^3
+mu <- \(x, z) f(x) * g(z)
+
+X <- runif(n <- 3e2, -1, 1)
+Z <- runif(n, -1, 1)
+Y <- mu(X, Z) + rnorm(n, sd = 0.05)
+pd2 <- data.frame(Y = Y, X = X, Z = Z)
+
+mY <- .ranger(Y ~ Z, data = pd2)
+pd2$RY <- residuals.ranger(mY)
+mX <- .ranger(X ~ Z, data = pd2)
+pd2$RX <- residuals.ranger(mX)
+pd2$TX <- (f(X) - mean(f(X))) * g(Z)
+mTX <- .ranger(TX ~ Z, data = pd2)
+pd2$RTX <- residuals.ranger(mTX)
+
+ggplot(pd2 |> mutate(id = 1:n) |> pivot_longer(X:Z, names_to = "Predictor",
+                                               values_to = "value"),
+       aes(y = Y, x = value, color = Predictor, group = Predictor)) +
+  geom_point(alpha = 0.5) +
+  geom_smooth(method = "gam", se = FALSE) +
+  scale_color_brewer(palette = "Dark2") +
+  theme(legend.position = "bottom", legend.box.spacing = unit(-0.3, "cm")) +
+  labs(y = "Response Y", subtitle = "Y depends on X given Z",
+       x = element_blank())
+
+ggsave("inst/figures/pcmgcmdata.pdf", height = 3.3, width = 3.3)
+
+ggplot(pd2, aes(x = X, y = RY)) +
+  geom_point(alpha = 0.5) +
+  geom_smooth(method = "lm", se = FALSE, col = "gray80") +
+  geom_smooth(method = "gam", se = FALSE, col = "#1E88E5", linetype = 2, linewidth = 1.8) +
+  labs(x = "X", y = "Residual Y on Z",
+       subtitle = "GCM (fails to reject)") +
+  geom_label(aes(x = -0.5, y = 0.05,
+                 label = paste0("R==", round(cor(RX, RY), 2))),
+             inherit.aes = FALSE, parse = TRUE)
+
+ggsave("inst/figures/gcm.pdf", height = 3.3, width = 3.3)
+
+ggplot(pd2, aes(x = TX, y = RY)) +
+  geom_point(alpha = 0.5) +
+  geom_smooth(method = "lm", se = FALSE, col = "gray80") +
+  labs(x = "Transformed X", y = "Residual Y on Z",
+       subtitle = "PCM (correctly rejects)") +
+  geom_label(aes(x = 0.3, y = 0.1,
+                 label = paste0("R==", round(cor(TX, RY), 2))),
+             inherit.aes = FALSE, parse = TRUE)
+
+ggsave("inst/figures/pcm.pdf", height = 3.3, width = 3.3)
+
+# Modality selection ------------------------------------------------------
+
+pd3 <- data.frame(MSE1 = runif(3e1, min = 1, max = 3),
+                  MSE2 = runif(3e1, min = 0.3, max = 2.5))
+ggplot(pivot_longer(pd3, MSE1:MSE2), aes(x = name, y = value)) +
+  geom_violin(width = 0.5) +
+  geom_boxplot(width = 0.3) +
+  ggbeeswarm::geom_quasirandom(alpha = 0.3, width = 0.2) +
+  labs(x = element_blank(), y = "MSPE") +
+  scale_x_discrete(labels = c("MSE1" = "Y | X", "MSE2" = "Y | X, Z"))
+
+ggsave("inst/figures/modsec.pdf", height = 2, width = 3)
+
+# GCM ---------------------------------------------------------------------
+
+X <- matrix(runif(3e2), ncol = 1)
+colnames(X) <- c("X1")
+Z <- matrix(runif(6e2), ncol = 2)
+colnames(Z) <- c("Z1", "Z2")
+Y <- exp(X[, 1]) + Z[, 2]^3 + exp(Z[, 1]) + rnorm(3e2)
+(gcm1 <- ranger_gcm(Y, X, Z))
+plot(gcm1) + labs(x = "Transformation of X and Z")
+
+ggsave("inst/figures/simple-gcm.pdf", height = 2, width = 2.5)
diff --git a/man/gcm.Rd b/man/gcm.Rd
new file mode 100644
index 0000000..8037965
--- /dev/null
+++ b/man/gcm.Rd
@@ -0,0 +1,35 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/gcm.R
+\name{gcm}
+\alias{gcm}
+\title{GCM using random forests}
+\usage{
+gcm(Y, X, Z, alternative = c("two.sided", "less", "greater"), ...)
+}
+\arguments{
+\item{Y}{Response}
+
+\item{X}{Covariates}
+
+\item{Z}{Covariates}
+
+\item{alternative}{Alternative}
+
+\item{...}{Additional arguments to ranger}
+}
+\value{
+Object of class htest
+}
+\description{
+GCM using random forests
+}
+\examples{
+X <- matrix(rnorm(2e3), ncol = 2)
+colnames(X) <- c("X1", "X2")
+Z <- matrix(rnorm(2e3), ncol = 2)
+colnames(Z) <- c("Z1", "Z2")
+Y <- rnorm(1e3) # X[, 2] + Z[, 2] + rnorm(1e3)
+(gcm1 <- gcm(Y, X, Z))
+plot(gcm1)
+
+}
diff --git a/man/pcm.Rd b/man/pcm.Rd
new file mode 100644
index 0000000..240c2ca
--- /dev/null
+++ b/man/pcm.Rd
@@ -0,0 +1,57 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/pcm.R
+\name{pcm}
+\alias{pcm}
+\title{Conditional mean independence test}
+\usage{
+pcm(
+  Y,
+  X,
+  Z,
+  rep = 1,
+  est_vhat = TRUE,
+  reg = c("pcm_ranger", "pcm_lasso"),
+  ghat_args = NULL,
+  mtry = identity,
+  do.check = FALSE,
+  ...
+)
+}
+\arguments{
+\item{Y}{Numeric; response}
+
+\item{X}{Numeric; covariates}
+
+\item{Z}{Numceric; covariates}
+
+\item{rep}{Number of repetitions}
+
+\item{est_vhat}{Estimate variance functional}
+
+\item{reg}{Character; regression method}
+
+\item{ghat_args}{Arguments passed to reg}
+
+\item{mtry}{Argument passed to \code{ranger}}
+
+\item{do.check}{Save check data}
+
+\item{...}{Additional arguments passed to \code{reg}}
+}
+\value{
+Object of class '\code{htest}'
+}
+\description{
+Conditional mean independence test
+}
+\examples{
+X <- matrix(rnorm(2e3), ncol = 2)
+colnames(X) <- c("X1", "X2")
+Z <- matrix(rnorm(2e3), ncol = 2)
+colnames(Z) <- c("Z1", "Z2")
+Y <- rnorm(1e3) # X[, 2] + Z[, 2] + rnorm(1e3)
+(pcm1 <- pcm(Y, X, Z, ghat_args = list(mtry = NULL, max.depth = NULL),
+    est_vhat = TRUE, do.check = TRUE))
+plot(pcm1)
+
+}
diff --git a/tests/testthat.R b/tests/testthat.R
new file mode 100644
index 0000000..6ef959c
--- /dev/null
+++ b/tests/testthat.R
@@ -0,0 +1,12 @@
+# This file is part of the standard setup for testthat.
+# It is recommended that you do not modify it.
+#
+# Where should you do additional test configuration?
+# Learn more about the roles of various files in:
+# * https://r-pkgs.org/tests.html
+# * https://testthat.r-lib.org/reference/test_package.html#special-files
+
+library("testthat")
+library("comet")
+
+test_check("comet")
diff --git a/tests/testthat/test_main.R b/tests/testthat/test_main.R
new file mode 100644
index 0000000..81edfc3
--- /dev/null
+++ b/tests/testthat/test_main.R
@@ -0,0 +1,28 @@
+
+test_that(".ranger works", {
+  set.seed(12)
+  dat <- data.frame(bin = factor(sample(0:1, 1e3, TRUE)),
+                    ord = ordered(sample(7:10, 1e3, TRUE)),
+                    mcc = factor(sample(11:15, 1e3, TRUE)),
+                    num = rnorm(1e3),
+                    x = rnorm(1e3))
+  lapply(colnames(dat)[-ncol(dat)], \(resp) {
+    fm <- reformulate("x", resp)
+    rf <- .ranger(fm, data = dat)
+    rr <- residuals.ranger(rf)
+    expect_lt(abs(mean(rr)), 0.005)
+  })
+})
+
+test_that("gcm and pcm work", {
+  expect_no_error({
+    set.seed(12)
+    X <- matrix(rnorm(2e3), ncol = 2)
+    colnames(X) <- c("X1", "X2")
+    Z <- matrix(rnorm(2e3), ncol = 2)
+    colnames(Z) <- c("Z1", "Z2")
+    Y <- rnorm(1e3)
+    gcm1 <- gcm(Y, X, Z)
+    pcm1 <- pcm(Y, X, Z)
+  })
+})