[R] Problems with boot and optim
Angelo Canty
canty at math.mcmaster.ca
Thu Sep 23 22:25:08 CEST 2004
Hi Luke,
I think your problem is in the function lik.hetprobit and not
remembering that R is case sensitive so X and x are not the same.
The parameters passed in are called X, Y and Z which change for each
bootstrap dataset. Within the function, however, your first three
lines are
Y <- as.matrix(y)
X <- as.matrix(x)
Z <- as.matrix(z)
since x, y, z (lowercase) do not exist in the function, they are being
sought in the global workspace which remains the same for each bootstrap
dataset so that after these three lines your X, Y and Z (uppercase) take
on these values no matter what was input to the function.
Replace x, y and z in the function by X, Y and Z and it should work.
HTH, Angelo
On Tue, 21 Sep 2004, Luke Keele wrote:
>
> I am trying to bootstrap the parameters for a model that is estimated
> through the optim() function and find that when I make the call to boot,
> it runs but returns the exact same estimate for all of the bootstrap
> estimates. I managed to replicate the same problem using a glm() model
> but was able to fix it when I made a call to the variables as data frame
> by their exact names. But no matter how I refer to the variables in the
> het.fit function (see below) I get the same result. I could bootstrap
> it with the sample command and a loop, but then the analysis in the next
> step isn't as nice
>
> The code for the likelihood and the call to boot is below. I have tried
> numerous other permutations as well.
>
> I am using R 1.9.1 on Windows XP pro.
>
> Thanks
>
> Luke Keele
>
>
> #Define Likelihood
> lik.hetprobit <-function(par, X, Y, Z){
>
> #Pull Out Parameters
> Y <- as.matrix(y)
> X <- as.matrix(x)
> Z <- as.matrix(z)
> K <-ncol(X)
> M <-ncol(Z)
> b <- as.matrix(par[1:K])
> gamma <- as.matrix(par[K+1:M])
>
> mu <- (X%*%b)
> sd <- exp(Z%*%gamma)
>
> mu.sd <-(mu/sd)
>
> #Form Likelihood
>
> log.phi <- pnorm(ifelse(Y == 0, -1, 1) * mu.sd, log.p = TRUE)
> 2 * sum(log.phi)
> }
>
> y <- as.matrix(abhlth)
> x <- as.matrix(reliten)
> ones = rep(1, nrow(x))
> x = cbind(ones,x)
> z = as.matrix(abinfo)
>
> data.het <- as.matrix(cbind(y,x,z))
>
> het.fit <- function(data){
> mod <- optim(c(1,0,0), lik.hetprobit, Y=data$y, X=data$x, Z=data$z,
> method="BFGS",
> control=list(fnscale=-1), hessian=T)
> c(mod$par)
> }
> case.fun <- function(d,i)
> het.fit(d[i,])
>
> het.case <- boot(data.het, case.fun, R=50)
> Luke Keele
> Post-Doctoral Fellow in Quantitative Methods
> Nuffield College, Oxford University
> Oxford, UK
>
>
> [[alternative HTML version deleted]]
>
> ______________________________________________
> R-help at stat.math.ethz.ch mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
>
--
------------------------------------------------------------------
| Angelo J. Canty Email: cantya at mcmaster.ca |
| Mathematics and Statistics Phone: (905) 525-9140 x 27079 |
| McMaster University Fax : (905) 522-0935 |
| 1280 Main St. W. |
| Hamilton ON L8S 4K1 |
More information about the R-help
mailing list