I'm trying to program a simple Monte Carlo experiment that after generating sample data within each MC rep then runs a GLM within it, accumulating the coefficients from each rep in a matrix. As the GLM command needs a dataframe I convert all my matrix held simulated data into this form and then try and run a for loop (see end of included code below). It does not like the $i notation. Can anyone please come up with a workaround? I had a quick look in the glm code to see if I could figure out how to create a modelmatrix myself and call the glm directly on that - it defeated me (as I'm not the greatest of programmers). If anyone could offer advice on that too I'd be grateful - would save the glm command having to convert the dataframe on each call - which on a large number of MC reps might be worth while.
Any help much appreciated (and apologies for all the code comments - my memories terrible so need to include). Best, Steve using DataFrames, GLM, Distributions reps=100 # number of MC sims n=50 # sample size in each sim X0 = rep(1,n) # the constant nLevelsX1 = 5 # distinct values in continuous X1 - make it a factor of n X1=gl(nLevelsX1,int(n/nLevelsX1)) # gl() is in DataFrames X2 = gl(2,1,n) - 1 # the treatmnt dummy X = [X0 X1 X2] # design matrix T=1.0 # the treatment effect B = [10, 1, T] # Beta vector mu = float(X*B) # varying mu parameter in lognormal distribution sigma = rep(1.0,n) # constant sigma parameter in lognormal distribution Obs = Array(Float64, reps,n) # holds data for analysis - each row new dataset # loop below creates the Y (dependent) data - one row for each MC sim for i = 1:n Obs[:,i]= rand(LogNormal(mu[i],sigma[i]),reps) end # create dataframe where first two cols hold the unvarying X # and the remaining cols hold each Y sim df= convert(DataFrame, [X[:,2:3] Obs']) rename!(df.colindex, [(symbol("x$i")=>symbol("Y$(i-2)")) for i in 3:(n+2)]) names(df) # x1 x2 Y1 Y2 ...Yreps GLMresults = Array(Float64, reps,3) # will hold all coeff #next loop fails - please help - won't accept $i for i = 1:reps GLMresults[i,:] = coef(glm(Y$i ~ x1 + x2,df,Normal(),LogLink())) end