Rotate a xyz pointcloud to fit xy plane

Hi!

I am trying to calculate a regression plane for a given pointcloud and rotate the pointcloud (and regression plane) in such a way that it aligns with the xy plane. I managed to use the lm() function to fit a linear model through the pointcloud, but I am stuck at using the output (the coefficients?) to calculate the rotations necessary to align the pointcloud with the xy plane.

I looked up rotation matrices (which I think i need), but the math behind going from the coefficients calculated by lm() to such a matrix is beyond my abilities. Can anyone help me on my way? Or is there an easier way to align a pointcloud with the xy plane?

Thanks in advance!

This is the code with a simplified dataset:

library(stats)
library(rgl)

test_pc <- data.frame(x = c(-13.237859,-24.852585,-15.430247,-1.503417,-11.97386,-16.003471,-24.604477,-22.399782,-6.105618,-20.054321,-23.751991,-19.948997,-23.001833,-7.849604,-18.049824,-3.800245,-15.693258,-18.56345),
y = c(-29.660378,-9.196059,-16.974348,-34.623558,-25.25021,-21.308338,-14.051781,-27.175217,-35.1679,-14.564184,-23.321609,-23.961811,-19.686213,-26.698538,-30.610306,-31.007727,-25.055323,-18.545834),
z = c(-31.4743,-33.636326,-33.938457,-33.12648,-32.718159,-32.621037,-32.157677,-29.278849,-31.54184,-33.30085,-29.605991,-30.708748,-30.841011,-33.205978,-30.3218,-33.0793,-31.904222,-32.556953))

test_rp = lm(test_pc$z ~ test_pc$x + test_pc$y)

plot3d(test_pc$x,test_pc$y,test_pc$z, type="s", rad=0.1)
planes3d(test_rp$coef[2], test_rp$coef[3], -1, test_rp$coef[1], col = 4, alpha = 0.3)

This is the plot with the simplified data and regression plane that I would like to orient in such a way that it runs parallel with the xy plane:

pc