diff --git a/R/pcm.R b/R/pcm.R
index 88674cc..356b2a6 100644
--- a/R/pcm.R
+++ b/R/pcm.R
@@ -99,10 +99,14 @@ pcm <- function(df=NULL,items=NULL,grp=NULL,dif.items=NULL,type.dif=NULL,weights
       mod <- olmm(resp ~ 0 + ce(item) + re(0|id),data=df.long,family = adjacent(link = "logit"),weights = df.long$weights,na.action = na.action)
     }
     comod <- coef(mod)
+    comod.se <- sqrt(diag(vcov(mod)))
     # output results
     restab <- t(sapply(1:nbitems,function(x) comod[seq(x,length(comod)-1,nbitems)]))
+    restab.se <- t(sapply(1:nbitems,function(x) comod.se[seq(x,length(comod.se)-1,nbitems)]))
     rownames(restab) <- paste0("item",1:nbitems)
     colnames(restab) <- paste0("delta_",1:maxmod)
+    rownames(restab.se) <- paste0("item",1:nbitems)
+    colnames(restab.se) <- paste0("delta_",1:maxmod)
     restab.dif <- NULL
     beta <- NULL
   }
@@ -158,12 +162,15 @@ pcm <- function(df=NULL,items=NULL,grp=NULL,dif.items=NULL,type.dif=NULL,weights
         mod <- olmm(formudif,data=df.long,family = adjacent(link = "logit"),control=olmm_control(fit=fit),weights = df.long$weights,na.action = na.action)
       }
       comod <- coef(mod)
+      comod.se <- sqrt(diag(vcov(mod)))
       # output results
       nbcoef <- nbitems+length(difvar.nonunif)
       if (is.null(weights)) {
         restab <- t(sapply(1:nbcoef,function(x) comod[seq(x,length(comod)-2-length(difvar.unif),nbitems+length(difvar.nonunif))]))
+        restab.se <- t(sapply(1:nbcoef,function(x) comod.se[seq(x,length(comod.se)-2-length(difvar.unif),nbitems+length(difvar.nonunif))]))
       } else {
         restab <- t(sapply(1:nbcoef,function(x) comod[seq(x,length(comod)-2-length(difvar.unif),nbitems+length(difvar.nonunif))]))
+        restab.se <- t(sapply(1:nbcoef,function(x) comod.se[seq(x,length(comod.se)-2-length(difvar.unif),nbitems+length(difvar.nonunif))]))
       }
       difcoef.unif <- NULL
       if (length(difvar.unif)>0) {
@@ -197,6 +204,8 @@ pcm <- function(df=NULL,items=NULL,grp=NULL,dif.items=NULL,type.dif=NULL,weights
       restab <- restab[1:nbitems,]
       rownames(restab) <- items_o
       colnames(restab) <- paste0("delta_",1:maxmod)
+      rownames(restab.se) <- items_o
+      colnames(restab.se) <- paste0("delta_",1:maxmod)
       restab.dif <- rbind(difcoef.nonunif,difcoef.unif)
       restab.diftype <- matrix(ifelse(type.dif==1,"HOMOGENEOUS","NON-HOMOGENEOUS"))
       restab.diftype <- noquote(restab.diftype)
@@ -250,14 +259,19 @@ pcm <- function(df=NULL,items=NULL,grp=NULL,dif.items=NULL,type.dif=NULL,weights
         mod <- olmm(resp ~ 0 + ge(grp) + ce(item) + re(0|id),data=df.long,family = adjacent(link = "logit"),control=olmm_control(fit=fit),weights=df.long$weights,na.action = na.action)
       }
       comod <- coef(mod)
+      comod.se <- sqrt(diag(vcov(mod)))
       # output results
       if (is.null(weights)) {
         restab <- t(sapply(1:nbitems,function(x) comod[seq(x,length(comod)-2,nbitems)]))
+        restab.se <- t(sapply(1:nbitems,function(x) comod.se[seq(x,length(comod.se)-2,nbitems)]))
       } else {
         restab <- t(sapply(1:nbitems,function(x) comod[seq(x,length(comod)-2,nbitems)]))
+        restab.se <- t(sapply(1:nbitems,function(x) comod.se[seq(x,length(comod.se)-2,nbitems)]))
       }
       rownames(restab) <- items_o
       colnames(restab) <- paste0("delta_",1:maxmod)
+      rownames(restab.se) <- items_o
+      colnames(restab.se) <- paste0("delta_",1:maxmod)
       restab.dif <- NULL
       beta <- comod[length(comod)-1]
       se.beta <- (confint(mod)["grp",2]-beta)/1.96