Улучшение производительности модели параметров оценки в R - PullRequest
Новая
       25

Улучшение производительности модели параметров оценки в R

0 голосов
/ 15 апреля 2019

Я разработал подход к корреляции угроз в смешанную модель между случайными эффектами и переменными.Я делаю процесс с nlmixed под SAS.Мой набор данных содержит 9 переменных для 8834 строк.

У меня есть 6 параметров для оценки. 

proc nlmixed data=naisga tech=NMSIMP seed=101;
parms beta0=-2.40 beta1=-0.28 delta0=0 delta1=0 sx=1 sb=1;
eta= beta0 + beta1*mat_age_dec + b;
expeta = exp(eta);
p = expeta/(1+expeta);
py = log((p**petitpoids) * ((1-p)**(1-petitpoids)));
sqrt2pi = sqrt(2*constant('pi'));
px = -log(sx) - log(sqrt2pi*mat_age_dec) - (log(mat_age_dec)-delta0 + delta1*b)**2  / (2*sx**2);
LL = px + py;
model petitpoids~general(LL);
random b ~ normal(0,sb) subject=matid;
run;

И у меня есть этот результат для моих фиксированных параметров. 

beta0   -3.1908
beta1   0.4374

С очень хорошим временем вычислений.

Итак, я хотел бы реализовать это в R, и я сделал следующий код.Я попытался воспроизвести методологию SAS в nlmixed proc.

library(fastGHQuad)
library(lme4)

set.seed(101)

data.naisga = read.csv2(file.choose(), sep=",", header = TRUE)

#Creation of my variable to predict
#Here is my Y
data.naisga$poidsnais.dummy = data.naisga$poidsnais<= 2500

#Mother age in decades
#Here is my X
data.naisga$mat_age_dec = data.naisga$mat_age / 10

#Deletion of missing rows for my predictor
data.naisga2 <- data.naisga[!is.na(data.naisga$mat_age),]

Мои данные выглядят так.

matid poidsnais dprenat mat_age sexe matrace fumeuse poidsnais.dummy mat_age_dec  
1    39      3720      NA      15    2       2      NA   FALSE         1.5
2    39      3260       4      17    1       2      NA   FALSE         1.7
3    39      3910       0      19    1       2      NA   FALSE         1.9
4    39      3320       0      24    1       2       2   FALSE         2.4
5    39      2480       6      25    1       2       2   TRUE         2.5
6    62      2381       5      17    2       2      NA   TRUE         1.7

init.par <- c(-2.40,-0.28,0,0,1,1)

#Function with join density
YX.terme = function(par,y,x,b)
{
  #Parameters to estimate
  beta = par[1:2]
  #Parameters of conditional density of X by b 
  delta = par[3:4]
  logsigmax = par[5]
  #Sigma of b
  logsigmab = par[6]

  py = ifelse(y==1,exp(cbind(1,x)%*%beta),1)/(1+exp(cbind(1,x)%*%beta))
  px = exp((-(x-delta[1]-delta[2]*b)^2)/(2*exp(logsigmax)))/exp(logsigmax/2)
  dG = exp(-b^2/(2*exp(logsigmab)))/exp(logsigmab/2)

  prod(py*(px+log(1.000001))*dG)

}

YX.vec = Vectorize(YX.terme,"b")

rule100 <- gaussHermiteData(100)

#Function which compute my likelihood
YX.mixte = function(par,y,x,id)
{
  uid = unique(id)
  ll = 0

  ll = sapply(uid, FUN = function(i) ll + log(ghQuad(YX.vec,rule100,par=par,y=y[id==i],x=x[id==i])))
  return(-sum(ll))
}

#optim to estimate my parameters 

mat.age.fit = optim(init.par,fn=YX.mixte,method="BFGS",y=data.naisga2$poidsnais.dummy,x=data.naisga2$mat_age_dec,id=data.naisga2$matid,control=list(trace=10,reltol=1e-2))

Эта оптимизированная часть чрезвычайно медленная, поэтому я снизил параметр reltol до 1e-2.Даже если я использую CG Solver, который работает быстрее, Optim работает почти целый день, чтобы вычислить мои параметры с крайне нестабильными результатами.

Итак, я хотел бы знать, в порядке ли мой настоящий код или его можно улучшить различными способами, чтобы повысить его производительность.Если вы знаете несколько хороших трюков, чтобы уменьшить время вычислений для моей части оптимизации, я буду благодарен.

...