diff --git a/NAMESPACE b/NAMESPACE
index 51e48cb6..af195796 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -124,6 +124,7 @@ S3method(family,rq)
 S3method(getCall,summ)
 S3method(get_formula,brmsfit)
 S3method(get_formula,default)
+S3method(get_formula,panelmodel)
 S3method(glance,summ.glm)
 S3method(glance,summ.merMod)
 S3method(glance,summ.svyglm)
diff --git a/R/make_new_data.R b/R/make_new_data.R
index 9ea04bbc..72922ba9 100644
--- a/R/make_new_data.R
+++ b/R/make_new_data.R
@@ -386,6 +386,14 @@ get_formula.brmsfit <- function(model, resp = NULL, dpar = NULL, ...) {
   }
 }
 
+#' @rdname get_formula
+#' @export
+get_formula.panelmodel <- function(model, ...) {
+  f <- formula(model)
+  class(f) <- class(f) %not% "Formula"
+  f
+}
+
 get_family <- function(model, ...) {
   UseMethod("get_family")
 }
@@ -416,6 +424,11 @@ get_family.brmsfit <- function(model, resp = NULL, ...) {
   }
 }
 
+#' @importFrom stats gaussian
+get_family.plm <- function(model, ...) {
+  gaussian(link = "identity")
+}
+
 # formerly built into make_new_data, but I want to use it for other times
 # when I just want the values of non-focal predictors
 get_control_values <- function(model, data, preds, at, center, design = NULL,
diff --git a/R/make_predictions.R b/R/make_predictions.R
index e998ab99..e115e6f4 100644
--- a/R/make_predictions.R
+++ b/R/make_predictions.R
@@ -152,7 +152,7 @@ make_predictions.default <- function(model, pred, pred.values = NULL, at = NULL,
   } else {pm <- new_data}
   
   
-  link_or_lm <- ifelse(family(model)$link == "identity",
+  link_or_lm <- ifelse(get_family(model)$link == "identity",
                        yes = "response", no = "link")
   
   # Do the predictions using built-in prediction method if robust is FALSE
@@ -211,10 +211,10 @@ make_predictions.default <- function(model, pred, pred.values = NULL, at = NULL,
   # Back-convert the predictions to the response scale
   if (outcome.scale == "response") {
     pm[[get_response_name(model)]] <-
-      family(model)$linkinv(pm[[get_response_name(model)]])
+      get_family(model)$linkinv(pm[[get_response_name(model)]])
     if (interval == TRUE) {
-      pm[["ymax"]] <- family(model)$linkinv(pm[["ymax"]])
-      pm[["ymin"]] <- family(model)$linkinv(pm[["ymin"]])
+      pm[["ymax"]] <- get_family(model)$linkinv(pm[["ymax"]])
+      pm[["ymin"]] <- get_family(model)$linkinv(pm[["ymin"]])
     }
   }
   
diff --git a/man/get_formula.Rd b/man/get_formula.Rd
index ce90a3b7..027a3b54 100644
--- a/man/get_formula.Rd
+++ b/man/get_formula.Rd
@@ -4,6 +4,7 @@
 \alias{get_formula}
 \alias{get_formula.default}
 \alias{get_formula.brmsfit}
+\alias{get_formula.panelmodel}
 \title{Retrieve formulas from model objects
 
 This function is primarily an internal helper function in \code{jtools} and
@@ -15,6 +16,8 @@ get_formula(model, ...)
 \method{get_formula}{default}(model, ...)
 
 \method{get_formula}{brmsfit}(model, resp = NULL, dpar = NULL, ...)
+
+\method{get_formula}{panelmodel}(model, ...)
 }
 \arguments{
 \item{model}{The fitted model object.}