[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