diff --git a/DESCRIPTION b/DESCRIPTION
index 0704598..f333629 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -1,11 +1,11 @@
-Package: flipscores
-Version: 1.3.0
-Date: 2023-01-31
-Title: Robust Score Testing in GLMs, by Sign-Flip Contributions 
-Author: Livio Finos, Jelle Goeman and Jesse Hemerik, with contribution of Riccardo De Santis.
-Maintainer: Livio Finos <livio.finos@unipd.it>
-Description: Provides robust tests for testing in GLMs, by sign-flipping score contributions. The tests are robust against overdispersion, heteroscedasticity and, in some cases, ignored nuisance variables. See Hemerik, Goeman and Finos (2020) <doi:10.1111/rssb.12369>.
-Imports: MASS, methods 
-Encoding: UTF-8
-License: GPL-2
-RoxygenNote: 7.1.2
+Package: flipscores
+Version: 1.3.0
+Date: 2023-01-31
+Title: Robust Score Testing in GLMs, by Sign-Flip Contributions 
+Author: Livio Finos, Jelle Goeman and Jesse Hemerik, with contribution of Riccardo De Santis.
+Maintainer: Livio Finos <livio.finos@unipd.it>
+Description: Provides robust tests for testing in GLMs, by sign-flipping score contributions. The tests are robust against overdispersion, heteroscedasticity and, in some cases, ignored nuisance variables. See Hemerik, Goeman and Finos (2020) <doi:10.1111/rssb.12369>.
+Imports: MASS, methods 
+Encoding: UTF-8
+License: GPL-2
+RoxygenNote: 7.2.3
diff --git a/R/compute_scores.R b/R/compute_scores.R
index 603d3dc..58d4f73 100644
--- a/R/compute_scores.R
+++ b/R/compute_scores.R
@@ -68,7 +68,7 @@ compute_scores <- function(model0, model1,score_type = "standardized"){
     .get_1score_basic <- function(X){
       B=t(t(X * D_vect) %*% (diag(sqrtinvV_vect**2, nrow = length(sqrtW))))
       scores = B * (sqrtinvV_vect_times_residuals) #* (1/sum(!is.na(model0$y))^0.5)
-      scale_objects=list(nrm=sqrt(sum(B^2)*sum((sqrtinvV_vect_times_residuals)^2)/sum(!is.na(model0$y))))
+      scale_objects=list(nrm=sqrt(sum(B^2)*sum((sqrtinvV_vect_times_residuals)^2)))
       list(scores=scores, scale_objects=scale_objects)
     }
 
@@ -87,7 +87,7 @@ compute_scores <- function(model0, model1,score_type = "standardized"){
         dinv=deco$d^-1
         Xt=(B%*%deco$u%*%diag(deco$d))
         scores=(diag(dinv)%*%t(deco$u)%*%(residuals))[,]#*(1/sum(!is.na(model0$y))**0.5)
-        nrm=sqrt(sum(B^2)*sum(residuals^2)/sum(!is.na(model0$y)))
+        nrm=sqrt(sum(B^2)*sum(residuals^2))
         scale_objects=list(U=deco$u,B=B,nrm=nrm,Xt=Xt)
        
         list(scores=scores, scale_objects=scale_objects)
@@ -114,7 +114,7 @@ compute_scores <- function(model0, model1,score_type = "standardized"){
       .get_1score_effective <- function(X){
         B<-X*(sqrtW)-t(crossprod(crossprod(A,X*(sqrtW)),solve(crossprod(A),t(A))))
         scores=B*sqrtinvV_vect_times_residuals#/(sum(!is.na(model0$y))**0.5)
-        scale_objects=list(nrm=sqrt(sum(B^2)*sum((sqrtinvV_vect_times_residuals)^2)/sum(!is.na(model0$y))))
+        scale_objects=list(nrm=sqrt(sum(B^2)*sum((sqrtinvV_vect_times_residuals)^2)))
         list(scores=scores, scale_objects=scale_objects)
       }
       
@@ -134,7 +134,7 @@ compute_scores <- function(model0, model1,score_type = "standardized"){
           # we divide it by sqrt(m) which is the sd scaling factor of the observed test stat (i.e. effective and standardized have the same observed test stat)
           A=b[,]*U/sqrt(m)
           scores=b*sqrtinvV_vect_times_residuals#/(sum(!is.na(model0$y))**0.5)
-          nrm=sqrt(sum(b^2)*sum((sqrtinvV_vect_times_residuals)^2)/sum(!is.na(model0$y)))
+          nrm=sqrt(sum(b^2)*sum((sqrtinvV_vect_times_residuals)^2))
           scale_objects=list(A=A,nrm=nrm)
           list(scores=scores, scale_objects=scale_objects)
         }
@@ -152,7 +152,7 @@ compute_scores <- function(model0, model1,score_type = "standardized"){
           .get_1score_orthogonalized <- function(X){
             B=(t(X*sqrtW)%*%OneMinusH*(sqrtinvV_vect))
             scores=t(B%*%deco$u)*(t(deco$u)%*%(residuals))[,]#*(1/sum(!is.na(model0$y))**0.5)
-            nrm=sqrt(sum(B^2)*sum(residuals^2)/sum(!is.na(model0$y)))
+            nrm=sqrt(sum(B^2)*sum(residuals^2))
             scale_objects=list(U=deco$u,B=B,nrm=nrm)
             list(scores=scores, scale_objects=scale_objects)
           }
diff --git a/R/flipscores.R b/R/flipscores.R
index 1386456..14209f7 100644
--- a/R/flipscores.R
+++ b/R/flipscores.R
@@ -141,7 +141,7 @@ flipscores<-function(formula, family, data,
     model <- eval(mf, parent.frame())
   } else { # input is a model
     param_x_ORIGINAL <- TRUE
-    model <- update(model,x=TRUE,family=model$family)
+    model <- update(model,x=TRUE)
   }
   
   if(is.null(model$y)) model$y=model$model[,1]
diff --git a/R/utility_functions.R b/R/utility_functions.R
index 7c5095b..0a4c420 100644
--- a/R/utility_functions.R
+++ b/R/utility_functions.R
@@ -53,7 +53,7 @@ summary.flipscores <- function (object, ...) {
   sum_model$coefficients[names(object$p.values),5]=(sum_model$coefficients[names(object$p.values),2]/attributes(object$scores)$nrm)[]
   sum_model$coefficients[names(object$p.values),6]=object$p.values
   # sum_model$coefficients=sum_model$coefficients[,c(1,4)]
-  colnames(sum_model$coefficients)[c(2,4,5)]=c("Score","z value","eff_size")
+  colnames(sum_model$coefficients)[c(2,4,5)]=c("Score","z value","Part. Cor")
   
   sum_model$aliased=rep(FALSE,length(sum_model$aliased))
   structure(sum_model, heading = get_head_flip_out(object), class = c("data.frame"))