# Estimating Variance as a Function of Treatment Rank Class

March 24, 2014
By

(This article was first published on Econometrics by Simulation, and kindly contributed to R-bloggers)

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 groupsnobs.group <- 500 # Simulate our datagrp1 <- 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 objectmydata <- 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 hessianstand.error <- sqrt(diag(solve(ans$hessian))) # This will create a nice table of resultscbind(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     CIUpperbeta1   4.9894067  0.04038367 123.550112 0.0000000  4.9229567  5.05585658beta2   2.0009055  0.10955198  18.264440 0.0000000  1.8206415  2.18116942beta3   3.7640531  0.19407912  19.394427 0.0000000  3.4447027  4.08340355beta4   6.4818562  0.23879420  27.144111 0.0000000  6.0889287  6.87478375beta5  -0.5547626  0.29730735  -1.865957 0.9689176 -1.0439715 -0.06555377gamma1 -0.4849449  0.05684407  -8.531142 1.0000000 -0.5784798 -0.39140993gamma2  1.3867000  0.04520519  30.675682 0.0000000  1.3123164  1.46108352


R-bloggers.com offers daily e-mail updates about R news and tutorials on topics such as: Data science, Big Data, R jobs, visualization (ggplot2, Boxplots, maps, animation), programming (RStudio, Sweave, LaTeX, SQL, Eclipse, git, hadoop, Web Scraping) statistics (regression, PCA, time series, trading) and more...

If you got this far, why not subscribe for updates from the site? Choose your flavor: e-mail, twitter, RSS, or facebook...

## Recent popular posts

Contact us if you wish to help support R-bloggers, and place your banner here.

# Never miss an update! Subscribe to R-bloggers to receive e-mails with the latest R posts.(You will not see this message again.)

Click here to close (This popup will not appear again)