# Estimating Variance as a Function of Treatment Rank Class

[This article was first published on

Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.

**Econometrics by Simulation**, and kindly contributed to R-bloggers]. (You can report issue about the content on this page here)Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.

Imagine that we have a treatment that we give to five different groups of individuals. Each individual has a variable response which as a unique mean and variance based on the treatment. We do not know how the means will change but we believe the variance of responses will expand depending upon what level of treatment the individual gets. We would like to expressly model both the differences in means and that of the variances.

This code was formulated in response to a question posted on CrossValidated.

We want to solve

$$ \max_{`\bf {`

\hat\beta,\hat\gamma}} (\sum_{i=1}(ln(D(x_i, \hat\mu, \hat\gamma_0+\hat\gamma_1 rank))) $$

# Specify how many individuals are in each of our groups nobs.group <- 500 # Simulate our data grp1 <- data.frame(values=rnorm(nobs.group,5,1), grp=1) grp2 <- data.frame(values=rnorm(nobs.group,3,2), grp=2) grp3 <- data.frame(values=rnorm(nobs.group,6,3), grp=3) grp4 <- data.frame(values=rnorm(nobs.group,5,4), grp=4) grp5 <- data.frame(values=rnorm(nobs.group,1,5), grp=5) # Group our data into a single object mydata <- rbind(grp1,grp2,grp3,grp4,grp5) # Speficy the function to maximize (minimize) lnp <- function(gamma, x, rank) # I include a negative here because the default option with optim is minimize -sum(log(dnorm(x,gamma[1]*(rank==1)+ gamma[2]*(rank==2)+ gamma[3]*(rank==3)+ gamma[4]*(rank==4)+ gamma[5]*(rank==5), gamma[6]+gamma[7]*rank))) ans <- optim(c( # Specify initial values for parameters to be estimated beta1=1,beta2=1,beta3=1,beta4=1, beta5=1, gamma1=1,gamma2=1), # Specify the function to minimize (maximize) lnp, # Input dependent variable as x and the explanatory variable as rank x=mydata$values, rank=mydata$grp, # Be sure to inlcude the hessian in the return for # calculating standard errors hessian=T) # The standard erros can be estimated using the hessian stand.error <- sqrt(diag(solve(ans$hessian))) # This will create a nice table of results cbind(par.est=ans$par, stand.error, tstat=ans$par/stand.error, pvalue=1-pt(ans$par/stand.error, nrow(mydata)-length(ans$par)), CILower=ans$par+stand.error*qt(.05,nrow(mydata)-length(ans$par)), CIUpper=ans$par+stand.error*qt(.95,nrow(mydata)-length(ans$par)) )

par.est stand.error tstat pvalue CILower CIUpper beta1 4.9894067 0.04038367 123.550112 0.0000000 4.9229567 5.05585658 beta2 2.0009055 0.10955198 18.264440 0.0000000 1.8206415 2.18116942 beta3 3.7640531 0.19407912 19.394427 0.0000000 3.4447027 4.08340355 beta4 6.4818562 0.23879420 27.144111 0.0000000 6.0889287 6.87478375 beta5 -0.5547626 0.29730735 -1.865957 0.9689176 -1.0439715 -0.06555377 gamma1 -0.4849449 0.05684407 -8.531142 1.0000000 -0.5784798 -0.39140993 gamma2 1.3867000 0.04520519 30.675682 0.0000000 1.3123164 1.46108352

To

**leave a comment**for the author, please follow the link and comment on their blog:**Econometrics by Simulation**.R-bloggers.com offers

**daily e-mail updates**about R news and tutorials about learning R and many other topics. Click here if you're looking to post or find an R/data-science job.

Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.