Edit: This answer is not correct. It's not clear to me what OP is looking for, but take a look at the discussion at https://chat.stackexchange.com/rooms/107527/discussion-between-jeremy-miles-and-user7148318 before answering.
There is no package that can do this (as far as I know).
The lavaan package for sem does not give standard errors for latent variable estimates (and I don't know you can predict a model with new data).
You might be able to do this with merTools::predictInterval(), however I'm having trouble understanding your data. What identifies the question?
You see to have very large number of questions (~ 16000?) - each of which is not answered by very many people, so your standard errors will be high. Have I understood that correctly?
My guess is that the will look like:
dat <- read.csv('https://raw.githubusercontent.com/izeh/n/master/g.csv', stringsAsFactors = F)
dat$item_id <- factor(item_id)
dat_no_1 <- dat[dat$person_id != 1, ]
dat_no_1 <- dat[dat$person_id > 1, ]
fit1 <- lme4::lmer(item_score ~ as.factor(item_id) +
(1 | person_id),
data = dat_no_1)
newdata = dat_no_1,
n.sims = 999)
I can't test this though, as I run out of memory.
Note that this ignores the non-normality issue.