# Insights for post-pandemic pedagogy across one CS department

Maxwell Bigman, Yosefa Gilon, Jenny Han, John Mitchell  mbigman, ygilon,
<EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract.

Adaptive remote instruction has led to important lessons for the future,
including rediscovery of known pedagogical principles in new contexts and new
insights for supporting remote learning. Studying one computer science
department that serves residential and remote undergraduate and graduate
students, we conducted interviews with stakeholders in the department (n=26)
and ran a department-wide student survey (n=102) during the four academic
quarters from spring 2020 to spring 2021. Our case study outlines what the
instructors did, summarizes what instructors and students say about courses
during this period, and provides recommendations for CS departments with
similar scope going forward. Specific insights address: (1) how instructional
components are best structured for students; (2) how students are assessed for
their learning; and (3) how students are supported in student-initiated
components of learning. The institution is a large U.S. research university
that has a history of online programs including online enrollment in regular
on-campus courses and large-scale open enrollment courses. Our recommendations
to instructors across the scope of this department may also be applicable to
other institutions that provide technology-supported in-person instruction,
remote enrollment, and hybrid courses combining both modalities.

CS education, remote learning, CS pedagogy CS departments, COVID-19, pandemic

††copyright: acmcopyright††conference: ; ; ††booktitle: UNPUBLISHED††ccs:
Applied computing Distance learning††ccs: Applied computing Interactive
learning environments

##  1\. Introduction

Colleges and computer science departments worldwide creatively responded to
the pressing need for remote instruction spurred by the global pandemic and
its consequences. As students and student expectations have been altered by a
traumatic and life-changing period, instructors have been open-minded and
resourceful in serving students successfully while confronting their own
challenges.

We believe several cultural and pedagogical changes spurred by the pandemic
will provide better learning opportunities for students as they return in
hybrid or in full. Hoping to help important shifts take hold on a department-
wide rather than course-by-course basis, we conducted a longitudinal best-
effort inventory of pandemic-era CS course experiences, combined with ongoing
department-wide pedagogical support and information exchange. While
significant changes have been observed at many institutions (Siegel et al.,
2021), the CS department represented in this case study both serves
approximately 20% of residential students and has a strong background in
online courses. Housed in a large research university with a computer science
teaching program at both undergraduate and graduate levels, a significant
remote student population enrolls in the same on-campus classes as residential
students. Over two decades, hybrid enrollment in many larger classes led to
widespread residential use of asynchronously accessible video. Many
instructors were also previously active in the MOOC movement from 2011 onward.
Thus both instructors and students were familiar with many tools of remote
instruction, including pre-recorded video, video recording of synchronous
class sessions, and online platforms that support submission, staff grading,
and automated grading of student work. Nonetheless, the pandemic led to many
shocking realizations.

Perhaps surprisingly, the instructors with online teaching experience who
approached pandemic instruction with the view that “we know how to do this”
encountered some of the greatest difficulties. Those who approached remote
instruction of displaced students as a new design problem appear to have fared
best, especially those who took the most care to monitor student conditions
and adapt their teaching accordingly. Looking forward from a disruptive two-
year period, our primary question is: What have we learned that can improve
our future programs, for enrolled students on and off campus?

In overview, Section 2 describes related work and Section 3 details our
approach. In Section 4, we offer a retrospective of teaching innovations from
2020 and 2021, organized into three main categories: (1) new formats for
instructor-directed learning, (2) alternative types of student assessment, and
(3) student-initiated learning, such as office hours or peer learning. This
section also highlights key changes in student expectations. In Sections 5 and
6, we shift toward the future, suggesting trajectories for post-pandemic
teaching and recommendations for each of the three categories of activities
highlighted in the paper.

##  2\. Related Work

Scholarship around online teaching and learning took off in the previous
decade with the advent of MOOCs, followed by a period of initial
disappointment as studies suggested online learning was not scaling
democratically nor reaching levels of engagement comparable to face-to-face
instruction (Kizilcec et al., 2013). When the COVID-19 pandemic provided a new
impetus to revisit online CS pedagogy (Kizilcec et al., 2021), instructors
explored interactive, human-centered online teaching practices (Latulipe,
2021; Piech et al., 2021), rethought the format and integrity of high-stakes
assessments (Ko, 2021), and redesigned their courses with student input (Bai
et al., 2020). A recent ITiCSE working group report underscored the fact that
we will likely never return to a pre-pandemic educational landscape (Siegel et
al., 2021). Because that working group provides an extensive literature
review, we refer the reader to their report (Siegel et al., 2021) rather than
provide a less extensive review in the limited space available here.

##  3\. Approach

Our case study for the purpose of instructional improvement focused on four
remote-learning academic quarters, Spring 2020 through Spring 2021
(inclusive). We carried out semi-structured interviews with stakeholders in
the department, tabulated course material, and conducted a department-wide
student survey in June 2021.

In total, we interviewed 19 instructors, 4 student teaching assistants (TAs)
or peer advisers, and representatives from 3 undergraduate interest groups
(Women in CS, American Indian Science and Engineering Society, Society of
Latinx Engineers). We sent out a general call to the department and selected
interviews for breadth among the instructors and TAs who volunteered their
time. We also conducted follow-up longitudinal interviews with 6 of the
instructors who had the most varied sets of insights and innovations.

The courses surveyed typically see enrollment of 50-400 students in any term,
including a small number of students (5–10%) who would have been remote
regardless of the pandemic. Notably, the courses represented various tracks
within the department, ranged from introductory core classes such as CS1 to
large-enrollment artificial intelligence classes to smaller human-computer
interaction design studios. In addition, we sought out interviews with student
groups, student TAs, and peer advisers, who advocated on behalf of a diverse
set of students in CS classes.

Semi-structured interviews lasted 40–60 minutes, depending on interviewee
interest. Questions were asked regarding course logistics, pedagogical
experiments, ongoing challenges, and initial reactions from students. The
final section of each interview followed the interest of the interviewee.

The department-wide survey to solicit students’ experiences with remote
learning received 102 participants. All undergraduate and graduate students
who had taken at least once CS class between Spring 2020 and Spring 2021 were
invited to fill out the survey, which was disseminated via email listservs.
The survey consisted of 5 Likert-scale or multiple choice questions, with
optional space to elaborate, and 2 two open response questions. We made clear
that participation was voluntary, optional, confidential, and uncompensated.
Although participation represented a small fraction of the current students,
we gained important insights from the instrument.

##  4\. Collected innovations

While our findings are consistent with pedagogical principles that were
understood pre-pandemic, several pedagogical changes have been normalized
fruitfully by a year of necessary changes. We have grouped the innovations
mentioned in the interviews into three categories and included quotations from
student survey responses where relevant.

###  4.1. Innovations in instructor-directed learning

####  4.1.1. Pre-pandemic

Most classes before the pandemic used a traditional lecture format, meeting
two to three times per week for 50-80 minutes each. Typically, student TAs
were assigned at a ratio of 20-30 students per assistant, with an additional
head TA allocated for larger courses over a few hundred students. TAs would be
involved in homework preparation, discussion or lab sections, office hours
providing direct student contact, and grading. Some courses were also
partially supported by a professional course coordinator who oversaw
management aspects of the course.

For decades, select large-classroom lectures from this department have been
recorded on video and made available online. With lecture material available
on demand, residential student attendance in lecture often dropped to around
15% after the third week of class. Office hours were often over-crowded, with
long lines.

####  4.1.2. Pandemic era

The main innovations explored were:

  * •

Disseminating content asynchronously through pre-recorded videos, available on
demand

  * •

Designing synchronous class time to complement pre-recorded material,
including

    * –

Fireside chats with guest speakers, high-level overviews,and Q&A

    * –

Problem sessions, often dividing the class in smaller groups or smaller
sections

    * –

Real-time questions answered by course assistants over chat or online
discussion platforms

  * •

Making synchronous lectures interactive with tools such as ohyay (ohy, [n.d.])

####  4.1.3. Instructor and student response

Students gave synchronous sessions mixed reviews, as shown in Figure 1. When
synchronous sessions were recorded, we often saw the same low attendance as in
our pre-pandemic in-person lectures. A limitation of our survey is that it
combines students from different classes, and therefore not all students
surveyed were exposed to the same course components in the same degree.

Figure 1. Survey responses to ”What aspects of synchronous class meetings did
you appreciate? Check all that apply.”

A few instructors who were redesigning synchronous online lectures began to
think of the experience as “more akin to going to a music performance.” The
instructors who used ohyay found that their attendance was higher in the
latter weeks of the course than it normally was when they were teaching in-
person, suggesting that there is “still value to the synchronous lecture” and
in fact “there are some aspects to better educating people in this [online]
space.” (Anonymized, 2021)

Eighty-five percent of students responding to the survey said that it was
“helpful” or “very helpful” to have lecture material pre-recorded and
available at any time. This suggests strong student demand for video recording
in the future but does not necessarily suggest that courses should be entirely
asynchronous.

In fact, students commented on the value of combining asynchronous and
synchronous learning: ”I really enjoyed having some asynchronous material and
then having synchronous lectures that delved deeper through real-world
examples once the concepts were out of the way.” Other students spoke about
the ability to ask online/chat question during class: ”TA and peer q and a
during lectures was infinitely helpful.” Another echoed the importance of
active learning: ”The in-class problem solving activities in [course redacted]
were great because we could anonymously message [name redacted] through zoom
chat and then through what he said see what others are thinking/feeling.”

These remarks highlight the importance of using synchronous class meetings as
a dynamic time for students to interact with their instructors to dive deeper
into topics that they were already somewhat familiar with based on the pre-
recorded lectures. The ability for students to ask more questions of their
instructors (particularly anonymously via Zoom chat) in real time, potentially
using the TAs to facilitate the Q&A, was an important innovation that moves
away from the previous model of largely one-way lectures.

For some courses, TA-led small group sections provided a crucial opportunity
for synchronous learning in the virtual setting. This builds on previous
research that underscores the importance of near-peer mentors as approachable
instructors (Forbes et al., 2017). Student remarks about these discussion
sections included: ”The only synchronous parts of class were section, which
was helpful for reviewing content and/or getting elaborations on stuff I
didn’t understand.” Additionally, another student said: ”I think it is
valuable to have sections meet synchronously. It was really valuable to
interact with my peers and digest the material learned together….As people
started to know each other, the activities and interactions became more and
more meaningful.”

###  4.2. Innovations in student assessment

####  4.2.1. Pre-pandemic.

Most classes before the pandemic used a combination of graded homework,
quizzes, programming projects, midterm exams and cumulative high-stakes final
exams. Grading has traditionally been competitive at this institution, with
students feeling considerable pressure and anxiety even though grade curves
often lead to a substantial number of high final grades.

####  4.2.2. Pandemic era.

University policy in spring 2020 required that all classes be taken pass/fail,
and the 2020-2021 school year prohibited high-stakes final exams and
eliminated the final exam week, although it should be recognized that some
instructors replicated their traditional exam structures despite official
policies. The primary innovations observed during 2020-21 were:

  * •

Multiple lower-stakes exams in lieu of one final exam

  * •

Revise and resubmit, allowing students to revise submitted work to receive a
higher grade (also known as mastery-based learning)

  * •

Pass/Fail grading policies

  * •

Partner or group assignments and exams

  * •

Concept checks, used regularly or weekly to track student academic progress
and personal wellbeing

####  4.2.3. Instructor and student response

Below, we consider the responses to each type of innovation:

Rethinking final exams. The remote learning environment prompted instructors
to recognize the unequal access to stable learning environments and resources
while students were under lockdown in their respective homes. For that reason,
instructors who traditionally relied on final exams pivoted towards lower-
stakes exams or a greater emphasis on assignments and projects.

One instructor explained that he would not use any exams during the pandemic
so as not to amplify disadvantages; he updated the grading policy to focus on
assignment grades instead. As he said: ”Right now… there are students who are
not outfitted to be able to take exams … Assuming that the people who are at a
disadvantage and taking exams right now are probably not the people who would
normally get top 5% on an exam, that they’re bringing a lot of their
shortcomings, even to campus,…to actually amplify that disadvantage is… just
not right.” He provided optional exams, which still gave students the
opportunity to self-assess their knowledge.

When surveyed, 70% of students considered it beneficial or very beneficial to
have multiple smaller exams in lieu of one final exam, and 69% of students
surveyed found partner / group assignments beneficial or very beneficial.
Students voiced appreciation for the changes above, which reduced stress: ”It
was helpful to not have exams in [course redacted] since I often struggle with
exams but definitely understood the material based on the self-assessments and
the assignments, so it just removed a lot of unnecessary stress.” Some
instructors who continued to give exams chose to accommodate more students by
giving a 24-hour window to complete exams or by giving shorter exams. By
giving shorter exams, some instructors also hoped to dissuade students from
cheating, which remained a concern in the virtual learning environment. Two
CS1 instructors reflected: ”To mitigate [cheating] we only gave them 25
minutes so it’s …not enough time for you to cheat.”

Revise and resubmit. For other instructors, the remote learning environment
proved to be an opportunity to pilot mastery-based learning, which includes
policies that allow students to revise and resubmit until they felt that they
had achieved mastery of what they learned. As one instructor explained: ”I’ve
always wanted to experiment with some form of mastery learning, and this was
my chance! We have 3 take-home exams, and after each take-home exam is graded,
students whose work was unsatisfactory have an opportunity (with extensive
feedback and staff support, though not the actual solutions) to revise and
resubmit their work a week later. Students whose work is still unsatisfactory
can revise & resubmit again, with more staff support.”

Students responded favorably to these policies: 75% of students who were in a
course that offered the option to revise and resubmit work considered it
beneficial or very beneficial to have the option to do so. The following
comments from students spoke to the advantages of a mastery-based approach
using a revise and resubmit policy: ”[Instructor redacted]’s revise and
resubmit policy in [course redacted] was fantastic. It meant I could engage
with topics until I reached mastery, rather than just settling for ’good
enough.’” The sentiment was echoed from other students: ”My learning
experience feels so much more whole now that I’ve worked out the correct
answer to every single midterm problem; I know this method of revision might
be controversial because it doesn’t provide a memorization-based meritocracy
grade distribution, but I feel like the CS program here has explained it’s
goal to be to make all of us excellent programmers through hard work - which I
think rewards revision, getting us students back on their feet after what
might’ve been a demoralizing fall from grace.

Pass / Fail Grading Policies. The ability to take courses for a Pass/Fail
grade instead of a letter grade changed students’ attitudes towards
assessment. Teaching team members of the CS1 courses reflected on why Pass-
Fail might be a helpful policy to level the playing field for introductory CS
courses: ”I think it’s not completely fair for intro CS classes to not be pass
fail because there are people who come in with different backgrounds … there
are people who start [CS1] who had never coded before in their life and I have
students who … had coded since they were in middle school.”

Weekly concept checks. A number of courses created short quizzes each week for
students to assess their understanding. 56% of students surveyed who had
weekly ”concept checks” in their courses considered it beneficial or very
beneficial to do them. When asked to elaborate on the use of “concept checks”
to accompany the asynchronous lectures in a flipped classroom environment, one
student stated: ”Weekly concept checks in [course redacted] also helped me
stay on track with the asynchronous lecture videos, and I thought that was a
great idea because I never got behind on lecture which I often do.”

Student appreciation for each of the assessment methods is reflected in Figure
2. Because some of the innovations were not widely deployed, a number of
respondents were not able to rate some of the methods, leading to lower scores
than might have been seen otherwise.

| Very beneficial | Beneficial | Somewhat ben. | Slightly ben. | Not ben.
---|---|---|---|---|---
Option to revise and resubmit work | 51.5% | 23.53% | 7.4% | 11.8% | 5.9%
Weekly concept checks | 30.6% | 25.88% | 20.0% | 11.8% | 11.8%
Multiple smaller exams in lieu of one final exam | 39.8% | 30.11% | 11.8% | 10.8% | 7.5%
Partner / group assignments or exams | 50.6% | 18.52% | 12.3% | 3.7% | 14.8%
Figure 2. Student response to assessment options

###  4.3. Innovations in student-initiated learning

####  4.3.1. Pre-pandemic.

Prior to the pandemic, office hours for larger lower-division classes were
traditionally held in a designated location on campus multiple times
throughout the week. The space, referred to as “the basement,” was often full
of student activity, with spontaneous group formation, peer help and small
group instruction with skilled TAs. Many students also had peers they could
rely on to work with through class connections or other relationships they had
built (e.g. student clubs).

####  4.3.2. Pandemic era.

In the shift to remote instruction, these student-initiated learning options
were arguably the hardest area to replace. Many of the ways that residential
life supports student learning were taken for granted before the pandemic.
Because the organizational structure of student peer collaboration was largely
invisible to instructors, many instructors began the pandemic with little
concrete understanding of how to replace residential peer learning. Because
there was no single solution for this complex problem, iterative experiments
with a range of platforms and techniques continued through the 2020-21
academic year, including:

  * •

Discussion forums, allowing students to ask questions and other students or
course staff to answer online

  * •

One-on-one office hours with course staff, as either drop-in, or 15 minute
appointments

  * •

Working office hours, allowing students to join online rooms to work with
other students, generally without course staff

  * •

”Homework Parties,” in which students can gather remotely to meet and work
together on Nooks(noo, [n.d.]), with course staff present at each session to
provide help

Unlike in the instructor-initiated activities and the assessments, experiments
in student-initiated learning struggled to replicate many of the in-person
benefits of office hours and informal campus study options. However, a number
of these formats were highly beneficial to the students that chose to take
advantage of them. This suggests that hybrid structures supporting student
collaboration might lead to productive innovations in the future.

####  4.3.3. Instructor and student response

Despite the challenges of hosting online student-initiated activities, it is
evident that students want to keep the option of online supports in the
future: ”Online OH [office hours] are so much more accessible and I actually
feel like I got a lot more help this year than in previous years. [The
basement] is chaotic and CAs are often running around like crazy and attention
is always divided. [Course redacted] group office hours functioned really well
online and my other classes … had individual online OH that worked super,
super well.” One graduate student compared the online option to the in-person
experience: ”Virtual OH with queues are, in my opinion, the best thing about
remote courses and I hope continue when things go back to normal… the virtual
queues on Nooks and Zoom made the process way more equitable … And probably
most crucially, if OH were crowded and there was a long queue, you could just
stay on the queue and comfortably do other work from your own home or location
as you waited to get helped.”

A number of students called out the benefits of using the Nooks platform as a
preferable alternative to Zoom for meeting and collaborating with groups of
peers during office hours. Finally, other students expressed appreciation for
the active role that instructors took on in online discussion spaces:
”[Instructor redacted] had a Slack workspace and he and the TA’s were super
active on it. I felt like I could receive and offer help to classmates, and I
really felt like a part of a community.”

##  5\. Synthesis and trajectories

In this section, we synthesize one or two main departmental priorities for
each of the three categories and connect them to trends reported in the
broader CS education literature.

###  5.1. Instructor-directed learning insights

The case study of our department suggests four key takeaways:

  1. (1)

Student appreciate when instructors make asynchronous lessons available prior
to synchronous meetings because it allows students to get better acquainted
with the material and maximizes the learning in synchronous sessions.

  2. (2)

Synchronous class time is most beneficial to students when designed around
“active learning” activities such as live problem practice, question and
answers with the teaching team, fireside chats with guest speakers, and
extending course material to apply it to real-world scenarios.

  3. (3)

Using technology to enable more opportunities for students to interact with
the teaching team and with peers makes learning more visible and creates more
opportunities for important discussions.

  4. (4)

Instructors who collect frequent student feedback using surveys and check-ins
have a better pulse on how their students are doing in the class and in their
personal lives, and this empathy translates to a better learning experience
for students.

Two of these changes speak to shifting trends over the past decade of CS
education that feel useful to elaborate upon.

Making asynchronous videos available. One clear trend is that students
strongly value the scheduling flexibility of asynchronous video presentations.
The ability to watch short lecture videos prior to synchronous course meetings
increases the value of the time spent as a group. Also, other institutions
have reported that the ability to rewatch and refer back to recordings of
synchronous meetings is useful for students (Culbert, 2021). These practices
are a nod to previous experiments with flipped classroom models (Horton et
al., 2014). While the research suggests that there is often low compliance to
watch asynchronous videos (Dazo et al., 2016), the recordings might still be
effective if they are paired with concept checks and/or if the videos’
relevance to synchronous meetings is properly communicated.

Recommendation: Make (short) pre-recorded videos available so that students
can cover the material at their own pace and can be better prepared for
synchronous meetings. This can also reduce the number of synchronous class
meetings per week.

New models of active learning in CS. Active learning teaching practices in
Physics (Crouch and Mazur, 2001), CS (Porter and Simon, 2013), and other STEM
disciplines have been around for quite a long time. While many instructors
were hesitant to change their practices and move away from traditional
lectures prior to the pandemic, the shift to remote instruction caused many to
rethink how to best use synchronous meetings. The experiments at our
institution did not reveal one “best way” to engage students during
synchronous class meetings. However, the four quarters of instruction did
reveal that students preferred classes in which they were (1) actively solving
practice problems, (2) asking questions of the instructor and the TAs (via
chat or audio), and (3) able to learn about real-world applications of
relevant course topics, including from guest speakers who are professionals in
the field.

These findings suggest that “active learning” might need to be more explicitly
defined in CS to encompass a broad range of activities in which students are
engaged with the instructor, the TAs, their classmates and/or guest speakers
during synchronous meetings. Making the course responsive to student needs and
questions, while creating opportunities to take advantage of the precious
human resources and opportunities for social interaction appears to greatly
benefit learning (Porter and Simon, 2013; Sarma and Yoquinto, 2020).

Recommendation: Instructors should structure synchronous class time to engage
students in activities, dialogue and practice problems that make learning more
social and provide opportunities for human interaction. Instructors should
survey their students and collect frequent feedback so that they can continue
to experiment and iterate with new formats for “active learning” methodologies
that are best suited to CS instruction in their courses.

###  5.2. Assessment insights

The choice by many institutions to require Pass/Fail grading in spring 2020 in
response to the pandemic forced instructors and students to approach grades
and assessment differently than they had in the past. Some students found it
liberating, while others were stressed about not knowing how to get an A in
the new class formats. The experiments in formative and summative assessments
highlighted a few key points of conversation within our department:

  1. (1)

Formative and summative assessments can take many forms. High-stakes exams are
not the only way to assess student learning, and more frequent exams and
concept checks might help students to assess how they are doing throughout a
course. Projects, group assignments and homeworks can all help to illuminate
student understanding.

  2. (2)

Revise and resubmit policies that support mastery of course material reduce
student stress and provide powerful opportunities for student learning and
success. A forced grading curve may not be necessary, particularly in
introductory courses.

  3. (3)

Some students appreciate grades and need graded transcripts to apply for jobs
and future degrees, while others are more focused on learning.

  4. (4)

Online exams are susceptible to cheating and misuse.

Formative and summative assessments take many forms. Important innovations
from the past experiments at our institution included the use of frequent
“concept checks” for students to self-assess their understanding of the
material as they moved through the course, rather than having to wait for a
higher-stakes midterm or final to find out how well they understood the
content. Similarly, more frequent but lower-stakes summative exams seemed to
reduce the stress on students and allow for faster learning cycles, although
more research is warranted. At other institutions, alternative assessments
proliferated as well, such as open-book exams (Quille et al., 2021), group
assessments (Ham and Myers, 2021) and oral exams (Sabin et al., 2021) with
generally positive results.

Recommendation: Use frequent assessment to reduce student stress and scaffold
learning effectively.

Mastery-Based Learning policies such as “revise and resubmit” support student
agency & learning. Students who had the opportunity to “revise and resubmit”
homework and assessments found that they were better able to learn the content
because they felt supported in continuing to learn the material until they
proved they had mastered it. Similar approaches, often termed mastery-based or
competency-based learning in education research (Voorhees, 2001), has been
shown to improve student outcomes by shifting the focus away from exams to
learning. The incredible response by students and instructors alike to these
policies suggest that both groups see tremendous benefits in this approach,
with opportunities to support learners without putting too much of a burden on
instructors (Lee, 2021).

Recommendation: Consider adopting policies that support mastery learning
methods.

###  5.3. Student-initiated learning insights

Shifting office hours and other essential yet often overlooked peer learning
structures online was a major challenge at our institution. The ability to
recreate the types of support systems that informally and spontaneously arose
in in-person settings was not readily available online. However, new systems
did emerge to reach students, especially those struggling in courses. Key
insights included:

  1. (1)

Both in-person and online office hours benefit students, and it is important
to have both options going forward as certain formats benefits certain types
of students.

  2. (2)

Creating opportunities for students to collaboratively work together and help
each other should be more deliberately designed for going forward. Leaving it
to chance advantages certain students (e.g. students with more friends in the
major, students who are more comfortable asking questions in class or in
office hours, etc.).

Looking to the future, it is evident that students at our institution are
excited to be back on campus in proximity to their peers and with access to
in-person academic spaces. However, the role of online support systems should
not be overlooked as a key resource for many students.

Recognize that students have different prior experiences; meet students where
they are. The strong desire to both return to in-person office hours and to
keep online office hours, as well as other new structures such as “homework
parties” suggests that different students learn in different way and need
different supports. For some who don’t feel comfortable in hectic environments
or in-person interactions, online office hours will continue to be their
preferred method for seeking help. For students that thrive in group
environments, an in-person setting might be preferred. For students without
strong social connections in their classes, facilitated “homework parties” can
provide invaluable access to peers. This takeaway echoes Universal Design for
Learning principles (Rose and Meyer, 2002) and other peer-learning models from
the pandemic (Blanton, 2021).

Recommendation: Meet students where they are by creating multiple formats for
students to get support. Consider different learning needs in both online and
in-person formats.

##  6\. Conclusions and Future Work

A descriptive case study of four academic quarters from spring 2020 to spring
2021 revealed energetic innovation and iterative experimentation, and student
comments strongly suggest that student expectations and demands have been
changed dramatically by this period. Remote learning at this institution has
at least partially continued to the time of writing. Given instructor
investment in the work they have done over this period, and the receptiveness
of students to many of these innovations, the most promising steps forward
will draw productively on the 2020-21 experience. Specific suggestions are
given in section 5, and future research is necessary to support the following
trends:

Trends in instructor-directed learning. Based on current and prior experience,
it is virtually certain that pre-recorded video and class session recordings
will continue to be valuable. However, with the low student attendance that
occurs when synchronous online sessions are recorded and accessible on demand,
there is a need to thoughtfully design synchronous course activities to foster
productive student engagement.

Trends in student assessment. During the pandemic, some instructors provided
more flexible deadlines to accommodate students in challenging circumstances,
replaced their single high-stakes final exam with alternatives, and explored
mastery-based learning. However, we predict that a nebulous transition period
will be necessary for all instructors and students to co-create and fully
adopt a new academic culture surrounding assessment.

Trends in student-initiated learning. Many of the ways that residential life
supports student learning were taken for granted before the pandemic but
became shockingly apparent as soon as we were confronted with their loss.
Because the organizational structure of student peer collaboration was largely
invisible to instructors, instructors began the pandemic with little concrete
understanding of how to replace residential peer learning. As there was no
single solution for this complex problem, iterative experiments with a range
of platforms and techniques continue.

A clear challenge for this department is to solidify and effectively expand
the use of pedagogy that was ushered in by the pandemic. We hope that the
simplicity of the recommendations highlighted in this paper will be useful in
focusing attention on a few selected incremental changes that are consistent
with current instructor and student views. Beyond those incremental steps, we
hope further that by encouraging instructors to monitor how well they support
students during the pandemic, we may have opened the door to longer-term
interest in adopting research-driving pedagogical practices that are
consistent with the local culture and values of the institution.

##  7\. Acknowledgements

We extend our deepest gratitude to the department’s instructors who responded
valiantly to the challenges of the past years with tremendous dedication to
their students. We are indebted to the instructors and students who
participated in our case study.

## References

  * (1)
  * noo ([n.d.]) [n.d.].  Nooks.  https://www.nooks.in/
  * ohy ([n.d.]) [n.d.].  Ohyay.  https://ohyay.co/
  * Anonymized (2021) Anonymized. 2021\.  Anonymized.  Anonymized
  * Bai et al. (2020) Yan Bai, Chunming Gao, and Bryan Goda. 2020.  Lessons Learned from Teaching Cybersecurity Courses During Covid-19. In _Proceedings of the 21st Annual Conference on Information Technology Education_. 308–313.
  * Blanton (2021) Maggie Blanton. 2021\.  Adapting STEM Peer Education Efforts Amidst a Pandemic. (2021).  https://whova.com/web/pprs_202105/Agenda Pandemic Pedagogy Research Symposium.
  * Crouch and Mazur (2001) Catherine H Crouch and Eric Mazur. 2001.  Peer instruction: Ten years of experience and results.  _American journal of physics_ 69, 9 (2001), 970–977.
  * Culbert (2021) Patrick D Culbert. 2021\.  When This Is All Over, Keep Recording Your Lectures.  https://www.chronicle.com/article/when-this-is-all-over-keep-recording-your-lectures
  * Dazo et al. (2016) Suzanne L. Dazo, Nicholas R. Stepanek, Robert Fulkerson, and Brian Dorn. 2016.  An Empirical Analysis of Video Viewing Behaviors in Flipped CS1 Courses. In _Proceedings of the 2016 ACM Conference on Innovation and Technology in Computer Science Education_ (Arequipa, Peru) _(ITiCSE ’16)_. Association for Computing Machinery, New York, NY, USA, 106–111.  https://doi.org/10.1145/2899415.2899468
  * Forbes et al. (2017) Jeffrey Forbes, David J. Malan, Heather Pon-Barry, Stuart Reges, and Mehran Sahami. 2017.  Scaling Introductory Courses Using Undergraduate Teaching Assistants. In _Proceedings of the 2017 ACM SIGCSE Technical Symposium on Computer Science Education_ (Seattle, Washington, USA) _(SIGCSE ’17)_. Association for Computing Machinery, New York, NY, USA, 657–658.  https://doi.org/10.1145/3017680.3017694
  * Ham and Myers (2021) Yeajin Ham and Brandon Myers. 2021.  Learning from Team Quizzes in CS2. In _Proceedings of the 52nd ACM Technical Symposium on Computer Science Education_ (Virtual Event, USA) _(SIGCSE ’21)_. Association for Computing Machinery, New York, NY, USA, 362–368.  https://doi.org/10.1145/3408877.3432412
  * Horton et al. (2014) Diane Horton, Michelle Craig, Jennifer Campbell, Paul Gries, and Daniel Zingaro. 2014.  Comparing Outcomes in Inverted and Traditional CS1. In _Proceedings of the 2014 Conference on Innovation & Technology in Computer Science Education_ (Uppsala, Sweden) _(ITiCSE ’14)_. Association for Computing Machinery, New York, NY, USA, 261–266.  https://doi.org/10.1145/2591708.2591752
  * Kizilcec et al. (2021) Rene F Kizilcec, Christos A Makridis, and Katharine C Sadowski. 2021.  Pandemic response policies’ democratizing effects on online learning.  _Proceedings of the National Academy of Sciences_ 118, 11 (2021).
  * Kizilcec et al. (2013) René F Kizilcec, Chris Piech, and Emily Schneider. 2013\.  Deconstructing disengagement: analyzing learner subpopulations in massive open online courses. In _Proceedings of the third international conference on learning analytics and knowledge_. 170–179.
  * Ko (2021) Melissa Ellen Ko. 2021\.  Revolutionizing Grading: Implications on Power, Agency, and Equity. In _2021 ASEE Virtual Annual Conference Content Access_.
  * Latulipe (2021) Celine Latulipe. 2021\.  A CS1 Team-Based Learning Space in Gather. Town. In _Proceedings of the 52nd ACM Technical Symposium on Computer Science Education_. 1245–1245.
  * Lee (2021) Cynthia Lee. 2021\.  Mastery Learning: Marrying High Standards and Flexibility Using a “Revise and Resubmit” Exam Policy.  https://ctl.stanford.edu/experiments-learning-series Experiments in Learning Series.
  * Piech et al. (2021) Christopher Piech, Ali Malik, Kylie Jue, and Mehran Sahami. 2021\.  Code in Place: Online Section Leading for Scalable Human-Centered Learning. In _Proceedings of the 52nd ACM Technical Symposium on Computer Science Education_. 973–979.
  * Porter and Simon (2013) Leo Porter and Beth Simon. 2013.  Retaining Nearly One-Third More Majors with a Trio of Instructional Best Practices in CS1. In _Proceeding of the 44th ACM Technical Symposium on Computer Science Education_ (Denver, Colorado, USA) _(SIGCSE ’13)_. Association for Computing Machinery, New York, NY, USA, 165–170.  https://doi.org/10.1145/2445196.2445248
  * Quille et al. (2021) Keith Quille, Keith Nolan, Brett A. Becker, and Seán McHugh. 2021.  _Developing an Open-Book Online Exam for Final Year Students_.  Association for Computing Machinery, New York, NY, USA, 338–344.  https://doi.org/10.1145/3430665.3456373
  * Rose and Meyer (2002) David H Rose and Anne Meyer. 2002.  _Teaching every student in the digital age: Universal design for learning._ ERIC.
  * Sabin et al. (2021) Mihaela Sabin, Karen H. Jin, and Adrienne Smith. 2021\.  _Oral Exams in Shift to Remote Learning_.  Association for Computing Machinery, New York, NY, USA, 666–672.  https://doi.org/10.1145/3408877.3432511
  * Sarma and Yoquinto (2020) Sanjay Sarma and Luke Yoquinto. 2020.  _Grasp: The Science Transforming How We Learn_.  Doubleday.
  * Siegel et al. (2021) Angela A Siegel, Mark Zarb, Bedour Alshaigy, Jeremiah Blanchard, Tom Crick, Richard Glassey, John R Hott, Celine Latulipe, Charles Riedesel, Mali Senapathi, et al. 2021\.  Teaching through a Global Pandemic: Educational Landscapes Before, During and After COVID-19.  In _Proceedings of the 2021 Working Group Reports on Innovation and Technology in Computer Science Education_. 1–25.
  * Voorhees (2001) Richard A Voorhees. 2001\.  Competency-Based learning models: A necessary future.  _New directions for institutional research_ 2001, 110 (2001), 5–13.


§ INTRODUCTION

When modeling spatial processes, it may be inappropriate to assume that
dependence is both stationary and isotropic, especially when dealing
with large domains. For example, when modeling US rainfall, we might
expect a different dependence structure over mountainous regions from
over plains. This work allows nonstationarity in dependence by deforming
the space on which a process is typically defined to one in which
stationarity and isotropy are more appropriate. Formally consider
\(\boldsymbol{x} \in \mathbb{R}^p\) and a mapping
\(\boldsymbol{g} : \mathbb{R}^p \mapsto \mathbb{R}^q\) such that
\(\boldsymbol{x}^* = \boldsymbol{g}(\boldsymbol{x})\) for
\(\boldsymbol{x}^* \in \mathbb{R}^q\). <cit.> introduced
this approach in the context of spatial modeling, so that \(p=q=2\):
then \(\boldsymbol{x} \in \mathbb{R}^2\), a coordinate in geographic
space (henceforth \(G\)-space), is mapped to
\(\boldsymbol{x}^* \in \mathbb{R}^2\), a coordinate in dispersion space
(henceforth \(D\)-space), where
\(\boldsymbol{x}^* = \boldsymbol{g}(\boldsymbol{x})\) for
\(\boldsymbol{g} : \mathbb{R}^2 \mapsto \mathbb{R}^2\). Stationarity and
isotropy are assumed for \(D\)-space. The remainder of this work focuses
on \(\boldsymbol{x} \in \mathbb{R}^2\).

<cit.> represent \(\boldsymbol{g}\) as a pair of thin plate
splines, which are fitted to coordinates obtained through non-metric
multidimensional scaling (NMDS). Subsequent works by <cit.>
and <cit.> adopt a Bayesian approach to inference and assume
that \(\boldsymbol{g}\) is a random function and data are realizations
of a Gaussian process (GP); i.e.,

\begin{equation} \label{gp} Y_t(\boldsymbol{x}) \mid \boldsymbol{g} \sim GP\big(\mu(\boldsymbol{x}), v\big(\boldsymbol{g} (\boldsymbol{x}), .\big)\big)\end{equation}

for time \(t=1, \ldots, n\), location
\(\boldsymbol{x} \in \mathbb{R}^2\), and mean and covariance functions
\(\mu\) and \(v\), respectively. For a fixed set of locations,
\(\boldsymbol{x}_1, \ldots, \boldsymbol{x}_m\), this allows inference to
be based on the likelihood

\begin{equation}L(\boldsymbol{\Sigma}) = |2 \pi \boldsymbol{\Sigma}|^{-(n - 1) / 2} \exp\Big(-\dfrac{n}{2}\text{tr}(\boldsymbol{\Sigma}^{-1} \boldsymbol{V})\Big),\label{gplik}\end{equation}

\(\boldsymbol{V} = (n - 1)^{-1}\sum_{t=1}^n(\boldsymbol{y}_t - \hat{\boldsymbol{\mu}})(\boldsymbol{y}_t - \hat{\boldsymbol{\mu}})^\text{T}\)
\(\boldsymbol{y}_t = (y_t(\boldsymbol{x}_1), \ldots, y_t(\boldsymbol{x}_m))\),
\(m \times m\) matrix \(\boldsymbol{\Sigma}\) has \((i, j)\)th element
\(\Sigma_{i, j} = v(\boldsymbol{g} (\boldsymbol{x}_i), \boldsymbol{g} (\boldsymbol{x}_j))\)
and \(\hat{\boldsymbol{\mu}} = (\hat \mu_1, \ldots, \hat \mu_m)\) with
\(\hat \mu_i = n^{-1}\sum_{t=1}^n y_t(\boldsymbol{x}_i)\).
<cit.> and <cit.> place thin plate spline and
Gaussian process priors on \(\boldsymbol{g}\), respectively, and use
Markov chain Monte Carlo to sample from the posterior distribution of
\(\boldsymbol{g}\), which readily allows its uncertainty to be

Spatial deformation models suffer the intuitively undesirable flaw of
allowing \(D\)-spaces that `fold', i.e., non-bijective mappings
\(\boldsymbol{g}\) or, more conceptually, mappings such that for every
\(\boldsymbol{x}\) in \(G\)-space there is not a single
\(\boldsymbol{x}^*\) in \(D\)-space. <cit.> hinder folding by
considering the `bending energy' of \(D\)-spaces, and increasingly
penalize spaces that require less energy to bend. Alternatively,
<cit.> propose to represent \(\boldsymbol{g}\) as a
multivariate GP, and state that `the GP formulation for
[\(\boldsymbol{g}\)] tends to eliminate the kind of non-injective
mappings that were noted by <cit.>'.

Various approaches have explicitly addressed avoiding mappings that
fold. <cit.> use a Delaunay triangulation of \(G\)-space to
identify and eliminate mappings that give folds in \(D\)-space.
<cit.> derive conditions on deformations based on radial
basis functions that avoid folds. Nonstationarity in dependence is also
considered when emulating computer models and referred to as input
warping (IW): the computer model's inputs are transformed to a scale on
which dependence is stationary; see, e.g., <cit.>.
<cit.> propose deep compositional spatial models for
representing \(\boldsymbol{g}\) in which the compositional formulation
can ensure bijectivity. <cit.> propose to represent the
compositions through IW GPs and deep stochastic processes (DSPs), both
of which are based on basis representations with weights and unknown
basis function parameters. The IW GPs have unknown weights, which are
estimated by maximum likelihood, whereas the DSPs have random weights,
which are assumed to be of log-Gaussian form and estimated by
variational Bayes. <cit.> then propose three approaches to
warping: axial warping units, with positive weights and monotonic basis
functions; radial basis functions, employing the constraints of
<cit.>; and Möbius transformation units, which make analogy
between mapping from \(\mathbb{C}\) to itself with mapping from
\(\mathbb{R}^2\) to itself. Each can be used with IW GPs and DSPs and
ensure bijectivity by virtue of the compositional structure.
<cit.> represent \(D\)-space as the tensor product of B-splines
and derive constraints on the knots that, when imposed, give
bijectivity. For a wide and recent review of statistical modeling of
nonstationary covariance – including deformation-based approaches –
see <cit.>.

Folding could be considered a consequence of a two-dimensional
\(D\)-space being insufficient to bring isotropy. <cit.> and
<cit.> propose extending \(D\)-space to \(2 + r\) dimensions
for \(r \geq 1\) so that
\(\boldsymbol{g} : \mathbb{R}^{2} \mapsto \mathbb{R}^{2 + r}\).
<cit.> refer to this approach as dimension expansion.
<cit.> place a GP prior on \(\boldsymbol{g}\), allow
covariates in \(v\) and then base \(v\) on Mahalanobis distance to
generalize the usual Euclidean distance. <cit.> estimate the
latent dimensions in a two-stage procedure that finds interim values
using a least squares fit between empirical and model-based variograms,
which are then approximated using thin plate splines.

This article introduces the deform R package, which is
designed to fit nonstationary spatial GPs using either spatial
deformations or dimension expansions to represent \(\boldsymbol{g}\). In
particular, it facilitates the visualization of fitted models; for
example, we may want to inspect \(\boldsymbol{g}\) for any physical

Various R packages exist that fit geostatistical models. A
selection include packages fields <cit.>, spatial
<cit.> and LatticeKrig <cit.>, which can
perform kriging, and geoR <cit.>, gstat
<cit.>, spmodel <cit.> and
RandomFields <cit.> give functions
for fitting stationary GPs. Extensions for fitting nonstationary GPs
include FRK <cit.>, tgp <cit.> and
enviroStat <cit.>, which implements <cit.>'s
method (but is no longer supported by CRAN). Alternatively, Gaussian
Markov random fields can be fit in mgcv with basis specification
(see <cit.>) and with
R-INLA, <cit.>,
which uses the integrated nested Laplace approximation of
<cit.>. R packages for non-Gaussian geostatistical
data also exist, such as SpatialExtremes <cit.>
and CARBayes <cit.>. The fitting of nonstationary GPs to
spatial data can be achieved through various software packages for DGPs,
such as dgpsi <cit.> for R, and various
python options, such as GPy <cit.> and
GPflow <cit.>.

The next section of this article introduces low-rank representations for
\(\boldsymbol{g}\) based on thin plate regression splines
<cit.>. These apply to both the spatial deformation approach
of <cit.>, which is extended to ensure bijectivity, and to
the dimension expansion approach of <cit.>. Section 3
introduces objective methods for inference for such models. Section 4
introduces the key functions in package deform. Section 5
demonstrates deform's functions on the solar radiation data
originally used in <cit.>. Section 6 covers the modeling of
censored data and then section 7 summarizes the work presented.

§ METHODOLOGY

In the following section, consider again \(Y_t(\boldsymbol{x})\), a
process at time \(t=1, \ldots, n\) and location
\(\boldsymbol{x} = (x_1, x_2) \in \mathcal{G}\). Specifically,
deform only considers the fitting of zero-mean GPs, so
\(\mu(\boldsymbol{x}) = 0\) in (<ref>). Forms assumed for
\(v(\boldsymbol{g}(\boldsymbol{x}), \boldsymbol{g}(\boldsymbol{x}'))\)
in deform are given in <ref>. In the two-dimensional
case, \(x_1\) and \(x_2\) are longitude and latitude coordinates,
respectively. Then consider the transformation
\(\boldsymbol{x}^* =\boldsymbol{g}(\boldsymbol{x})\), where
\(\boldsymbol{x}\) exists in \(G\)-space and \(\boldsymbol{x}^*\) exists
in \(D\)-space. The methodology presented readily extends to
\(G\)-spaces defined over any number of dimensions, as in
<cit.>. Independence over time will be assumed to focus on
spatial dependence.

§.§ Spatial deformation

A spatial deformation is defined here as the transformation
\(\boldsymbol{x}^* =\boldsymbol{g}(\boldsymbol{x})\), given a mapping
\(\boldsymbol{g}\), where
\(\boldsymbol{g}:\mathbb{R}^2 \mapsto \mathbb{R}^2\),
\(\boldsymbol{x}^* = (x_1^*, x_2^*) = (g_1(\boldsymbol{x}), g_2(\boldsymbol{x}))\)
and \(g_d : \mathbb{R}^2 \mapsto \mathbb{R}\) for \(d = 1, 2\). For
spatial deformations, <cit.>'s parametrization is adopted, so
that for
\(\boldsymbol{x}^* = (x_1^*, x_2^*) = (g_1(\boldsymbol{x}), g_2(\boldsymbol{x}))\),

\begin{align}
\label{g1} g_1(\boldsymbol{x}) = \exp(\alpha_1) x_1 + \alpha_3 x_2 + \sum_{i = 1}^m \delta_{1i} \eta(\norm{\boldsymbol{x} - \boldsymbol{x}_i}),\\
\label{g2} g_2(\boldsymbol{x}) =  \alpha_3 x_1 + \exp(\alpha_2) x_2 + \sum_{i = 1}^m \delta_{2i} \eta(\norm{\boldsymbol{x} - \boldsymbol{x}_i}),
\end{align}
where \(\eta(l) = l^2 \log l\) with
\(\norm{\mathbf{x} - \mathbf{x}'} = \sqrt{(x_1 - x_1')^2 + (x_2 - x_2')^2}\)
and the above equations are subject to the constraints
\(\mathbf{T}^\text{T} \boldsymbol{\delta}_1 = \mathbf{T}^\text{T} \boldsymbol{\delta}_2 = \mathbf{0}_{3 \times 1}\),
where \(m \times 3\) matrix \(\mathbf{T}\) has \(i\)th row
\((1, x_{i1}, x_{i2})\), for \(i = 1, \ldots, m\),
\(\boldsymbol{\delta}_1 = (\delta_{11}, \ldots, \delta_{1m})^\text{T}\),
\(\boldsymbol{\delta}_2 = (\delta_{21}, \ldots, \delta_{2m})^\text{T}\)
and \(\mathbf{0}_{n \times m}\) denotes a \(n \times m\) matrix
comprizing only zeros. These constraints avoid over-parametrization,
while the coefficients of \(x_1\) and \(x_2\) in equations (<ref>)
and (<ref>) avoid rotational invariance. Following <cit.>,
let matrix \(\mathbf{E}\) have \((i, j)\)th element
\(E_{ij} = \eta(\norm{\boldsymbol{x}_i - \boldsymbol{x}_j})\) and
\(\mathbf{E} = \mathbf{U} \boldsymbol{\Lambda} \mathbf{U}^\text{T}\).
Then let
\(\mathbf{E}_k = \mathbf{U}_k \boldsymbol{\Lambda}_k \mathbf{U}_k^\text{T}\),
where \(\mathbf{U}_k\) denotes the first \(k\) columns of \(\mathbf{U}\)
and \(\boldsymbol{\Lambda}_k\) denotes the upper \(k \times k\) block of
\(\boldsymbol{\Lambda}\), with
\(\boldsymbol{\Lambda} = \text{diag}(e_1, \ldots, e_m)\) comprizing the
\(m\) eigenvalues of \(\boldsymbol{\Lambda}\) arranged as
\(e_1 \geq \ldots \geq e_m\). Attention can then be restricted to the
space spanned by \(\mathbf{U}_k\), i.e., to
\(\boldsymbol{\delta}_1^{(k_1)}\) such that
\(\boldsymbol{\delta}_1 = \mathbf{U}_{k_1} \boldsymbol{\delta}_1^{(k_1)}\)
and \(\boldsymbol{\delta}_2^{(k_2)}\) such that
\(\boldsymbol{\delta}_2 = \mathbf{U}_{k_2} \boldsymbol{\delta}_2^{(k_2)}\),
where \(k_1\) and \(k_2\) are the finite ranks chosen to represent
\(g_1\) and \(g_2\), respectively (although often \(k_1 = k_2\) will be
used in practice). An unconstrained optimization problem can be formed
from the constraints
\(\mathbf{T}^\text{T} \boldsymbol{\delta}_1 = \mathbf{T}^\text{T} \boldsymbol{\delta}_2 = \mathbf{0}_{3 \times 1}\)
for the rank-\(k\) case by choosing some matrix \(\mathbf{Z}_k\) such
\(\mathbf{T}^\text{T} \mathbf{U}_k \mathbf{Z}_k = \mathbf{0}_{k \times (k - 3)}\).
This can be achieved via the QR-decomposition,
\(\mathbf{T}^\text{T} \mathbf{U}_k = \mathbf{QR}\), say, taking
\(\mathbf{Z}_k\) as the final \(k - 3\) columns of \(\mathbf{Q}\). For
the unconstrained problem, \(\tilde{\boldsymbol{\delta}}_1\) and
\(\tilde{\boldsymbol{\delta}}_2\) satisfying
\(\boldsymbol{\delta}_1^{(k_1)} = \mathbf{Z}_{k_1} \tilde{\boldsymbol{\delta}}_1\)
\(\boldsymbol{\delta}_2^{(k_2)} = \mathbf{Z}_{k_2} \tilde{\boldsymbol{\delta}}_2\)
can be used. The unknown parameters that define the deformation are
given by the \((k_1 + k_2 - 3)\)-vector
\(\boldsymbol{\beta} = (\alpha_1, \alpha_2, \alpha_3, \tilde{\boldsymbol{\delta}}_1^\text{T}, \tilde{\boldsymbol{\delta}}_2^\text{T})^\text{T}\).

Typically an additive wiggliness penalty is imposed on \(g_1\) and
\(g_2\) of the form
\(\lambda_1 \boldsymbol{\delta}_1^\text{T} \mathbf{E} \boldsymbol{\delta}_1 + \lambda_2 \boldsymbol{\delta}_2^\text{T} \mathbf{E} \boldsymbol{\delta}_2\).
This becomes
\(\lambda_1 (\boldsymbol{\delta}_1^{(k_1)})^\text{T} \boldsymbol{\Lambda}_{k_1} \boldsymbol{\delta}_1^{(k_1)} + \lambda_2 (\boldsymbol{\delta}_2^{(k_2)})^\text{T} \boldsymbol{\Lambda}_{k_2} \boldsymbol{\delta}_2^{(k_2)}\),
allowing for the finite-rank representations of \(g_1\) and \(g_2\), and
\(\lambda_1 (\tilde{\boldsymbol{\delta}}_1^{(k_1)})^\text{T} \mathbf{Z}_{k_1}^\text{T} \boldsymbol{\Lambda}_{k_1} \mathbf{Z}_{k_1} \tilde{\boldsymbol{\delta}}_1^{(k_1)} + \lambda_2 (\tilde{\boldsymbol{\delta}}_2^{(k_2)})^\text{T} \mathbf{Z}_{k_2}^\text{T} \boldsymbol{\Lambda}_{k_2} \mathbf{Z}_{k_2} \tilde{\boldsymbol{\delta}}_2^{(k_2)}\),
allowing for the unconstrained optimization problem. Let
\(\boldsymbol{S}_{\boldsymbol{\lambda}} = \lambda_1 \boldsymbol{S}_1 + \lambda_2 \boldsymbol{S}_2\),
\(\boldsymbol{S}_1 = \text{diag}(\mathbf{0}_{3 \times 3}, \allowbreak \mathbf{Z}_{k_1}^\text{T} \boldsymbol{\Lambda}_{k_1} \mathbf{Z}_{k_1}, \allowbreak \mathbf{0}_{(k_2 - 3) \times (k_2 - 3)})\)
\(\boldsymbol{S}_2 = \text{diag}(\mathbf{0}_{k_1 \times k_1}, \mathbf{Z}_{k_2}^\text{T} \boldsymbol{\Lambda}_{k_2} \mathbf{Z}_{k_2})\),
where \(\boldsymbol{\lambda} = (\lambda_1, \lambda_2)\). Given
\(\boldsymbol{\beta}\) above, the wiggliness penalty can be written
\(\boldsymbol{\beta}^\text{T} \boldsymbol{S}_{\boldsymbol{\lambda}} \boldsymbol{\beta}\),
which is of the form presented in <cit.> and hence used in
Section <ref>. Fitting the spatial deformation model is achieved
by estimating \(\boldsymbol{\beta}\) subject to the penalty
\(\boldsymbol{\beta}^\text{T} \boldsymbol{S}_{\boldsymbol{\lambda}} \boldsymbol{\beta}\).

§.§ Bijective spatial deformation

<cit.> demonstrate how smoothing parameters associated with
thin plate splines, or equivalently that control the bending energy, may
be fixed to ensure bijectivity. This work aims to maintain objectivity
by allowing optimal estimation of smoothing parameters under the
condition of bijectivity. <cit.> ensure bijectivity by
representing \(G\)-space as a Delaunay triangulation, which, when
transformed to \(D\)-space, is bijective if none of the vertices lie
within any of the triangles. The approach of <cit.> applies
to any form for \(\boldsymbol{g}\), which is a criterion that the
approach proposed here also satisfies. The approaches of
<cit.> and <cit.> require specific—albeit
seemingly rather flexible—forms for \(\boldsymbol{g}\).

Clockwise areas of triangles. Left: A triangle represented as clockwise points $(x_{11}, x_{12})$, $(x_{21}, x_{22})$, $(x_{31}, x_{32})$ with clockwise area 0.5. Right: A transformation to the left-hand triangle, in which $(x_{31}, x_{32}) \mapsto (x_{31}^*, x_{32}^*)$, giving a negative clockwise area of -0.5, based on the clockwise ordering in the left-hand triangle. Such negative areas are used to identify grids that have folded.

A related approach to <cit.> is proposed here in which
\(\mathcal{G}\), the domain of interest, is represented as a triangular
tiling. The clockwise area of each triangle is computed, which, based on
Figure <ref>, is given by
\((x_{21} x_{12} + x_{31} x_{22} + x_{11} x_{32} - x_{11} x_{22} - x_{21} x_{32} - x_{31} x_{12})/2\),
where \((x_{i1}, x_{i2})\), \(i=1, 2, 3\), are vertices of a triangle
defined in clockwise order. Subject to the triangular tiling's finite
representation of \(\mathcal{G}\), a change in ordering can be used to
identify non-bijective \(\boldsymbol{g}\), which is equivalent to
\(\boldsymbol{g}\) turning a triangle's clockwise area negative. This is
illustrated in Figure <ref> in which the left-hand triangle
has clockwise area 0.5 whereas the right-hand triangle has clockwise
area \(-0.5\). A space represented by a triangular tiling (see Figure
<ref>, row 1, column 1) with a mixture of positive and negative
areas must have folded (see Figure <ref>, row 2, column 2); all
positive areas corresponds to a fold-free space; and all negative areas
corresponds to a fold-free space that has `flipped'. Flipped spaces can
be eliminated without loss of generality since equivalent distances for
such spaces can be achieved if the space is flipped back.

Representations of spatial deformations using triangular tilings. Row 1, column 1: The original triangular tiling of the domain-spanning grid for calculating clockwise triangle areas. Row 1, column 2: An example of an affine transform of the form $\boldsymbol{x} \mapsto {\bf M} \boldsymbol{x}$, for a $2 \times 2$ matrix ${\bf M}$ and $\boldsymbol{x} = (x_1, x_2)'$. Row 2: Examples of deformations achieved by $\boldsymbol{x} \mapsto \boldsymbol{g}(\boldsymbol{x})$, where $\boldsymbol{g}$ comprises two thin plate regression splines, demonstrating fold-free (column 1) and folded (column 2) deformations.

Consider the triangular tiling
\(\mathcal{G} = \cup_{l=1}^L \mathcal{W}_l\), where each
\(\mathcal{W}_l\), for \(l=1, \ldots, L\), is a triangle with clockwise
area \(A(\mathcal{W}_l)\). Imposing the penalty
\(\delta I\big[\min\{A(\mathcal{W}_1), \ldots, A(\mathcal{W}_L)\} < 0\big]\)
could be used to eliminate folds, for large \(\delta > 0\) and where
\(I(\,)\) is the indicator function. In practice this can work poorly
with the inference approach proposed in <ref> as it can result
in parameter estimates at a non-differentiable part of parameter space.
Instead, a penalty of the form
\begin{equation}
\label{foldpen}
\frac{\delta}{2} \big[A_-(\mathcal{W}_1 / \epsilon) + \ldots + A_-(\mathcal{W}_L / \epsilon)\big]^2
\end{equation}
is chosen, where \(A_-(W) = \min\{A(W), 0\}\) and
\(\epsilon > 0\), which increasingly penalizes deformations with more
folding and is twice-differentiable, which aids inference. User-friendly
choice of \(\delta\) and \(\epsilon\) is discussed in

§.§ Dimension expansion

Dimension expansions have
\(\boldsymbol{x}^* =\boldsymbol{g}(\boldsymbol{x})\) with
\(\boldsymbol{g}:\mathbb{R}^2 \mapsto \mathbb{R}^{2 + r}\) where
\(\boldsymbol{x}^* = (g_1(\boldsymbol{x}), \ldots, g_{2 + r}(\boldsymbol{x}))\),
with \(g_d : \mathbb{R}^2 \mapsto \mathbb{R}\) for
\(d = 1, \ldots, r + 2\). The parametrization for dimension expansions
is a slight extension to that of <cit.>, so that
\(\boldsymbol{x}^* = (g_1(\boldsymbol{x}), g_2(\boldsymbol{x}), \ldots, g_{r + 2}(\boldsymbol{x}))\),
\begin{align*}
g_1(\boldsymbol{x}) &= \exp(\alpha_1) x_1 + \alpha_3 x_2,\\
g_2(\boldsymbol{x}) &=  \alpha_3 x_1 + \exp(\alpha_2) x_2,\\
g_{d + 2}(\boldsymbol{x}) &= \sum_{i = 1}^m \delta_{di} \eta(\norm{\boldsymbol{x} - \boldsymbol{x}_i}),
\end{align*}
for \(d = 1, \ldots, r\), with \(\eta\) as in
<ref>. As in <ref>, matrix
\(\mathbf{E}\) with \((i, j)\)th element
\(E_{ij} = \eta(\norm{\boldsymbol{x}_i - \boldsymbol{x}_j})\) can be
formed with eigen-decomposition
\(\mathbf{E} = \mathbf{U}^\text{T} \boldsymbol{\Lambda} \mathbf{U}\),
which can be used to give a rank-\(k_d\) approximation to
\(g_{d + 2}(\boldsymbol{x})\). For dimension expansions, the usual
affine transformation of the thin plate spline is omitted from
\(g_{d + 2}\) here, for \(d = 1, \ldots, r\); instead a scale and
rotation transformation is captured with \(g_1\) and \(g_2\). As a
result, parameters of \(g_{d + 2}\) for \(d = 1, \ldots, r\) are
\(\tilde{\boldsymbol{\delta}}_d\) satisfying
\(\boldsymbol{\delta}_d = \mathbf{U}_{k_d} \tilde{\boldsymbol{\delta}}\),
with \(\mathbf{U}_{k_d}\) the first \(k_d\) columns of \(\mathbf{U}\)
and where
\(\boldsymbol{\delta}_d = (\delta_{d1}, \ldots, \delta_{dm})^\text{T}\).
The unknown parameters that define the dimension expansion are given by
the \((3 + K_r)\)-vector
\(\boldsymbol{\beta} = (\alpha_1, \alpha_2, \alpha_3, \tilde{\boldsymbol{\delta}}_1^\text{T}, \ldots, \tilde{\boldsymbol{\delta}}_r^\text{T})^\text{T}\),
where \(K_d = \sum_{j = 1}^d k_j\). A wiggliness penalty can be imposed
on each \(g_{d + 2}\) of the form
\(\lambda_d \tilde{\boldsymbol{\delta}}_d^\text{T} \boldsymbol{\Lambda}_{k_d} \tilde{\boldsymbol{\delta}}_d\)
such that \(\lambda_d \to \infty\) gives \(g_d(\boldsymbol{x}) \to 0\),
for \(d = 1, \ldots, r\); see <cit.>. Large values of
\(\lambda_d\) can be used to identify and then eliminate redundant
dimensions. Penalties can be combined as
\(\boldsymbol{\beta}^\text{T} \boldsymbol{S}_{\boldsymbol{\lambda}} \boldsymbol{\beta}\),
\(\boldsymbol{S}_{\boldsymbol{\lambda}} = \lambda_1 \boldsymbol{S}_1 + \ldots + \lambda_r \boldsymbol{S}_r\),
where \(\boldsymbol{\lambda} = (\lambda_1, \ldots, \lambda_r)\) and
\(\boldsymbol{S}_d\) is a \((3 + K_r) \times (3 + K_r)\) matrix
\(\boldsymbol{S}_d\) comprizing entirely zeros except for
\((e_1, \ldots e_{k_d})\) on its leading diagonal at positions
\(S_{4 + K_{d - 1}, 4 + K_{d - 1}}, \ldots, S_{3 + K_d, 3 + K_d}\),
respectively, with \(K_0 \equiv 0\).

§.§ Notes on finite-rank basis representations

<cit.>, <cit.> and <cit.> have
previously used thin plate splines to define locations in \(D\)-space.
This work proposes to use thin plate regression splines, which
are a finite-rank approximation based on a truncated eigen-decomposition
of corresponding thin plate splines, with various optimality criteria:
see <cit.>. Instead of thin plate splines, the tensor product
of two one-dimensional splines could be used to form the two-dimensional
transformation from \(G\)- to \(D\)-space. More generally, the proposed
framework allows any \(g(\boldsymbol{x})\) function to be treated as a
smooth, and therefore can be characterized similarly to smooths in a
generalized additive models; see, e.g., <cit.>. Such
specifications would need careful consideration in order to avoid
over-parametrization and rotational invariance.

§ INFERENCE

§.§ Data and notation

Estimation of the spatial deformation or dimension expansion models will
be presented for a fixed set of locations,
\(\mathcal{X} = \{\boldsymbol{x}_1, \ldots, \boldsymbol{x}_m\}\), and a
fixed set of time points, \(\mathcal{T} = \{1, \ldots, n\}\).
Corresponding data are
\(\boldsymbol{y} = (\boldsymbol{y}_1, \ldots, \boldsymbol{y}_n)\), with
\(\boldsymbol{y}_t\) for \(t = 1, \ldots, n\) as defined in
<ref>. Note that deform absorbs any parameters required
by the covariance function \(v\) in relation (<ref>) into
\(\boldsymbol{\beta}\) and pads
\(\boldsymbol{S}_{\boldsymbol{\lambda}}\) with zeros accordingly.
Fitting models in deform therefore involves estimating
coefficients \(\boldsymbol{\beta}\), which determine the \(G\)- to
\(D\)-space mapping \(\boldsymbol{g}\), and smoothing parameters,
\(\boldsymbol{\lambda}\). Finally let \(\ell(\boldsymbol{\beta})\)
denote the log-likelihood of a zero-mean GP, as defined in equation

§.§ Restricted maximum likelihood

For \(\boldsymbol{\beta}\) given \(\boldsymbol{\lambda}\), let
\(\ell_p(\boldsymbol{\beta}, \boldsymbol{\lambda})\) denote the
penalized log-likelihood given by

\begin{equation*} \label{penalized} \ell_p(\boldsymbol{\beta}, \boldsymbol{\lambda}) = \ell(\boldsymbol{\beta}) - \frac{1}{2}\boldsymbol{\beta}^T \boldsymbol{S}_{\boldsymbol{\lambda}} \boldsymbol{\beta}.\end{equation*}

In the case of bijective spatial deformations, this may be further
penalized with the penalty of equation (<ref>) to restrict
smoothing parameters \(\boldsymbol{\lambda}\) to those that lead to
bijective \(\boldsymbol{g}\).

Smoothing parameters, \(\boldsymbol{\lambda}\), are estimated
objectively using restricted maximum likelihood (REML), as outlined in
<cit.> and <cit.>. This results from recognizing
that the penalized likelihood's penalty is proportional to the exponent
of a MVN(\({\bf 0}, \boldsymbol{S}_{\boldsymbol{\lambda}}^-)\)
distribution, where \(\boldsymbol{S}_{\boldsymbol{\lambda}}^-\) denotes
the generalized inverse of \(\boldsymbol{S}_{\boldsymbol{\lambda}}\).
Then \(\boldsymbol{\beta}\) is treated as a vector of random effects and
integrated out by Laplace's method. The resulting restricted
log-likelihood takes the form

\begin{equation*} \label{ell} \ell(\boldsymbol{\lambda}) = \ell_p(\boldsymbol{\beta}_{\boldsymbol{\lambda}}, \boldsymbol{\lambda}) + \frac{1}{2} \log |\boldsymbol{S}_{\boldsymbol{\lambda}}|_+ - \frac{1}{2} \log|{\bf H}(\hat{\boldsymbol{\beta}}_{\boldsymbol{\lambda}})| + \frac{M_p}{2} \log(2 \pi),\end{equation*}

where \(|\boldsymbol{S}_{\boldsymbol{\lambda}}|_+\) denotes the product
of positive eigenvalues of \(\boldsymbol{S}_{\boldsymbol{\lambda}}\),
\(\hat{\boldsymbol{\beta}}_{\boldsymbol{\lambda}}\) maximises
\(\ell_p(\boldsymbol{\beta}, \boldsymbol{\lambda})\) w.r.t.
\(\boldsymbol{\beta}\),
\({\bf H}(\hat{\boldsymbol{\beta}}_{\boldsymbol{\lambda}})\) is the
negative Hessian of \(\ell_p(\boldsymbol{\beta}, \boldsymbol{\lambda})\)
evaluated at \(\hat{\boldsymbol{\beta}}_{\boldsymbol{\lambda}}\) and
\(M_p\) is number of zero eigenvalues in
\(\boldsymbol{S}_{\boldsymbol{\lambda}}\).

Let \(\hat{\boldsymbol{\lambda}}\) denote the value of
\(\boldsymbol{\lambda}\) that maximises \(\ell(\boldsymbol{\lambda})\)
w.r.t. \(\boldsymbol{\lambda}\). Obtaining
\(\hat{\boldsymbol{\lambda}}\) is an iterative procedure in which each
evaluation of \(\ell(\boldsymbol{\lambda})\) involves obtaining
\(\hat{\boldsymbol{\beta}}_{\boldsymbol{\lambda}}\), which is achieved
in deform using Newton's method with first and second derivatives
of \(\ell_p(\boldsymbol{\beta}, \boldsymbol{\lambda})\) w.r.t.
\(\hat{\boldsymbol{\beta}}_{\boldsymbol{\lambda}}\) calculated
analytically. Higher-order derivatives in deform are approximated
by finite differencing due to their analytical complexity, so that
quasi-Newton methods are used to estimate \(\boldsymbol{\lambda}\).

Uncertainty in \(\boldsymbol{g}\) is quantified in <cit.>
using samples of parameters from their posterior distributions. Here,
however, its uncertainty is readily calculated once
\(\hat{\boldsymbol{\lambda}}\) has been obtained by assuming that the
sampling distribution of
\(\hat{\boldsymbol{\beta}}_{\hat{\boldsymbol{\lambda}}}\) is
\(MVN(\hat{\boldsymbol{\beta}}_{\hat{\boldsymbol{\lambda}}}, {\bf H}(\hat{\boldsymbol{\beta}}_{\hat{\boldsymbol{\lambda}}})^{-1})\);
see <ref>.

§ FUNCTIONS AND MODEL SPECIFICATIONS

The package deform mainly relies on the functions aniso(),
deform() and expand(). These all fit zero-mean GPs.
Specifically, aniso() fits a conventional anisotropic model,
i.e. where \((x_1^*, x_2^*) = (\alpha_1 x_1, \alpha_2 x_2)\), with
\(\alpha_1, \alpha_2 > 0\). Then deform() and expand() fit
the spatial deformation and dimension expansion models of
<ref> and <ref>, respectively.

Their core usage is given by

aniso(x, z, n)
deform(x, z, n, k)
expand(x, z, n, k)

where x is an \(m \times 2\) matrix of coordinates, with each row
giving longitude and then latitude, z is an \(m \times m\)
empirically-calculated variance-covariance matrix, n \(= n\) is
the number of realizations from which z has been calculated and
k is a vector that specifies the ranks of latent dimensions.
Alternatively, x can be supplied as a list with elements
x, z and n, as described above.

Each model fitting function by default fits the powered exponential
covariance function, so that
\(v(\boldsymbol{g}(\boldsymbol{x}), \boldsymbol{g}(\boldsymbol{x}')) = \sigma^2 \rho(||\boldsymbol{g}(\boldsymbol{x}) - \boldsymbol{g}(\boldsymbol{x}')||)\),
\[
\rho(l) =
\left\{\begin{array}{ll}
(1 - \kappa) \exp(-l^\gamma) & \text{if} ~~ l \neq 0,\\
1 & \text{if} ~~ l = 0,
\end{array}\right.
\]
for \(0 \leq \kappa < 1\) and \(0 < \gamma \leq 2\), and
\(\norm{\cdot}\) denotes Euclidean distance. The powered exponential
form is chosen for its greater flexibility over the exponential form and
analytical tractability over the Matérn form. Specifying
correlation = TRUE fixes \(\sigma^2 = 1\) and specifying
cosine = TRUE takes
\(v(\boldsymbol{g}(\boldsymbol{x}), \boldsymbol{g}(\boldsymbol{x}')) = \sigma^2 \rho(||\boldsymbol{g}(\boldsymbol{x}) - \boldsymbol{g}(\boldsymbol{x}')||) \cos(||\boldsymbol{g}(\boldsymbol{x}) - \boldsymbol{g}(\boldsymbol{x}')|| / \phi)\),
for \(\phi > 0\).

The default values in deform() are k = c(10, 10), so that
\(g_1(\boldsymbol{x})\) and \(g_2(\boldsymbol{x})\) are both represented
as rank-10 thin plate regression splines. These ranks can be changed.
The default value in expand() is k = 10, so that \(r = 1\)
and \(g_3(\boldsymbol{x})\) is also represented as a rank-10 thin plate
regression spline. Various illustrations below show aniso(),
deform() and expand() in action.

§ ILLUSTRATIONS

§.§ Model fitting and visualization

This section demonstrates the methods introduced in <ref> and
<ref> on solar radiation data for British Columbia. These data
were used in <cit.>'s original paper on spatial
deformations, and originated from <cit.>. The radiation data
serve as proof-of-concept data, due to their popularization in
subsequent related works, such as <cit.> and
<cit.>. Similarly to <cit.>, here the \(n=732\)
spring-summer measurements (22 March 1980 – 20 September 1983,
excluding 21 September to 21 March each year) on solar radiation from 12
monitoring stations are studied. These data are available in
deform as dataset solar. The following accesses the data

R> library(deform)
R> data(solar)

and then the following shows the locations of the 12 monitoring

R> library(mapdata)
R> maps::map('worldHires',
+           xlim = range(pretty(solar$x[, 1])),
+           ylim = range(pretty(solar$x[, 2])))
R> box()
R> maps::map.cities(canada.cities, minpop = 2e4)
R> points(solar$x, pch = 19)
\end{CodeInput}
\begin{figure}

{\centering \includegraphics[width=0.7\linewidth]{deform_files/figure-latex/solar-1}

\caption[Locations of the 12 monitoring stations]{Locations of the 12 monitoring stations.}\label{fig:solar}
\end{figure}
\end{CodeChunk}

\subsubsection{Anisotropic modeling} \label{S:solar:aniso}

To start, we'll fit the conventional anisotropic model. This serves as a
benchmark against which subsequent nonstationary covariance models will
be compared. The dataset \code{solar} is a \code{list} in the form
described in \S\ref{deform:fct}, so the anisotropic model is fit simply
with the following.

\begin{CodeChunk}
\begin{CodeInput}
R> m0 <- aniso(solar)
\end{CodeInput}
\end{CodeChunk}

Once a model has been fit, \code{plot()} can be used to see a
representation of the model.

\begin{CodeChunk}
\begin{CodeInput}
R> plot(m0, asp = 1)
\end{CodeInput}
\begin{figure}

{\centering \includegraphics[width=0.7\linewidth]{deform_files/figure-latex/aniso2-1}

\caption[Anisotropic model representation]{Anisotropic model representation.}\label{fig:aniso2}
\end{figure}
\end{CodeChunk}

This takes a regular grid and then shows how this changes as a result of
the fitted model. In the case of the conventional anisotropic model,
this is rather trivial, as coordinates are just scaled.

\subsubsection{Spatial deformation} \label{S:solar:deform}

The following code then fits the default spatial deformation model.

\begin{CodeChunk}
\begin{CodeInput}
R> m1 <- deform(solar)
\end{CodeInput}
\end{CodeChunk}

We can supply \code{plot()} with longitude and latitude coordinates that
define the plotting grid through \code{xp} and \code{yp}, respectively.

\begin{CodeChunk}
\begin{CodeInput}
R> x_plot <- seq(-123.3, -122.25, by = .05)
R> y_plot <- seq(49, 49.4, by = .05)
R> plot(m1, xp = x_plot, yp = y_plot, asp = 1)
\end{CodeInput}
\begin{figure}

{\centering \includegraphics[width=0.7\linewidth]{deform_files/figure-latex/deform2-1}

\caption[Deformation model representation]{Deformation model representation.}\label{fig:deform2}
\end{figure}
\end{CodeChunk}

The resulting representation of the fitted deformation shows a space
very different from that achieved by the simple scaling possible with
the conventional anisotropic model. Variograms can be used in
\pkg{deform} to compare models, which are covered in
\S\ref{S:solar:variogram}.

\subsubsection{Bijective spatial deformation} \label{S:solar:biject}

The representation of the deformation shown in Figure \ref{fig:deform2}
is seen to be non-bijective. We can ask \code{deform()} to numerically
eliminate bijective deformations by specifying \code{bijective = TRUE}.
This implements the numerical method described in \S\ref{S:meth:biject}.

\begin{CodeChunk}
\begin{CodeInput}
R> m3 <- deform(solar, bijective = TRUE)
R> plot(m3, xp = x_plot, yp = y_plot, asp = 1)
\end{CodeInput}
\begin{figure}

{\centering \includegraphics[width=0.7\linewidth]{deform_files/figure-latex/bijective1-1}

\caption{\label{fig:biject1}Bijective deformation model representation.}\label{fig:bijective1}
\end{figure}
\end{CodeChunk}

The argument \code{bijective.args} can be supplied to \code{deform()} as
a \code{list} to specify the parameters that control identifying
bijectivity, i.e.~\(\delta\) and \(\epsilon\) in equation
\eqref{foldpen}. Its default values are
\code{list(mult = 1e3, scl = 1, nx = 40, ny = 40)}, which uses a
triangular tiling of dimension \code{nx} \(\times\) \code{ny} to cover
the original data points, where \(\delta =\) \code{mult} and
\(\epsilon =\) \code{scl}
\(\times~ \hat \alpha_1 \hat \alpha_2 l_1 l_2\), with \(\hat \alpha_1\)
and \(\hat \alpha_2\) resulting from conventional anisotropic model
fits, and \(l_1\) and \(l_2\) are the resolutions of the heights and
widths, respectively, of the triangles in the tiling used to identify
bijectivity. Sensitivity to choice of \(\delta\) and \(\epsilon\) is
briefly studied in Appendix \ref{pensens}.

\subsubsection{One-dimensional expansion} \label{S:solar:expand1}

The following fits a one-dimensional expansion model, i.e.~with
\(r = 1\), based on the definition of \S\ref{S:meth:expand}.

\begin{CodeChunk}
\begin{CodeInput}
R> m2 <- expand(solar)
\end{CodeInput}
\end{CodeChunk}

Note that this model has one \emph{latent} dimension, although
\(\boldsymbol{g}(\boldsymbol{x}) = (g_1(\boldsymbol{x}), g_2(\boldsymbol{x}), g_3(\boldsymbol{x}))\),
and so technically uses three dimensions. The following plots the
model's three dimensions.

\begin{CodeChunk}
\begin{CodeInput}
R> par(mfrow = c(1, 3))
R> plot(m2)
\end{CodeInput}
\begin{figure}

{\centering \includegraphics{deform_files/figure-latex/expand2-1}

\caption[One-dimensional expansion model representation]{One-dimensional expansion model representation.}\label{fig:expand2}
\end{figure}
\end{CodeChunk}

\subsubsection{Two-dimensional expansion} \label{S:solar:expand2}

We can extend the number of latent dimensions by specifying a longer
vector for \code{k}. The following specifies two latent dimensions,
i.e.~\(r = 2\), each of which is represented as a rank-10 thin plate
regression spline, and then fits the model.

\begin{CodeChunk}
\begin{CodeInput}
R> m4 <- expand(solar, k = c(10, 10))
\end{CodeInput}
\end{CodeChunk}

When plotting such a model, we may only be interested in visualizing its
latent dimensions. The following starts plotting at the latent
dimensions, i.e.~omits plotting \(g_1(\boldsymbol{x})\) and
\(g_2(\boldsymbol{x})\). We can also use the \pkg{lattice} package to
produce plots (which gives colour scales), and we further specify that
all plots should be on one panel (\code{onepage = TRUE}).

\begin{CodeChunk}
\begin{CodeInput}
R> plot(m4, start = 3, graphics = 'lattice', onepage = TRUE)
\end{CodeInput}
\begin{figure}

{\centering \includegraphics{deform_files/figure-latex/expand4-1}

\caption[Two-dimensional expansion model representation]{Two-dimensional expansion model representation.}\label{fig:expand4}
\end{figure}
\end{CodeChunk}

\subsection{Variograms} \label{S:solar:variogram}

A simple way to assess the fit of a model is through a variogram,
comparing model-based estimates with empirical counterparts.
Specifically these plot semivariance against distance, where distance is
calculated in \(D\)-space. These can be achieved with function
\code{variogam()}, as shown below for the conventional anisotropic,
bijective deformation, and one- and two-dimensional expansion models.

\begin{CodeChunk}
\begin{CodeInput}
R> par(mfrow = c(2, 2))
R> variogram(m0); title('Conventional anisotropic model')
R> variogram(m3); title('Bijective deformation model')
R> variogram(m2); title('One-dimensional expansion model')
R> variogram(m4); title('Two-dimensional expansion model')
\end{CodeInput}
\begin{figure}

{\centering \includegraphics{deform_files/figure-latex/variogram1-1}

\caption[Variograms for conventional anisotropic, bijective deformation, and one- and two-dimensional expansion models]{Variograms for conventional anisotropic, bijective deformation, and one- and two-dimensional expansion models.}\label{fig:variogram1}
\end{figure}
\end{CodeChunk}

The variograms show that the bijective deformation model improves on the
conventional anisotropic model, and so too do the dimension expansion
models. However, there is little to distinguish between the one- and
two-dimensional expansions. Note that variograms for the non-bijective
and bijective deformation models are very similar, and hence only the
latter's variogram is shown.

\subsection{Prediction} \label{S:deform:predict}

Two forms of prediction are possible from objects of class \code{deform}
through \code{predict()}. The default, demonstrated here with the
bijective deformation model, is to give the coordinates in \(D\)-space
as a \code{data.frame}.

\begin{CodeChunk}
\begin{CodeInput}
R> predict(m3)
\end{CodeInput}
\begin{CodeOutput}
[,1]         [,2]
[1,] -0.01085683  0.202837264
[2,] -0.05234367  0.066223514
[3,] -0.05595379  0.038905276
[4,] -0.08727547 -0.023210978
[5,] -0.08045686 -0.082107777
[6,]  0.02509023  0.022495256
[7,]  0.07745855  0.022532967
[8,]  0.09015940 -0.003449057
[9,]  0.07610715 -0.014919774
[10,]  0.03105723 -0.014947634
[11,] -0.06635145 -0.002552110
[12,] -0.10746621 -0.005829157
\end{CodeOutput}
\end{CodeChunk}

The \code{predict()} argument \code{newdata} can be used to get
\(D\)-space coordinates for arbitrary \(G\)-space coordinates. If
argument \code{se.fit = TRUE} is given to \code{predict()}, then a
\code{list} is returned with elements \code{fitted}, which is as above,
and \code{se.fit}, a \code{data.fame} of pointwise standard errors for
estimated \(D\)-space coordinates. For example, the following gives the
predictions plotted in Figure \ref{fig:bijective1} together with
estimates of their standard error, with just the estimates corresponding
to the first six \(G\)-space coordinates shown below.

\begin{CodeChunk}
\begin{CodeInput}
R> grid_df <- as.matrix(expand.grid(x = x_plot, y = y_plot))
R> lapply(predict(m3, newdata = grid_df, se.fit = TRUE), head)
\end{CodeInput}
\begin{CodeOutput}
[,1]        [,2]
[1,] -0.13414235 -0.12866835
[2,] -0.11767010 -0.11295621
[3,] -0.10106640 -0.09763759
[4,] -0.08525571 -0.08502421
[5,] -0.07213916 -0.07950228
[6,] -0.06106455 -0.07919355

[,1]        [,2]
[1,] 0.013533467 0.013215403
[2,] 0.012245433 0.011957164
[3,] 0.010986218 0.010745496
[4,] 0.009807857 0.009726404
[5,] 0.008749314 0.009148672
[6,] 0.007733122 0.008869357
\end{CodeOutput}
\end{CodeChunk}

Alternatively, \code{predict()} can be used to give the estimated
variance-covariance from a model.

\begin{CodeChunk}
\begin{CodeInput}
R> predict(m3, type = 'vcov')
\end{CodeInput}
\begin{CodeOutput}
[,1]     [,2]     [,3]     [,4]     [,5]     [,6]     [,7]     [,8]
[1,] 55.53827 49.54464 48.35936 45.40355 43.10049 47.75753 47.02726 45.78450
[2,] 49.54464 55.53827 54.41427 51.57226 49.18786 51.87670 49.79774 48.85493
[3,] 48.35936 54.41427 55.53827 52.70079 50.38424 52.14364 49.90845 49.13754
[4,] 45.40355 51.57226 52.70079 55.53827 53.13184 50.47817 48.31829 47.98987
[5,] 43.10049 49.18786 50.38424 53.13184 55.53827 49.29108 47.51718 47.58486
[6,] 47.75753 51.87670 52.14364 50.47817 49.29108 55.53827 53.41883 52.68034
[7,] 47.02726 49.79774 49.90845 48.31829 47.51718 53.41883 55.53827 54.36137
[8,] 45.78450 48.85493 49.13754 47.98987 47.58486 52.68034 54.36137 55.53827
[9,] 45.57995 49.14600 49.55185 48.63880 48.34416 52.96570 54.02394 54.77055
[10,] 46.12323 50.69191 51.30008 50.59472 50.09272 54.00642 53.11693 53.09347
[11,] 46.51092 52.67456 53.81226 54.34254 52.22437 51.62408 49.40505 48.94672
[12,] 45.77364 51.79984 52.75733 54.44953 52.21913 49.85905 47.61900 47.16249
[,9]    [,10]    [,11]    [,12]
[1,] 45.57995 46.12323 46.51092 45.77364
[2,] 49.14600 50.69191 52.67456 51.79984
[3,] 49.55185 51.30008 53.81226 52.75733
[4,] 48.63880 50.59472 54.34254 54.44953
[5,] 48.34416 50.09272 52.22437 52.21913
[6,] 52.96570 54.00642 51.62408 49.85905
[7,] 54.02394 53.11693 49.40505 47.61900
[8,] 54.77055 53.09347 48.94672 47.16249
[9,] 55.53827 53.71861 49.53505 47.76149
[10,] 53.71861 55.53827 51.47838 49.71658
[11,] 49.53505 51.47838 55.53827 53.87270
[12,] 47.76149 49.71658 53.87270 55.53827
\end{CodeOutput}
\end{CodeChunk}

\subsection{Simulation}

The \code{simulate()} function can be used to simulate realizations of
zero-mean GPs from a fitted object of class \code{deform}. Put simply,
realizations of nonstationary GPs in \(G\)-space are given by simulating
stationary GPs in \(D\)-space. The following gives one realization of a
nonstationary GP based on the transformation to \(D\)-space of the
coordinates of the 12 monitoring stations given by the bijective
deformation model, which was represented in Figure \ref{fig:biject1}.

\begin{CodeChunk}
\begin{CodeInput}
R> simulate(m3)
\end{CodeInput}
\begin{CodeOutput}
[1,]  3.8218149
[2,] -0.1519009
[3,] -1.6879336
[4,] -3.6426413
[5,] -1.4220332
[6,] -3.1855555
[7,] -4.3454698
[8,] -2.4847252
[9,] -2.7502581
[10,] -2.3198710
[11,] -2.9696018
[12,] -3.2947852
\end{CodeOutput}
\end{CodeChunk}

The argument \code{nsim} specifies the number of realizations to
simulate, which defaults to one, and argument \code{newdata} can be used
similarly to \code{predict()} to specify the locations in \(G\)-space
for which simulations are sought. The following example shows three
realizations simulated for the grid used to represent fitted models

\begin{CodeChunk}
\begin{CodeInput}
R> par(mfrow = c(1, 3))
R> sims <- simulate(m3, nsim = 3, newdata = grid_df)
R> for (i in 1:3) {
+   sim_mat <- matrix(sims[, i], length(x_plot))
+   image(x_plot, y_plot, sim_mat)
+ }
\end{CodeInput}
\begin{figure}

{\centering \includegraphics{deform_files/figure-latex/simulate2-1}

\caption[Three simulated realizations of GPs based on the bijective deformation model]{Three simulated realizations of GPs based on the bijective deformation model.}\label{fig:simulate2}
\end{figure}
\end{CodeChunk}

Note that simulations are based on a Cholesky decomposition with
pivoting of the model-based variance-covariance matrix from the fitted
model. The pivoting allows for rank-deficient GPs to be simulated.
However, this is likely to be very slow and memory-hungry for
high-dimensional simulations.

\section{Censored data}

\subsection{Correlation and covariance matrix estimation}

Sometimes we might want to estimate a GP from censored data. When only
interested in the tail of some data, instead of modeling all data to
estimate a distribution, only the values of exceedances of threshold can
be considered, and the values of non-exceedances treated as censored.
\cite{bor-ct} introduced the multivariate Gaussian tail model in which a
generalized Pareto model represented exceedances of a threshold and then
a probability integral transformation gave Gaussian data. These tails
were then modeled as multivariate Gaussian, which was estimated using a
pairwise censored likelihood approach; see \cite{ledford1997}, for

The \pkg{deform} package includes three functions related to censored

\begin{Code}
fitcenmvn(x, left)
cencor(x, left, scale)
cencov(x, left, scale)
\end{Code}

Function \code{fitcenmvn()} takes \(n \times m\) matrices \code{x} and
\code{left} and estimates a \(m \times m\) variance-covariance matrix by
fitting bivariate Gaussian distributions to each pairwise combination of
columns of \code{x} by assuming that the \((t, i)\)th entry in \code{x}
is left-censored at the \((t, i)\)th entry in \code{left}, for
\(t = 1, \ldots, n\) and \(i = 1, \ldots, m\). Functions \code{cencor()}
and \code{cencov()} return correlation and covariance matrices,
respectively, which are scaled by the mean estimate if
\code{scale = TRUE}; otherwise \code{cencor()} and \code{cencov()}
assume a zero mean.

The functions above operate by first fitting marginal censored Gaussian
distributions, i.e.~assuming for data
\(y_1(\boldsymbol{x}), \ldots, y_n(\boldsymbol{x})\) a likelihood of the
form \[
L(\mu(\boldsymbol{x}), \tau(\boldsymbol{x})) = \prod_{t = 1}^n h(y_t(\boldsymbol{x}); u_t(\boldsymbol{x}), \mu(\boldsymbol{x}), \tau(\boldsymbol{x})),
\] where \[
h(y; u, \mu, \tau) = \left\{\begin{array}{ll} \Phi_1(u; \mu, \tau) & \text{if } y \leq u,\\ \phi_1(y; \mu, \tau) & \text{if } y > u,\end{array}\right.
\] where \(\Phi_1(\, ; \mu, \tau^2)\) and \(\phi_1(\, ; \mu, \tau^2)\)
denote the mean \(\mu\) variance \(\tau^2\) univariate Gaussian
cumulative distribution (cdf) and probability density (pdf) functions,
respectively. The fitted marginal Gaussian distributions are then used
to give
\(w_t(\boldsymbol{x}) = [y_t(\boldsymbol{x}) - \hat \mu(\boldsymbol{x})] / \hat \tau(\boldsymbol{x})\).
Pairwise combinations are modeled as unit bivariate Gaussian, so that
thecorrelation between \(W_t(\boldsymbol{x})\) and
\(W_t(\boldsymbol{x}')\), \(\rho(\boldsymbol{x}, \boldsymbol{x}')\), is
estimated by maximizing the likelihood
\[L\big(\rho(\boldsymbol{x}, \boldsymbol{x}')\big) = \prod_{t=1}^n f\big(w_t(\boldsymbol{x}), w_t(\boldsymbol{x}'); u_t(\boldsymbol{x}), u_t(\boldsymbol{x}'), \rho(\boldsymbol{x}, \boldsymbol{x}')\big)\]
where \(f\big(w, w'; u, u', \rho\big)\) is given by
\begin{equation*}\left\{\begin{array}{lll} \Phi_2\big(u, u'; \rho) & \text{if} & w \leq u, w' \leq u',\\ \phi_1(w'; 0, 1) \Phi_1(u; \rho w', 1 - \rho^2)  & \text{if} & w \leq u, w' > u',\\ \phi_1(w; 0, 1) \Phi_1(u'; \rho w, 1 - \rho^2)  & \text{if} & w > u, w' \leq u',\\ \phi_2\big(w, w'; \rho) & \text{if} & w > u, w' > u',\end{array}\right. \end{equation*}
and \(\phi_2(\, , \,; \rho)\) and \(\Phi_2(\, , \,; \rho)\) denote the
bivariate unit Gaussian cdf and pdf, respectively. The bivariate
Gaussian cdf is approximated using an \pkg{Rcpp} implementation of Alan
Gentz's MVTDST set of Fortran subroutines

\subsection{A simple stochastic rainfall model}

Here we demonstrate a simple use of the left-censored multivariate
Gaussian distribution by producing a basic stochastic weather generator
of daily rainfall. This uses the rainfall data of \cite{cooley2007} by
modeling the tail of rainfall measurements from 1st April to 31st
October in the years 1990 to 2019 over 10mm as multivariate Gaussian.
This model is for demonstration purposes only: a better approach might
be to model all rainfall amounts using a zero-inflated gamma or Weibull
distribution and then use a probability integral transformation so that
the multivariate Gaussian tail model can be used.

The following loads the data (which are given as datasets \code{COprcp}
and \code{COprcp_meta} in package \pkg{evgam}), combines the rainfall
data with its meta data, and then forms a matrix of the data in which
each column corresponds to station and each row to a daily rainfall

\begin{CodeChunk}
\begin{CodeInput}
R> library(evgam)
R> data(COprcp)
R> COprcp <- cbind(COprcp, COprcp_meta[COprcp$meta_row,])
R> COprcp <- subset(COprcp, as.integer(substr(date, 6, 7)) %in% 4:10)
R> tms <- as.factor(COprcp$date)
R> locs <- as.factor(COprcp$meta)
R> COprcp_mat <- matrix(NA, length(levels(tms)), length(levels(locs)))
R> COprcp_mat[cbind(as.integer(tms), as.integer(locs))] <- COprcp$prcp
\end{CodeInput}
\end{CodeChunk}

We then fit a bivariate Gaussian distribution to rainfall measurements
for each pair of stations, treating the measurements as left-censored at

\begin{CodeChunk}
\begin{CodeInput}
R> thresh <- 10
R> thresh_mat <- matrix(thresh, nrow(COprcp_mat), ncol(COprcp_mat))
R> mS <- fitcenmvn(COprcp_mat, thresh_mat)
\end{CodeInput}
\end{CodeChunk}

The pairwise covariance matrix estimate can then act as
\(\boldsymbol{V}\) in equation \eqref{gplik}. This is supplied below to
\code{expand()} and a rank-12 dimension expansion model is then fitted.

\begin{CodeChunk}
\begin{CodeInput}
R> n0 <- round(ncol(COprcp_mat) * mean(!is.na(COprcp_mat)))
R> m0 <- expand(as.matrix(COprcp_meta[, c('lon', 'lat')]), mS[[2]], n0, k = 12)
\end{CodeInput}
\end{CodeChunk}

The following visualises the dimension expansion

\begin{CodeChunk}
\begin{CodeInput}
R> par(mfrow = c(1, 3))
R> plot(m0)
\end{CodeInput}
\begin{figure}

{\centering \includegraphics{deform_files/figure-latex/prcp4-1}

\caption[Dimension expansion model representation for censored data]{Dimension expansion model representation for censored data.}\label{fig:prcp4}
\end{figure}
\end{CodeChunk}

and then we can also view its variogram, which appears to show an
adequate fit.

\begin{CodeChunk}
\begin{CodeInput}
R> variogram(m0)
\end{CodeInput}
\begin{figure}

{\centering \includegraphics[width=0.65\linewidth]{deform_files/figure-latex/prcp5-1}

\caption[Variogram for censored data]{Variogram for censored data.}\label{fig:prcp5}
\end{figure}
\end{CodeChunk}

Finally, we can use the fitted model to simulate rainfall amounts over
10mm. The following simulates ten realizations of daily rainfall at each
station and then shows the amounts for the first six stations, with any
simulated amounts below 10mm shown as \code{NA}.

\begin{CodeChunk}
\begin{CodeInput}
R> set.seed(5)
R> w_sim <- mS[[1]] + simulate(m0, nsim = 12)
R> y_sim <- replace(w_sim, w_sim <= thresh, NA)
R> head(y_sim)
\end{CodeInput}
\begin{CodeOutput}
[,1] [,2] [,3] [,4] [,5] [,6]     [,7]     [,8]     [,9] [,10] [,11] [,12]
[1,]   NA   NA   NA   NA   NA   NA 11.39131 16.85773       NA    NA    NA    NA
[2,]   NA   NA   NA   NA   NA   NA       NA       NA 46.18505    NA    NA    NA
[3,]   NA   NA   NA   NA   NA   NA       NA       NA       NA    NA    NA    NA
[4,]   NA   NA   NA   NA   NA   NA       NA       NA       NA    NA    NA    NA
[5,]   NA   NA   NA   NA   NA   NA       NA       NA       NA    NA    NA    NA
[6,]   NA   NA   NA   NA   NA   NA 11.15420       NA       NA    NA    NA    NA
\end{CodeOutput}
\end{CodeChunk}

The above example is a simple demonstration based on a fixed set of
stations. More generally, the multivariate Gaussian tail model of
\cite{bor-ct} can be coupled with continuous models for threshold
exceedances, such as in \cite{y2019}, to give a continuous probability
integral transformation to unit Gaussian, and then a nonstationary GP
fitted with \pkg{deform}. This would allow simulations for an entire
domain, as opposed to a fixed set of sites.

\section{Summary}

The \proglang{R} package \pkg{deform} is designed to give functions for
nonstationary GP modeling and subsequent visualization of spatial data.
It allows users to fit GPs with a conventional anisotropic stationary
dependence structure, and also with a nonstationary structure, based on
the methods of spatial deformation and dimension expansion. Deformations
and latent dimensions are achieved using thin plate regression splines
to take coordinates in \(G\)-space and give coordinates in \(D\)-space
such that covariance becomes closer to isotropic. Having fitted a model,
users can view spatial deformations and dimension expansions. For
example, points brought closer together in \(D\)-space have stronger
dependence, given their positions in \(G\)-space. Users can also use
variograms to assess the fit a \pkg{deform} model, which allows the
model-based representation of covariance to be compared against its
empirical counterpart. Predictions of points in \(D\)-space for
arbitrary points in \(G\)-space are possible, and simulations of
realizations of nonstationary GPs are also possible, given arbitrary
points in \(G\)-space.

The models in \pkg{deform} are likely to be very time-consuming to fit
for high dimensional data due to various \(O(m^3)\) calculations, such
as Cholesky decompositions. It is hoped that updates can be brought to
\pkg{deform} that facilitate the modeling of such data, in particular
gridded climate model output. Facilitating high-dimensional simulations
is likely to be implemented prior to high-dimensional model fitting. The
\pkg{deform} package has been introduced by considering
\(\boldsymbol{x}\) defined in terms of longitude and latitude. Although
\(G\)-space must be two-dimensional, it does not need to be defined
geographically: see \cite{cooley2007} for the notion of `climate space'.
Future developments may allow the possibility for latent dimensions in
the dimension expansion approach to be defined in terms of covariates.
For example, adding a latent dimension that is a function of elevation
may help bring isotropy for some environmental phenomena.

Uncertainty estimation in \(\boldsymbol{g}\) currently does not allow
for uncertainty in \(\boldsymbol{\lambda}\). This could be achieved
using the method of \citet[\S 6.11.1]{wood-book}, or the more general
method of \cite{rue2009}. Finite-differencing can be used for some of
the more complex derivatives.

\section*{Computational details}

The results in this paper were obtained using \proglang{R} 4.3.1 with
the \pkg{deform} 1.0.1 package. \proglang{R} itself and \pkg{deform} are
available from the Comprehensive \proglang{R} Archive Network (CRAN) at
\url{https://CRAN.R-project.org/}.

\bibliography{deform.bib}

\begin{appendix}

\section{Sensitivity to parameters that control bijectivity} \label{pensens}

The following shows the effects of changing $\delta$ and $\epsilon$, which control the folding penalty, and are defined in equation \eqref{foldpen}. Changes by one order of magnitude in either direction are considered.

\begin{CodeChunk}
\begin{CodeInput}
R> bij_pars <- cbind(mult = c(1e2, 1e4, 1e3, 1e3), scl = c(1, 1, .1, 10))
R> par(mfrow = c(2, 2), mar = c(2, 2, 3, 1))
R> for (i in 1:4) {
+   args_i <- list(mult = bij_pars[i, 1], scl = bij_pars[i, 2])
+   mod_i <- deform(solar, bijective = TRUE, bijective.args = args_i)
+   plot(mod_i, xp = x_plot, yp = y_plot, asp = 1)
+   title(paste('mult =', args_i$mult, ', ', 'scl =', args_i$scl))
+ }
\end{CodeInput}
\begin{figure}

{\centering \includegraphics[width=0.8\linewidth]{deform_files/figure-latex/bijective2-1}

\caption{Sensitivity to changes in bijectivity arguments \code{mult} and \code{scl}, wihch are related to $\delta$ and $\epsilon$, as detailed in \S\ref{S:solar:biject}.}\label{fig:bijective2}
\end{figure}
\end{CodeChunk}

In summary, too small $\delta$ or too high $\epsilon$ (which shrinks areas) fails to impose enough of a penalty to ensure bijectivity. Therefore if \code{deform()}'s default values fail to bring bijectivity, $\delta$ should be increased and/or $\epsilon$ should be decreased. Users are advised to try increasing $\delta$ first.

\end{appendix}

\end{document}