title: predict() - multinomial logistic regression

text: I come to you today because I face a huge problem that I cannot explain. I have run a multinomial logistic regression (using the mlogit package) on behavioral data. I prepare the data by doing                   mlogit <- mlogit.data(Merge, choice = "Choice", shape = "long", alt.var = "Comp",                                drop.index = TRUE)      on my Merge data. which gives me the following:                               Date     Time ActivityX ActivityY Temp Behavior Valley Age Month Year kid Individual Choice         1.F   01/05/2012 00:00:00        80        58   10        F  Fuorn   8     5 2012   Y         26   TRUE         1.R   01/05/2012 00:00:00        80        58   10        F  Fuorn   8     5 2012   Y         26  FALSE         1.M   01/05/2012 00:00:00        80        58   10        F  Fuorn   8     5 2012   Y         26  FALSE         1.RUN 01/05/2012 00:00:00        80        58   10        F  Fuorn   8     5 2012   Y         26  FALSE         2.F   01/05/2012 00:05:00        90        76   10        F  Fuorn   8     5 2012   Y         26   TRUE         2.R   01/05/2012 00:05:00        90        76   10        F  Fuorn   8     5 2012   Y         26  FALSE         2.M   01/05/2012 00:05:00        90        76   10        F  Fuorn   8     5 2012   Y         26  FALSE         2.RUN 01/05/2012 00:05:00        90        76   10        F  Fuorn   8     5 2012   Y         26  FALSE         3.F   01/05/2012 00:10:00        51        47   10        M  Fuorn   8     5 2012   Y         26  FALSE         3.R   01/05/2012 00:10:00        51        47   10        M  Fuorn   8     5 2012   Y         26  FALSE         3.M   01/05/2012 00:10:00        51        47   10        M  Fuorn   8     5 2012   Y         26   TRUE         3.RUN 01/05/2012 00:10:00        51        47   10        M  Fuorn   8     5 2012   Y         26  FALSE         4.F   01/05/2012 00:15:00         0         0   10        R  Fuorn   8     5 2012   Y         26  FALSE         4.R   01/05/2012 00:15:00         0         0   10        R  Fuorn   8     5 2012   Y         26   TRUE         4.M   01/05/2012 00:15:00         0         0   10        R  Fuorn   8     5 2012   Y         26  FALSE         4.RUN 01/05/2012 00:15:00         0         0   10        R  Fuorn   8     5 2012   Y         26  FALSE         5.F   01/05/2012 00:20:00         0         0    9        R  Fuorn   8     5 2012   Y         26  FALSE         5.R   01/05/2012 00:20:00         0         0    9        R  Fuorn   8     5 2012   Y         26   TRUE         5.M   01/05/2012 00:20:00         0         0    9        R  Fuorn   8     5 2012   Y         26  FALSE         5.RUN 01/05/2012 00:20:00         0         0    9        R  Fuorn   8     5 2012   Y         26  FALSE      then I ran my regression :               m1 <- mlogit(Choice ~ 1 |Temp + Valley + Age + kid + Month , mlogit)      and it gave me significant results :                                         Estimate  Std. Error  t-value  Pr(>|t|)             M:(intercept)      -4.2153e-01  5.7533e-02  -7.3268 2.358e-13 ***         R:(intercept)       6.2325e-01  3.4958e-02  17.8284 < 2.2e-16 ***         RUN:(intercept)    -1.2275e+01  4.0526e-01 -30.2895 < 2.2e-16 ***         M:Temp              1.5371e-02  9.8680e-04  15.5764 < 2.2e-16 ***         R:Temp             -3.9871e-02  6.7926e-04 -58.6975 < 2.2e-16 ***         RUN:Temp           -4.4532e-02  6.8696e-03  -6.4825 9.023e-11 ***         M:ValleyTrupchun   -3.6154e-01  1.6362e-02 -22.0968 < 2.2e-16 ***         R:ValleyTrupchun   -4.0186e-02  9.7968e-03  -4.1020 4.096e-05 ***         RUN:ValleyTrupchun  1.2895e+00  8.5357e-02  15.1066 < 2.2e-16 ***         M:Age              -1.1026e-02  2.6902e-03  -4.0985 4.158e-05 ***         R:Age               1.9465e-02  1.6479e-03  11.8119 < 2.2e-16 ***         RUN:Age             5.5473e-02  1.6661e-02   3.3294 0.0008703 ***         M:kidY              6.0686e-02  2.2638e-02   2.6807 0.0073460 **          R:kidY             -4.1638e-01  1.2391e-02 -33.6024 < 2.2e-16 ***         RUN:kidY            6.2311e-01  1.0410e-01   5.9854 2.158e-09 ***         M:Month            -2.0466e-01  8.4448e-03 -24.2346 < 2.2e-16 ***         R:Month             2.4148e-02  5.2317e-03   4.6157 3.917e-06 ***         RUN:Month           9.8715e-01  5.6209e-02  17.5622 < 2.2e-16 ***      those results were in line with what I expected to find in literature so I was quite happy. My next step was to plot my results and here is when I have some trouble. First of all when I plot my original data and compare it with the result of my regression I find some huge differences. For example, when I plot the %of time spend in a behavior (M for moving, F for feeding, R for resting and Run for running, in my regression F is the reference) in function of age, I find that the older an individual gets, the more they will rest and the more they will move, but the estimates I got from my regression shows that they should rest more (when they get older) but move less. So to summarize, my graph on the original data shows the opposite as what I got from the regression. I don't know if it is normal, in the sense that I don't know if I can compare my original data to the result of my regression in a way that my regression shows the probability from switching to a behavior from an other each time my variable grows of one unit. So I wanted to use the `predict()` function but I don't know how to do that. I was hoping to get some help here.


title: Statistical test for computer-systems performance analysis?

text: Suppose we do the following study. We have 2 web-servers (server A and server B) which are processing user requests. Both run on the same hardware, but one of them (server B) uses a slightly different algorithm for request processing. When a user request comes in we have 50/50 chance of processing it on each server. So the work being done by each server is the same in long term. We have 1 hour system-log containing information about all the requests in this hour ( _all the population, not sample_ ): what server processed the request and how long it took to process. For example:               21:00:00 SERVER-A 0.11s     21:00:00 SERVER-A 0.13s     21:00:00 SERVER-B 0.09s     21:00:00 SERVER-A 0.16s     21:00:00 SERVER-B 0.11s     21:00:01 SERVER-A 0.14s      The log is quite large (~250K requests from each server, ~500K in total). It turns out, that server B is slightly faster (mean response time is smaller). The question is which test should I use to show statistical significance in difference of response times of those two servers?


title: Expected number of trials to get 6 linearly independent vectors

text: **The Problem** There's a bag with six distinct tokens in it (e.g. a six-sided die). Each token has an unknown value of some real number. A "trial" consists of drawing five tokens (five die rolls) and being told the sum (to clarify, you do know the identifier for each token, just not its value). What's the expected number of trials to be able to determine the values of all six tokens? **Algebraic representation** There are six unknowns $x_1, x_2, ..., x_6$. Each trial gives you an equation of the form $x_{d_1} + x_{d_2} + x_{d_3} + x_{d_4} + x_{d_5} = r$, where each of the $d_i$s are chosen randomly (independently, each with 1/6 probability) from $\{1,2,3,4,5,6\}$. The logical way to solve it is to convert each trial into a linear equation of the form $$n_{1,i}x_1+n_{2,i}x_2+n_{3,i}x_3+n_{4,i}x_4+n_{5,i}x_5+n_{6,i}x_6=r_i$$ and then consider the resulting system of linear equations. In many (probably most) cases, six trials will suffice. But it's possible to get the same combination of dice more than once (either exactly or as a permutation), causing the six equations to not be linearly independent, which would necessitate another trial. And it's also possible for the latest trial to produce an equation that's a linear combination of some of the previous equations. (e.g. if $x_1+x_1+x_2+x_2+x_3=r_1$ and $x_3+x_4+x_4+x_5+x_5=r_2$ then $x_1+x_2+x_3+x_4+x_5=(r_1+r_2)/2$) So the problem becomes, how many trials are necessary to get 6 linearly independent equations? **My attempts** I thought I'd solve it by building it up one trial at a time, considering the probability that the next trial will produce a linearly independent equation. Then I can find the expected number of trials to get to 2, then expected number of trials to get to 3, etc. However I didn't know how to proceed. I thought I'd consider things casewise the cases where the first trial has no repeats / one pair / two pair / one triple / one four-of-a-kind / two pair / a pair and a triple, but that quickly becomes unfeasible.


title: How do factor categories with no variance influence logistic regression?

text: I'm modeling the effect of a categorical predictor on a binary dependent variable using logistic regression. I'm comparing models with/without the predictor using a likelihood-ratio test. Two categories of the predictor are associated with values of 1 only (no 0s) for the dependent variable. Regression coefficients for these categories (expressed as changes in log(odds) compared to a reference category) are very large and highly suspicious, as this reference category is always associated with response values of 1 (but for one case), and I would thus expect regression coefficients close to 0 for these two categories. Comparisons between the reference category and other categories having more balanced distribution of 1 and 0s matches what I'm expecting from visual inspection of the data. Removing cases associated with these two 'problematic' categories does not change the logLikelihood of the models, but because it changes the number of parameters it affects the results of the likelihood ratio test. Models are fitted using the glm function with binomial family and logit link in R. My question therefore is: what model (or procedure) should I use to: (1) test the global significance of the effect of the predictor on the dependent variable? Should I keep data from the 'problematic' categories in the model or not before conducted the likelihood ratio test? (2) compare these two 'problematic' categories with others? Any hint appreciated,


title: Evaluate $\lim_{n \to \infty} \sum_{j=0}^{n}{{j+n-1} \choose j}\frac{1}{2^{j+n}}$

text: > Using the Central Limit Theorem , Evaluate $$\lim_{n \to \infty} > \sum_{j=0}^{n}{{j+n-1} \choose j}\frac{1}{2^{j+n}}$$ My solution: Let $\{X_n\}$ be a sequence of iid R.V's each having $Geo(\frac{1}{2})$ Then $E(X_n)=1$ and $Var(X_n)=2 < \infty,n\in N$ Also $S_n=\sum_{k=1}^nX_k\sim NB(n,\frac{1}{2})$ By Lindeberg-Levy CLT, $$\begin{align} \lim_{n \to \infty}P[\frac{S_n-E(S_n)}{\sqrt{V(S_n)}}\leq x] &= \Phi(x), \forall x \in R\\ \implies \lim_{n \to \infty}P[\frac{S_n-n}{\sqrt{2n)}}\leq x] &= \Phi(x), \forall x \in R\end{align}$$ Now put $x=0$ $$\begin{align} \lim_{n \to \infty}P[S_n \leq n]&= \Phi(0) \\ \implies \lim_{n \to \infty} \sum_{j=0}^{n}{{j+n-1} \choose j}\frac{1}{2^{j+n}} &=\frac{1}{2}\end{align}$$ Am I right? Is there any other way to solve the limit? So that I can check the limit(Though it is completely unnecessary for the given question). Please help. Thank in advance.


title: What does "20/ln(2)" mean in logistic regression?

text: I am trying to understand Logistic Regression in relation to credit scoring model. I wish to understand the significance of "20/ln(2)" in logistic regression. Why and how is it used?


title: How to test the effect of various manipulations of an independent variable on a dependent variable in a timeseries?

text: I am currently monitoring the daily energy consumption of a house since January of this year. That daily consumption depends on several variables, including outdoor temperature, solar irradiance etc. At regular intervals (once every two weeks) I am changing the algorithm of the heating controller, to see if I can improve the energy savings. What is the methodology one should follow for testing whether such a change has any effect on a timeserie? I've found a short example in BHH2 (where they analyze the effect of a law on pollution levels) but I need a more thorough discussion. **Update** Here is a plot of the data I'm working with: ![Energy vs time](http://farm4.static.flickr.com/3301/5694800429_70e184727c_d.jpg) This shows the daily energy consumption. The green line indicates when we introduced the change to the algorithm. As you can see, there already was a trend towards lower energy consumption even before we introduced the change, probably due to outdoor temperatures getting warmer. What I need to know is whether the introduction of the change has significantly accentuated that trend or not. Note also that since that time, we have collected much more data than what is shown on this plot.


title: Estimating Lambda for Box Cox transformation for ANOVA

text: ### Assumptions: In an ANOVA where the normality assumptions are violated, the Box-Cox transformation can be applied to the response variable. The `lambda` can be estimated by the using maximum likelihood to optimize the normality of the model residuals. ### Question: When the estimates for `lambda` in the null model and the full model differ, how should `lambda` be estimated? ### My Data: In my data the lambda estimate for the null model is `-2.3` and the lambda estimate for the full model is `-2.8`. Transforming the response using these different parameters and preforming the ANOVA leads to different F-statistics. I have produced below a simplified version of the analysis using beta distributions with different parameters to simulate non-normality. Unfortunately, in this example the results of the ANOVA are insensitive to the different estimates of `lambda`. So, it doesn't fully illustrate the problem.               library(ggplot2)     library(MASS)     library(car)               #Generating random beta-distributed data     n=200     df <- rbind(       data.frame(x=factor(rep("a1",n)), y=rbeta(n,2,5)), # more left skewed       data.frame(x=factor(rep("a2",n)), y=rbeta(n,2,2))) # less left skewed          print(qplot(data=df, color=x, x=y, geom="density"))          print("Untransformed Analaysis of Variance:")     m.null <- lm(y ~ 1, df)     m.full <- lm(y ~ x, df)     print(anova(m.null, m.full))          # Estimate Maximum Liklihood Box-Cox transform parameters for both models     bc.null <- boxcox(m.null); bc.null.opt <- bc.null$x[which.max(bc.null$y)]     bc.full <- boxcox(m.full); bc.full.opt <- bc.full$x[which.max(bc.full$y)]          print(paste("ML Box-Cox estimate for null model:",bc.null.opt))     print(paste("ML Box-Cox estimate for full model:",bc.full.opt))          df$y.bc.null <- bcPower(df$y, bc.null.opt)     df$y.bc.full <- bcPower(df$y, bc.full.opt)          print(qplot(data=df, x=x, y=y.bc.null, geom="boxplot"))     print(qplot(data=df, x=x, y=y.bc.full, geom="boxplot"))               print("Analysis of Variance with optimial Box-Cox transform for null model")     m.bc_null.null <- lm(y.bc.null ~ 1, data=df)     m.bc_null.full <- lm(y.bc.null ~ x, data=df)     print(anova(m.bc_null.null, m.bc_null.full))          print("Analysis of Variance with optimial Box-Cox transform for full model")     m.bc_full.null <- lm(y.bc.null ~ 1, data=df)     m.bc_full.full <- lm(y.bc.null ~ x, data=df)     print(anova(m.bc_full.null, m.bc_full.full))


title: How to interpret results of dimensionality reduction/multidimensional scaling?

text: I performed both a SVD decomposition and a multidimensional scaling of a 6-dimensional data matrix, in order to get a better understanding of the structure of the data. Unfortunately, all the singular values are of the same order, implying that the dimensionality of the data is indeed 6. However, I would like to be able to interpret the values of the singular vectors. For instance, the first one seems to be more or less equal in each dimension (ie `(1,1,1,1,1,1)`), and the second also has an interesting structure (something like `(1,-1,1,-1,-1,1)`). How could I interpret these vectors? Could you point me to some literature on the subject?


title: What is more important in choosing parametric or nonparametric tests - skewness & kurtosis or normality test results?

text: There are often differences between results of Skewness & Kurtosis and normality tests, and I have always doubts if it is better to choose parametric or nonparametric tests (I use SPSS). Sometimes histograms show if distribution looks normal or not, and I noticed that most often S&K are better pointers but when I did analysis last time it was different and I really don't know what to do... I read that e.g. when groups are equinumerous, in choosing between t-Student test and nonparametric ones it is better to choose t-Student's even if distributions aren't normal. Is that true?


title: How to decide whether the distribution is unimodal or bimodal in grain size distribution?

text: In the examples given at this link, I am not able to decide whether the distribution is unimodal or bimodal. I think it is in between unimodal and bimodal, but I do not know if this kind of class exists or not. Can anyone please suggest something so that I can easily understand in which case I should say unimodal or bimodal distribution!


title: Plot normal probability for effect estimates in factorial design in R

text: I have 2 level Design (DOE) with 4 factors (A,B,C,D). I've already calculated the estimates for each main effect and all the interaction effects. How can I construct the normal probability plot to see which effects are significant? I have looked at many packages in R for Design of Experiment but cannot find a package that produces the plot. How can I construct the plot with all the estimates I have by code in R Thank you.


title: Control mean lower than treatment mean in a list experiment -- why?

text: I have a list experiment with N=150. The control group was given 3 factors and the treatment group was given 4 factors, with random assignment to these groups. Each group had to report how many factors apply to them. I've done 18 batches of this experiment. In most of them, the mean for the treatment group is higher than the mean for the control group or if it is the other way around, the difference is not significant.BUT in two of the batches, the mean for the treatment group is lower than the mean for the control group and that difference is significant...The questions have been pre-tested and seemed to work well. So why these results in the 2 off-batches? A million thanks in advance for any advice and suggestions, TP


title: Logistic regression: controlling variables not significant, what should I conclude/further test?

text: I ran annual logisitic regression on time-series datas. The most important independant variable have coefficient that are significant in a lot of years, that's a relief. But the "controlling variables", have non-significant coefficients. I'm far from an expert in stats. My sample is very small compared to the litterature that made that test, because I'm analysing a sub-industry that has very few companies in it. But I can't change my sample. In the literature, the authors that made this precise analysis on bigger samples found significant coefficients to controlling variables too, fact they use to state that: "the effect of the main variable is therefore a separate effect than the ones of controlling variables". As I can't say that, should I test the correlation between these variables and the main one, and eliminate some of them that present a very high correlation with it every year, stating that these are redundant variables ? If so, which level is considered high ? 0.75, 0.8 ? I'm planning to say:   * Main variable is significant, but not controlling variable   * There are two potential reasons: 1) My yearly samples are too small 2) The effect of the main variable is the same as the effect of some of these controlling variables   * To try to rule out 2), I tested correlation and found that bla bla bla (if low correlations, probably the samples are just too small, if high, I eliminate some variables and re-run regressions) Thanks ! (non-native english-speaker, sorry) ## ------------------------------------------------------------ Following Peter's answer (thank you !), I think I should add some details about what I'm testing my study: _**Is the propensity to distribute dividends (1 or 0) dependant of the firm's life-cycle stage (using a ratio) ?_** Controlling for profitability, asset growth rate, size (market cap.), and lagged dividends (1 or 0). In every year, the coefficient for lagged dividend is highly significant. And for some of the years (sufficiently for my student thesis), the coefficient of firm's life-cycle stage is significant. The signs of all the other controlling variables are as predicted (most of the time), but always far from significant. I'm hoping that doesn't make my model irrelevant. My pseudo-R square are pretty high. (By the way, does this pseudo-R take in account the variables that are non-significant, or is it the value of the model including only variables that have a significant coefficient ?)


title: Cointegration of a VAR(1) process

text: I am using a Johansen procedure to test for cointegration a vectorial 4-dim vector (timeserie). First I tested for differential stationarity of each individual vector, all of those have a unit root so we are good there. Then I have to find the optimal number of lags as to test for cointegration using Johansen, we assume that the underlying process is a `VAR` process. And that's where the problem occur, I usually find `K=1` (`K` minimizes several information criterias). I use the package `urca` and the function `ca.jo` to run the `trace` and `eigen` tests, but `ca.jo` only allows `K>1` (which makes sense because of the model). What does this `K=1` means ? I add that If I take `K=2` both Johansen tests comes back significant and the process indeed looks cointegrated...


title: Difference and connection between generative learning, discriminative learning and max-margin learning

text: I once heard that, generative learning, discriminative learning and max-margin learning can be separated in terms of their respective definition of loss function. I am not sure how to achieve that?


title: Highly correlated predictors in backward stepwise regression?

text: I know that it's not right to enter variables having multicollinearity (high correlation) into a regression analysis. But if I'm using **backward stepwise regression** could I add all the highly correlated predictors and expect the best predictors to remain in the model at the end or the analysis could go awry? My experience with my own data resulted in removal of the redundant (with multicollinearity) variables thus it seems there's no problem, but I want to report my method and I want to know whether I did it right or it resulted just by chance. So to say briefly: **Is it right to enter variables with multicollinearity into a backward stepwise regression analysis?**


title: Main challenges in data-mining

text: I am an occasional data-mining methods user and I simply wonder what are the main challenges and research directions for method designers. Thanks for your expertise, Peuhp


title: Need help deriving regression statistics

text: I am self-studying this book and need assistance on one of the exercises, specifically 3.9.2: Suppose that, for a certain population, we can predict log earnings from log height as follows:   * A person who is 66 inches tall is predicted to have earnings of $30K.   * Every increase of 1% in height corresponds to a predicted increase of 0.8% in earnings.   * The earnings of approximately 95% of people fall within a factor of 1.1 of predicted values. (a) Give the equation of the regression line and the residual standard deviation of the regression. For this part, I use the first and second bullets to get the equation $\ln(30,000) = .008\ln(66) + c$ and, solving for $c$, I get the regression equation $\ln(earnings) = .008 \ln(height) + 10.257$. Is this correct? For the second part of the question, I need help determining how to find the residual standard deviation of the regression. If $p$ is the value of the prediction line in the residual plot, then it follows that $p \pm 2\sigma = p \pm 1.1p$, based on the third bullet. However, I'm not sure how to actually calculate $\sigma$. (b) Suppose the standard deviation of log heights is 5% in this population. What then, is the $R^2$ of the regression model described here? No idea where to start on this one...could I get a hint, please?


title: Why are rlm() regression coefficient estimates different than lm() in R?

text: I am using rlm in the R MASS package to regress a multivariate linear model. It works well for a number of samples but I am getting quasi-null coefficients for a particular model:               Call: rlm(formula = Y ~ X1 + X2 + X3 + X4, data = mymodel, maxit = 50, na.action = na.omit)     Residuals:            Min         1Q     Median         3Q        Max      -7.981e+01 -6.022e-03 -1.696e-04  8.458e-03  7.706e+01           Coefficients:                  Value    Std. Error t value      (Intercept)    0.0002   0.0001     1.8418     X1             0.0004   0.0000    13.4478     X2            -0.0004   0.0000   -23.1100     X3            -0.0001   0.0002    -0.5511     X4             0.0006   0.0001     8.1489          Residual standard error: 0.01086 on 49052 degrees of freedom       (83 observations deleted due to missingness)      For comparison, these are the coefficients calculated by lm():               Call:     lm(formula = Y ~ X1 + X2 + X3 + X4, data = mymodel, na.action = na.omit)          Residuals:         Min      1Q  Median      3Q     Max      -76.784  -0.459   0.017   0.538  78.665           Coefficients:                   Estimate Std. Error t value Pr(>|t|)         (Intercept)  -0.016633   0.011622  -1.431    0.152         X1            0.046897   0.004172  11.240  < 2e-16 ***     X2           -0.054944   0.002184 -25.155  < 2e-16 ***     X3            0.022627   0.019496   1.161    0.246         X4            0.051336   0.009952   5.159  2.5e-07 ***     ---     Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1           Residual standard error: 2.574 on 49052 degrees of freedom       (83 observations deleted due to missingness)     Multiple R-squared: 0.0182, Adjusted R-squared: 0.01812      F-statistic: 227.3 on 4 and 49052 DF,  p-value: < 2.2e-16       The lm plot doesn't show any particularly high outlier, as measured by Cook's distance: ![lm diagnostic](http://i.stack.imgur.com/98Hxj.png) **EDIT** For reference and after confirming results based on the answer provided by Macro, the R command to set the tuning parameter, `k`, in the Huber estimator is (`k=100` in this case):               rlm(y ~ x, psi = psi.huber, k = 100)


title: Numerical stability of IWLS for Gamma models with log-link

text: The combination of a $\Gamma$-distribution with the log-link function in a generalized linear model can be a useful model. However, in my experience the iterative weighted least squares (IWLS) algorithm (specifically, the implementation in R via `glm`) does not always converge. The log-link is not the canonical link for the $\Gamma$-distribution, and for this reason you cannot use standard results from exponential families to discuss questions about existence and uniqueness of the solution to the score equation. An analysis specifically of the $\Gamma$-distribution with the log- link is fairly straight forward, but this question may have been addressed by others before me. I have just not been able to find anything in the literature. My question is therefore if there are any references on the numerical aspects of generalized linear models with the $\Gamma$-distribution and the log-link? I am specifically interested in discussions related to the IWLS algorithm and whether there are more stable alternatives. **Edit:** A small example. Suppose that $Y$ is log-normally distributed with $\log Y \sim \mathcal{N}(\beta_0 + \beta_1 x, \sigma^2)$. Then $$\log(E Y) = \beta_0 + \sigma^2/2 + \beta_1 x$$ and $VY = (e^{\sigma^2} - 1) (EY)^2$. Thus this model has the same mean-variance structure as a $\Gamma$-model with log-link and dispersion parameter $\psi = e^{\sigma^2} - 1$. If we keep the reparametrization of the intercept in mind, this model can be fitted using `glm` with `family = Gamma("log")`. The parameter estimates will not be the MLE, and the estimates are not efficient when compared to simply log- transforming $Y$ and using `lm`.               sigma <- 4     beta <- 1     beta0 <- 1     n <- 100     x <- rnorm(n)     y <- exp(beta0 + beta * x + rnorm(n, 0, sigma))     fit1 <- lm(log(y) ~ x)     sigmasqhat <- sum(fit1$residuals^2) / (n-2)     betastart <- coef(fit1) + c(sigmasqhat / 2, 0)     fit <- glm(y ~ x,             start =  betastart,            family = Gamma("log"))      The code above will need to be run a couple of times to discover a problematic data set. I get an error around 2-3% of the times. More frequently, the algorithm does not converge with the default maximal number of iterations (which is 25), but this problem can easily be fixed by increasing `maxit`. The errors I encounter are `NA/NaN/Inf in 'x'` and `inner loop 1; cannot correct step size`. The example above was considered for educational purposes to investigate differences and similarities between the log-normal and the $\Gamma$ models. I also have a couple of data sets where a $\Gamma$-model with a log-link seems to fit well. It is not a problem to fit the model to the data set itself, but once I bootstrap I start encountering similar convergence problems on a few of the bootstrapped data sets.


title: Can I compare ordinal rankings (and if so, how)?

text: I have 15 supervisors who have ranked 21 employees from 1 (best) to 21 (worst), but the rankings are strictly ordinal-there is no interval data available. Can I compare any two (or more) supervisors' rankings and determine how closely their evaluations are aligned? I'll be using Excel 2010 (or pen/paper) for my analysis, I don't have access to statistical software. Sample Data (Is supervisor 3 more similar to 1 or 2)?                        Ranking    Supervisor 1    Supervisor 2    Supervisor 3     Employee 1          1               3               1     Employee 2          2               2               3     Employee 3          3               1               2


title: Interpretation of multinomial logistic regression output from R

text: I have used `mlogit` package and I am trying to summarize the results I have from my model. I have a question regarding the reference value and will get to that in a moment.               redata.full <- mlogit(no.C~ 1| WR+age+age2+BP+noC.1yr, data=redata, reflevel="0", na.action=na.fail)          no.C = number of offspring         WR = risk     age+age2 = the non-linear relationship that as an individual ages their production decreases     BP = browsing pressure     noC.1yr = number of offspring produced the year before      I recognize that my data is ordinal in nature, but Im following other people's methods who have done this and used the reference based approach rather than ordinal logistic regression. However, I am still shakey on justification other than citing the other person and saying "he did it too!" If anyone has a suggestion I would appreciate it. My results for this model are:               Call:     mlogit(formula = no.C ~ 1 | WR + age + age2 + BP + noC.1yr, data = redata,          na.action = na.fail, reflevel = "0", method = "nr", print.level = 0)          Frequencies of alternatives:            0        1        2      0.233766 0.675325 0.090909           nr method     5 iterations, 0h:0m:0s      g'(-H)^-1g = 2.16E-07      gradient close to zero           Coefficients :                    Estimate Std. Error t-value Pr(>|t|)       1:(intercept) -0.281226   1.225763 -0.2294  0.81854       2:(intercept) -0.605312   1.997179 -0.3031  0.76183       1:WR           0.847273   0.518854  1.6330  0.10248       2:WR           1.347976   0.689916  1.9538  0.05072 .     1:age          0.314075   0.275486  1.1401  0.25425       2:age         -0.422368   0.395240 -1.0686  0.28523       1:age2        -0.018998   0.014446 -1.3151  0.18847       2:age2         0.022572   0.018949  1.1912  0.23359       1:BP          -0.143720   0.173585 -0.8280  0.40770       2:BP          -0.074553   0.331108 -0.2252  0.82185       1:noC.1yr      0.574304   0.377821  1.5200  0.12850       2:noC.1yr      1.251673   0.626033  1.9994  0.04557 *     ---     Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1          Log-Likelihood: -116.6     McFadden R^2:  0.079844      Likelihood ratio test : chisq = 20.236 (p.value = 0.0271)          exp(cbind(OddsRatio = coef(redata.full), ci))                   OddsRatio      2.5 %    97.5 %     1:(intercept) 0.7548580 0.06831155  8.341351     2:(intercept) 0.5459038 0.01089217 27.360107     1:WR          2.3332750 0.84394900  6.450831     2:WR          3.8496270 0.99577472 14.882511     1:age         1.3689929 0.79782462  2.349065     2:age         0.6554925 0.30209181  1.422317     1:age2        0.9811815 0.95379086  1.009359     2:age2        1.0228284 0.98553735  1.061530     1:BP          0.8661299 0.61634947  1.217136     2:BP          0.9281585 0.48504538  1.776078     1:noC.1yr     1.7758933 0.84686698  3.724076     2:noC.1yr     3.4961862 1.02497823 11.925441      I would like confirmation of my interpretations: The model is better than a null - obtained from the likelihood ratio test. Question: How do I test how well the model is actually working (i.e., goodness of fit)? Hosmer-Lemshow test? Ive read warnings about using the McFaddin's Pseudo R where they really aren't applicable to multinomial regressions. Ive found a HL test with `ResourceSelection` library and it says my model is NOT doing well at all. Now what? Interpretation: WR and noC.1yr are the only variables that are coming out as slightly significant. But this is only between the reference value of 0 and production of 2 calves. It is not significantly different between 0 or 1 for these variables. Question: Ive been trying to find somewhere in the vignette what the t-value is - it is just a t-test? How would I refer to the estimate as being significant? "The estimated odds for 2-offspring being produced versus 0 were 3.85 (95% CI = 1.0-14.88) which was significant (t= 1.99, P=0.05)" Referring to my statement regarding setting the reference value. When I run this exact same model using my other options of 0 or 1 offspring - I get completely different results of which variables are significant. If I use 2 as the reference value then Age+WR+noC.yr are significant. If I use 1, then Age only is sig. So, which one to use? I have read you want to pick one that is most relevant to your hypothesis, but in this case I could motivate any of the 3 levels.


title: ARIMA model improvement

text: Suppose an ARIMA(p,d,q) model fits the data very poorly. What are some ways to improve the fit? Is there any way to do it without guessing and checking? **Edit.** I am using auto.arima() to search for an ARIMA model. But it is a poor fit.


title: How to predict future reservations when data for the current day is incomplete?

text: I'm trying to build a model to predict reservations up to 15 days in advance. So, if I want to predict how many reservations there will be tomorrow, I use historical data of how many total reservations there were on two days prior since when using the model to forecast for tomorrow, today won't be finished... Does that make sense? I think it is an OK and unbiased model, however, it does not use all of the available data, namely, how many reservations for tomorrow have been made today. So, if I want my model to be more 'realtime' and account for current reervations as well, how do I do that? Do I look at reservation data over the course of the day, and just set up a proportion or something? For example, if I use my model at 13:00, find that I have 20 additional reservations, and I know that historically 40% of additional reservations are made by 13:00, do I just take 20+ 20*(60/40) = 50. So now I should 'expect' 50 more on top of what yesterdays total told me? This is all I can think of. **Edit:** " For example if it takes 1 hour to complete three innings at a baseball game you can rest assured that the next 6 innings is going to take a lot longer than 2 hours to complete." Well, what if you knew that historically the first 3 innings accounted for 25% of the time? Sticking with the ballgame example, that's good illustration of my problem which I think you understand Dave, but I'm not sure I've made it clear for the others... If you wanted to predict the duration of a ball game, you might have some formula based on the current teams, the pitchers, etc.. But, consider a game where you predicted a 2.5 hour game, and you're in the 4th inning, and it's already taken 3 hours, now what do you do? Dave, I've tried different techniques, involving the seasonal arima methods, etc. I'm not familiar with 'level shifts', although I've seen you mention them in many of your posts. I will do some further web surfing to understand this concept. I'd be very interested in a chat session. Please let me know when you're available, I'm +11 hours EST.


title: Add a column to a dataframe based on a probability distribution

text: Suppose I have two dependent categorical attributes A, B. I have a dataframe X that lists probabilities (or expected counts) for all combinations of all categories of A and B. Let's assume two categories in each attribute, the dataframe could then look like this:               X <- data.frame(A = c(1,1,2,2), B = 1:2, PROB = c(0.1,0.3,0.2,0.4))     X      Output:                 A B PROB     1 1 1   .1     2 1 2   .3     3 2 1   .2     4 2 2   .4      Now I would like to add a B column to another dataframe Y that has only an A column. Y exists and has many more rows than X, and also additional columns. A synthetic Y could look like this:               Y <- data.frame(A=sample(2, size=10000, replace=TRUE),                     B=NA,                     C=sample(5, size=10000, replace=TRUE),                     D=sample(3, size=10000, replace=TRUE))      The contents of the new B column should be sampled at random using the conditional probabilities for B given the value of the A column. I am new to R, and I was wondering whether this operation can be handled more elegantly than writing a loop that iterates over all A categories. Perhaps a dedicated model class for which `predict` could be applied? Also, I lack the knowledge of statistical terminology -- I would appreciate any hint on how this kind of imputation (?) is referred to in the literature.


title: Pseudorandom sequence whose values can be computed directly for given seed?

text: Wikipedia says that Blum Blum Shub is a pseudorandom generator whose ith value can be computed directly for a given seed. Wikipedia also says that this generator is too slow to be practical for simulation applications. Is there any standard generator with the direct computation property that is fast enough for Monte Carlo? Is there a reason this cannot exist?


title: Computer specs for R?

text: What should I look for when putting together a computer specifically for runnig R and similar programs. A lot of CPU power and RAM are given, but what about such things as an SSD and can R run on the graphics card? (like OpenCL..) Thanks in advance.


title: Sampling under normal distribution

text: The concentration of impurities in a semiconducter in the production of microprocessors for computer is a normally distributed random variable with mean 127 parts per million and standard deviation 22 parts per million. A semiconducter is acceptable only if its concentration of impurities is below 150 parts per million. What is the proportion of semiconducters that are acceptable for use?(the area under normal curve for the value of $Z=1.5$ is $0.668$).


title: significance of coefficients and significance of marginal effects

text: Suppose I have a non-linear model, say probit/logit, how can I understand the significance of a coefficient as opposed to the significance of it's marginal effect? Say, I just need to know the _importance_ of a variable, would the former suffice?


title: Knowledge based annotation vs. Corpus based annotation

text: What is the difference between knowledge based annotation and corpus based annotation in the context of training a Machine Learning algorithm?


title: Mixing dichotomous and polytomous items in one test

text: I have to evaluate a test paper which contains both dichotomous and polytomous items. I am currently using R. But in R IRT models for scoring of dichotomous and polytomous items are different and I have not found any way of combining both type of items so how can I mix both dichotomous and polytomous items into one. Specifically -   1. Is there any programme/function in R that handles both dichotomous and polytomous items    2. If no direct function can combine both dichotomous and polytomous items in R then is there any indirect way of doing it    3. If not 1 and 2 then any other open source IRT software which can handle dichotomous and polytomous items together


title: Values in predictors have same pattern for two possible value in dependent variable?

text: I am trying to build a model for auto rejection of crowd-sourced task (eg: human name transcription of census data). My data set has 5 predictors, and dependent variable is either correct or incorrect. However, there are few issues I have had.   1. Data set is unbalanced with ~5% (up to 10% sometimes) of incorrect instances and 95% of correct instances.   2. Also, the patterns in predictors for incorrect and correct tasks are similar. I tried logistic regression, but it didn't work. But, it is very important for me to be able to predict those 5% of incorrect tasks since it is very critical. In such situation what should be the direction to move forward? I would appreciate your help. Thanks.


title: How to find the pdf of a time-dependent parameter in a linear regression model with drift?

text: Suppose there exists a simple linear equation, where the dependent variable y(time series; measurements available) depends on x1 and x2. If the parameter multiplier of x2 is time dependent and the multiplier of x1 is a deterministic variable. How to find the pdf of the parameter multiplier of x2. Given the time series observations of y, x1 and x2. Please tell me how to solve this. Kalman/particle filter?


title: When do we use multinomial regression and Poisson regression?

text: They are both regression methods for discrete dependent variables. I have a discrete dependent variables like $\{A, B, C, D, E\}$, and $A$, $B$, $C$, $D$, $E$ can be cardinal data(10, 20, 30...) or non-cardinal data(like car, train, and ship). Which model is proper for each case?


title: Linear regression with time-varying parameters

text: I am attempting to estimate the following model in Stata (the only statistics package I'm familiar with): $$ Y_{t}=\alpha_{t} + \beta_{t}X_{t}+\varepsilon_{t} $$ where $X$ and $Y$ are series and $\alpha_{t}$ and $\beta_{t}$ are time-varying parameters defined in state equations: $$ \alpha_{t}=\alpha_{t-1}+\eta_{1t} \\ \beta_{t}=\beta_{t-1}+\eta_{2t} $$ which are assumed to follow a random walk. $\varepsilon$, $\eta_{1t}$, $\eta_{2t}$ are assumed to be gaussian independent white noise. I am interested in finding the value of $\beta$ over time. The paper from which I am working suggests estimation by the Kalman filter. I have attempted to understand and implement the `sspace` function here, but I believe it is incapable of handling the time variation. Is any of this achievable in Stata 11?


title: How to test that a particular parameter in a nonlinear model does not vary across subpopulations?

text: My original model is: $Y=\alpha(\delta*X_1^{-\rho}+(1-\delta)*X_2^{-\rho})^{-\eta/\rho}exp(e_t)$ Which I partially linearize to: $ln(Y)=ln(\alpha)-(\eta/\rho)*ln(\delta*X_1^{-\rho}+(1-\delta)*X_2^{-\rho})+e_t$ I've used a nonlinear least squares algorithm to estimate the parameters. Now, I am trying to determine if the "intercept" in this partly linearized model (in this case $ln(\alpha)$) is different for the different years of data that I have. I have 3 years of data, with about 24 observations in each year. My first thought is to just add a dummy variable for each year as $ln(\beta*Year)$ and just test if the two coefficients that would result (three years minus one comparison years) are jointly statistically different from zero using an F-test. But I doubt that this would be legit since I am trying to test if one parameter ($ln(\alpha)$) varies over years rather than if all the parameters vary over the years. My second thought is to just split up the sample by year, run three different regressions, and somehow test if all the $ln(\alpha)$ parameters are different. I don't what test I would use or how I would use it, though, to do this. I'm thinking maybe a t-test with unequal variance on each pair of the coefficients, but I feel like pairwise tests on each parameter does not get me a legitimate "joint" test. I have also read about the Chow test, but again I want to test if just one parameter is different across the subpopulations, not all the parameters, which is what the Chow test does. Any thoughts on what the right procedure is? Thanks!


title: What can be concluded from small sample size?

text: I had just conducted a test where the customer wants to see a 20% improvement in perforation length from a baseline perforator. The baseline perforator has not been tested in the customer's specific conditions in which he will eventually perforate in. The customer plans to shoot 40,000 perforations in the future. I'm not sure if that bit of information is important but I wanted to be sure that it is understood that the true population mean of perforation length does not exist yet. Under the specific conditions tailored for the customer the baseline perforator was shot 3 times with the following results: [7.68, 8.72, 8.08] inches. Three new perforator designs were constructed in hopes to out perform the baseline perforator and each design was shot 3 times with the following results:   * D1 [8.04, 8.93, 10.05] inches   * D2 [11.15, 6.91, 9.21] inches   * D3 [7.1, 7.6, 7.02] inches As I am relatively new to statistics I was wondering if some of you professionals can help me out with what information can be gained from these tests results. If you can point me in the right direction as in what methods I should refer to in gaining these insights you suggest would be of much help. I know that effort is encouraged on my part so I would like to include some of my initial attempts at inferential statistics and maybe you can correct me where I am wrong or include what more I should/can do. Estimating Sample Size Requirements: Solving for the sample size requires the population standard deviation, $\sigma$. Most often we do not know it so we have to use an estimate or “planning value” in its place. Here are a few options:   1. Estimate σ from previous studies using the same population of interest   2. Conduct a pilot study to select a preliminary sample. Use the sample standard deviation from the pilot study   3. Use a judgment or “best guess” for $\sigma$. A common “guess” is the data range (high value – low value) divided by 4 The question I asked is: “What minimum sample size is necessary to produce 95% confidence (assuming a normal distribution) that the sample mean is $\pm$ 1 inch of the true population mean, $\mu$?” Since we do not have an estimate of $\sigma$ from previous studies we will go with option 2 and consider our information as pilot studies. Test charge Design 2 perforation length data is as follows: [11.15, 9.21, 6.91] inches. The sample mean is 9.09 inches and sample standard deviation is 2.12 inches. How large a sample should be shot to provide a 95% confidence interval with a margin of error (E) of 1 inch, assuming the population standard deviation is equal to the sample standard deviation? $$n = \frac{(z_{\alpha /2})^2 (\sigma)^2}{E^2}$$ $$n = \frac{(1.96)^2 (2.12)^2}{(1)^2}$$ $$n = 17$$ Interpretation: To have 95% of our sample means contain the true population mean, $\mu$, we need a sample size of 17 shots. Or, 95 out of 100 tests made up of 17 perforation shots will contain within its average length $\pm$ 1 inch the perforation design’s true average perforation length. It is not clear to me, however, how to move forward with this result. Does it mean I need to make an additional 17 shots as it's own sample or do I need to make 14 more shots to add to the original 3 shots? How do I handle the new information? I would suspect that the standard error of the mean would narrow with more measurements and would thus adjust the required sample size? If we do not assume that the population standard deviation is equal to the sample standard deviation we can set the planning value for the population standard deviation as the Range of the perforation lengths we observed divided by 4 (option 3). From there we can estimate the minimum number of shots that should be taken with the margin of error of 1 inch to contain the true population mean 95% of the time. Test charge Design 2 perforation length data is as follows: [11.15, 9.21, 6.91] inches. $$\text{Planning Value for Std. Dev.} = \frac{\text{Range}}{4}$$ $$\text{PV} = \frac{(11.15-6.91)}{4} = 1.06$$ $$n = \frac{(z_{\alpha /2})^2 (\sigma)^2}{E^2}$$ $$n = \frac{(1.96)^2 (1.06)^2}{(1)^2}$$ $$n = 4.3$$ Interpretation: To have 95% of our sample means contain the true population mean, $\mu$, we need a sample size of 4 shots. Or, 95 out of 100 tests made up of 4 perforation shots will contain within its average length $\pm$ 1 inch the perforation design’s true average perforation length. So which would be the more appropriate method (option 2 or 3)? And how should that information be used? In both cases more shots were required. But what does it tell me if the opposite occurs? e.g. for charge Design 3, using both option 2 and 3 method, I got n = 0.37 shots and n = 0.08 shots respectively. Lastly, I tried to estimate what the %chance that the new charge design's true-population-mean-perforation-length will meet the required improvement. Taking the mean of the baseline shots to be 8.16 inches and adding 20% to it I found the goal length to be $\approx 9.792$ inches. By taking the sample means of each design we see that none of them averaged to meet the required length. However, I found the 95% confidence interval using a t-distribution for Design 2 to be [3.81 to 14.36] inches. Using excel Goal Seek I believe there is a 31% chance that Design 2's true perforation mean length will meet the 9.792 inch requirement. I'm not sure if my method was correct but maybe you can tell me if I am wrong or right or if interested, I will spell out how I went about that estimation.


title: Generate pairs of random numbers uniformly distributed and correlated

text: I would like to generate pairs of random numbers with certain correlation. However, the usual approach of using a linear combination of two normal variables is not valid here, because a linear combination of uniform variables is not any more an uniformly distributed variable. I need the two variables to be uniform. Any idea on how to generate pairs of uniform variables with a given correlation?


title: Weighted average of percentage increase

text: I believe I have quite a simple problem, but want clarification on whether it is the best method to use. Say I have have insurance line, and the net income (£) from this business for 2011, 2012, 2013 is 1,000,000; 4,500,000; and 5,000,000 respectively. Hence, the percentage increase from 11-12 is 350%, and from 12-13 is 11.1%. Now, the average of these is 180.5%, which is a representation of the average percentage increase across all years, but it has been suggested to me to instead represent the data by a “weighted” average. I am curious… would the weights be the difference in income from each year? i.e. Weighted Avg. = $((3.5)*350\% + (0.5)*11.1\%)/4 = 307\%$ or would the weights possibly the values themselves? Or maybe there is another option? Maybe a standard average is sufficient? Thanks very much for your inputs.


title: Handling borderline cases of the proportional hazards assumption

text: In order to verify the assumption of proportional hazards, I plotted the following: $$\log(-\log(S(t))) = \log(\Lambda(t))$$ If the hazards are indeed proportional for two groups, these curves should be parallel with the vertical distance being the log hazard ratio (which is constant if the proportional hazards assumption holds). I used the following R code to make such a plot (using the `survival` and `rms` packages):               foo <- survfit(Surv(TIME, HADEVENT)~CASE, data=bar)     survplot(foo, loglog=T, xlim=c(0,5*365.25), ylim=c(-10,-2), col=c("red","blue"))      I obtained the figure shown below for the groups in my study, and I am unsure what to make of it. For small values of $t$ ($t \leq 120$ days), the proportional hazards assumption does not hold at all based on visual inspection. For larger values ($ \gt 1$ year), it seems to be very reasonable. ![bla](http://i.stack.imgur.com/RT54b.png) I also fitted a Cox PH model and used `cox.zph` to assess the proportionality assumption, with the following results:               foocox <- coxph(Surv(TIME, HADEVENT)~CASE, data=bar)     print(cox.zph(foocox))              rho chisq        p     CASE -0.0372    46 1.18e-11      Based on the $p$-value, I conclude that the null hypothesis of proportionality is strongly rejected. Here is the `cox.zph` plot, which also doesn't look too promising (e.g. it should be constant): ![coxzph](http://i.stack.imgur.com/Nnm5e.png) Here's a frequency table for my data (`CASE=0` is a 1:1 matched sample for case `CASE=1`):               CASE  event no event       n     0    10.464   93.832 104.296     1    22.903   81.393 104.296      **Edit** : after filtering out very large times -- which were inaccurate -- I obtain a highly significant $p$-value via `cox.zph`, indicating the proportionality assumption is violated. What are my options given this situation? Is it justified to continue working with proportional hazards models based on the graph of `cox.zph` for larger values of $t$? Should I somehow refrain from using them for $t \lt 1$ year based on the visual inspection?


title: Converting (normalizing) very small likelihood values to probability

text: I am writing an algorithm where, given a model, I compute likelihoods for a list of datasets and then need to normalize (to probability) each one of the likelihood. So something like [0.00043, 0.00004, 0.00321] might be converted to may be like [0.2, 0.03, 0.77]. My problem is that the log likelihoods, I am working with, are quite small (for instance, in log space, the values are like -269647.432, -231444.981 etc). In my C++ code, when I try to add two of them (by taking their exponent) I get an answer of "Inf". I tried to add them in log-space (Summation/Subtraction of log), but again stumbled upon the same problem. Can anybody share his/her expert opinion on this?


title: How to calculate the confidence intervals for likelihood ratios from a 2x2 table in the presence of cells with zeroes

text: I am analysing a diagnostic test (against a gold standard, using a 2x2 table). I want to calculate likelihood ratios (sensitivity / (1-specificity) etc) however I have several sets of data with 0 false positives therefore a specificity of 1.... I have come across other data (e.g. in evidence bases for clinical guidelines for investigation of periprosthetic infection of hip or knee by AAOS) where I know that the authors had 0 false positives but have still presented confidence intervals for the likelihood ratios. Does anyone know how they have adjusted their calculations to allow this? Thanks


title: Skewness in residuals

text: I have skewness of -0.5 in the residuals which does not seem to be improve much after using logs/root. Is this a major concern? Tests relating to homoskedasticity, multicollinearity, outliers all seem fine. Perhaps interpreting the regression using robust standard errors would help..


title: Probability that two values drawn at random from a normal distribution are separated by at most T

text: I need a closed-form expression for the probability that 2 values drawn at random from a normal distribution are separated by at most T, as a function of T and the variance of the distribution. I can formulate the problem as an integration, but can't solve it. I'd be happy with a well-motivated approximation. I'd be amenable to making some simplifying assumptions about the distribution, if that would help.


title: goodness of fit

text: A random sample of size `n = 40` is taken. Conduct a goodness of fit test to see whether this sample appears to have been selected from a normal distribution. The number of intervals for the sample data is equal to:               4     5            8            10           20      How do you calculate this?


title: What kind of model is this and how do you estimate its parameters?

text: For the following model: $$T = A_0 + A_1M+ A_2M^2 + A_3\ln(P) + A_4\ln(P)^2+A_5M\ln(P)$$ where: > $T$ = Dependent variable > > $P$ = Independent variable > > $M$ = Independent variable > > $A_0 \dots A_5$ = Coefficients of the regression. I would like to know which class of regression model this equation corresponds to, when you use this kind of model, and how you find the parameters $A_0 \dots A_n$.


title: How to interpret this correspondence analysis plot with individuals as nominal variable?

text: **Background:** Although correspondence analysis is used mainly for visualizing similarities of categories of two or more nominal variables, I tried following: 39 students (from the same survey as in this question) were asked to ask their friends about smoking habits with four possible responses (I do smoke; I used to, but I do not anymore; I have tried and I have not even tried). Then I created matrix with two columns: First is nominal variable with 39 categories (my students, in further labeled by numbers) and second one is also nominal with 4 categories (smoking levels). This data were used as an input into correspondence analysis and resulting plot can be seen below. **Questions:** I am not sure about some interpretations of the plot.   * Can we say something about frequencies of categories of second (smoking) variable? It looks for example that there is higher number of respondents near the category "I have tried". Does it mean that most of participants have only tried smoking?   * If two respondents are close together (for example 6 and 18), does it mean that their friends have similar smoking "behavior"?   * Where in the map are for example those respondents whose friends mostly used to smoke, but they do not anymore? Are they those who are near the label "I used to, but I do not anymore" (respondents 3, 22, ...) or those whose label has similar angle with x -axis as label "I used to, but I do not anymore" (respondents 3, 22, 30, 29, 10, 21, ...) or those whose projection to any axis is similar as projection of label "I used to, but I do not anymore" (respondents 5, 8, 2, 25 when projected to x - axis)?   * Does the plot below say something about the degree of dependence between two variables? How should plot below look like if variables were independent? ![Correspondence analysis](http://i.stack.imgur.com/zpFUe.jpg)


title: IQ adaptive test items in 1pl, 2pl or 3pl IRT model?

text: Some adapative test systems (e.g. school assessment tools) use the 1pl IRT model, while others use the 2pl or the 3pl. When developing an adaptive IQ test, is there a rule of thumb about which model to choose in calibrating the item difficulty and test takers ability? I can't find any research that gives some insights in fit between IQ test items and different kinds of IRT models. Many thanks in advance!


title: How to average quantized and truncated data?

text: So I have data that has been quantized by an analogue to digital converter. (continuous data has been turned into discrete data and the values range from 0 to the saturation value , which is 127 in this case). This particular instrument I used to gather the data is quite noisy, let's say there is added Gaussian noise to the real value. Luckily, when taking single measurements, I have enough time to take multiple measurements and average them to reduce the noise. Note that sampling rate is not an issue here since the thing that I'm taking measurements of is completely stable. Obviously , taking the simple mean will produce a biased result because values cannot exceed 0 or 127 (for example, if you attempt to use a plain old averaging on something with a "real" value of 126, you will get an estimated value that is less than 126. This is because the added Gaussian noise will not give you any value higher than 127 because of the truncation). So how do I take the average so the result will give me an unbiased estimator of the real value?


title: Getting the bootstrap-validated AUC in R

text: In a paper by Faraklas et al, the researchers create a Necrotizing Soft-Tissue Infection Mortality Risk Calculator. They use logistic regression to create a model with mortality from necrotizing soft-tissue infection as the main outcome and then calculate the area under the curve (AUC). They use the bootstrap method to find the "bootstrap optimism-corrected ROC area." If I were to do this in `R`, how would it look like? The code I have been toying with looks something like below:               library(boot)     library(ROCR)          auc_calc <- function(data, indices, outcomes) {       d <- data[indices,]       # Using glm for logistic regression       # Do I recreate the glm model for each dataset?       fit <- glm(outcomes[indices,] ~ X1 + X2 + X3, data=d, family=binomial)       fit.predict <- predict(fit, type="response")            # Using ROCR to calculate AUC       pred <- prediction(fit.predict, outcomes[indices,])       perf <- performance(pred, "auc")            # Returning the AUC       return(perf@y.values[[1]])     }          boot.results <- boot(data=my.data, statistic=auc_calc, R=10000, outcomes=my.outcomes)      Is this correct? Or am I doing something wrong - namely should I be passing in a glm model rather than recalculating it each time? As always thanks for the help.


title: How to get exact distribution of estimated p for binomial distribution?

text: This question is kind of a follow up of another question I had: Asymptotic normal distribution via the central limit theorem There I had to calculate the estimator for $p$ (meaning $p$ for success) and approximate it's distribution by approximation with a normal distribution. Now I would like to get the exact distribution of $p$. I got already the following hint: "You have the functional form of $\hat{p}$. Look up how we derive the distribution of a function of a discrete random variable." Unfortunately that did not lead me to a solution...


title: Evaluation of classifiers: learning curves vs ROC curves

text: I would like to compare 2 different classifiers for a multiclass text classification problem that use large training datasets. I am doubting whether I should use ROC curves or learning curves to compare the 2 classifiers. On one hand, learning curves are useful for deciding the size of the training dataset, since you can find the size of the dataset at which the classifier stops learning (and maybe degrades). So the best classifier in this case might be the one reaching the highest accuracy with the smallest dataset size. On the other hand, ROC curves let you find a point with the right trade-off between between sensitivity/specificity. The best classifier in this case is just the one closer to the top-left part, with the highest TPR for any FPR. Should I use both evaluation methods? Is it possible for a method with a better learning curve to have a worse ROC curve, and vice-versa?


title: Robust version of Hotelling $T^2$ test

text: I am looking for a robust version of Hotelling's $T^2$ test for the mean of a vector. As data, I have a $m\ \times\ n$ matrix, $X$, each row an i.i.d. sample of an $n$-dimensional RV, $x$. The null hypothesis I wish to test is $E[x] = \mu$, where $\mu$ is a fixed $n$-dimensional vector. The classical Hotelling test appears to be susceptible to non-normality in the distribution of $x$ (just as the 1-d analogue, the Student t-test is susceptible to skew and kurtosis). what is the state of the art robust version of this test? I am looking for something relatively fast and conceptually simple. There was a paper in COMPSTAT 2008 on the topic, but I do not have access to the proceedings. Any help?


title: Estimating variability over time

text: I have a dataset that contains ~7,500 blood tests from ~2,500 individuals. I'm trying to find out if variability in the blood tests increases or decreases with the time between two tests. For example - I draw your blood for the baseline test, then immediately draw a second sample. Six months later, I draw another sample. One might expect the difference between the baseline and the immediate repeat tests to be smaller than the difference between the baseline and the six-month test. Each point on the plot below reflects the difference between two tests. X is the number of days between two tests; Y is the size of the difference between the two tests. As you can see, tests aren't evenly distributed along X - the study wasn't designed to address this question, really. Because the points are so heavily stacked at the mean, I've included 95% (blue) and 99% (red) quantile lines, based on 28-day windows. These are obviously pulled around by the more extreme points, but you get the idea. ![alt text](http://a.imageshack.us/img175/6595/diffsbydays.png) It looks to me like the variability is fairly stable. If anything, it's higher when the test is repeated within a short period - that's terribly counterintuitive. How can I address this in a systematic way, accounting for varying n at each time point (and some periods with no tests at all)? Your ideas are greatly appreciated. Just for reference, this is the distribution of the number of days between test and retest: ![alt text](http://a.imageshack.us/img697/6572/testsateachtimepoint.png)


title: Distribution. How to work out what distribution data fits from visual analogue scales

text: I am having some problems figuring out what distribution my data best fits into. My data is from a 100 mm visual analogue scale, so data points range from 0 to 100. Participants were asked to score an animal's 'lameness' using these scales. There were 4 animals, all with varying severity of lameness. Am I correct in assuming my data could be considered: 1. Continuous - as it's measured on a continuous scale OR 2. Binomial (with 100 binomial totals) as something either did or did not happen at each point on the scale I am using a GLMM to make my model. I have run the data through using gamma (had to recode my zero values), Poisson, and binomial. And bar eye balling the results on a bar chart of the data, I have no idea which is correct.


title: Cheaper/faster method to estimate uncertainties than bootstrap

text: I'm using a genetic algorithm (GA) to estimate the minimum value of a likelihood function $L[x]$ which is too complicated to evaluate mathematically. This likelihood function quantifies the goodness of fit between a model (which depends on the $x$ parameter) and a set of $M$ observations $O=\{a_1, a_2, ..., a_M\}$. After enough iterations (when the optimal value found by the GA hasn't changed for a while) I stop the GA and assume the value $x_0$ for the model parameter $x$ it settled on, is one that satisfies: $$L[x_0] \sim min(L[x])$$ which identifies the "best" model/observation fit. Now, I want to assign an uncertainty to this $x_0$ value for which I apply a bootstrap process with replacement. This means I run the GA again a number $N$ of times (each time using a random sample of $M$ observations taken from $O$) which results in $N$ values for the $x$ parameter: $(x_1, x_2, ..., x_N)$. Finally I assign the standard deviation of the $N$ values returned by the GA as the uncertainty (error) of $x_0$: $$e_{x_0} = std(x_1, x_2, ..., x_N)$$ This process works fine, but since I am forced to run the GA a large $N$ number of times, this can quickly become unmanageable in computational time. The question is: is there a faster/computationally cheaper way of assigning uncertainty to a parameter value estimated through a genetic algorithm?


title: How to test if reading from two devices are significantly different?

text: I am trying to compare readings from an old flowmeter and a new digital flowmeter to find out if they are significantly different. With the new meter I did readings under different conditions eg. calibrated at 30 seconds, calibrated at 60 seconds, uncalibrated at 30 seconds and uncalibrated at 60 seconds. I want to compare my old flowmeter (first column on dataset) to the other readings (columns 2:5 on dataset). I used anova to compare the groups and found some differences but I am stuck with TukeyHSD...My question is: Is this the correct way to do the comparisons...I would appreciate some advice on how to accomplish this task...Here is my dataset and code....                   meter <- structure(list(SiteNumber = structure(c(1L, 1L, 1L, 1L, 1L, 1L,      1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L,      2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 3L, 3L, 3L, 3L,      3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 4L, 4L, 4L,      4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L), .Label = c("site1",      "site2", "site3", "site4"), class = "factor"), old_dev = c(2.09,      1.927, 2.11, 2.114, 2.11, 1.995, 2.605, 2.563, 2.393, 2.286,      2.529, 2.371, 2.214, 2.128, 2.302, 2.204, 2.274, 3.513, 2.61,      3.502, 3.155, 3.131, 3.55, 3.656, 3.666, 3.618, 3.023, 3.226,      3.649, 3.519, 3.322, 3.352, 3.323, 3.144, 2.869, 3.383, 3.076,      2.978, 3.189, 2.934, 2.813, 3.015, 3.014, 3.056, 2.952, 2.988,      2.976, 3.049, 2.869, 2.965, 2.895, 2.064, 1.333, 1.432, 1.599,      1.52, 1.51, 2.734, 2.182, 2.71, 2.488, 2.242, 1.998, 1.175, 1.791,      2.152, 3.414, 3.278), new_dev_unc30 = c(2.6, 2.6, 2.7, 2.6, 2.3,      2.6, 3.2, 2.9, 2.1, 2.8, 1.8, 2.6, 2.5, 1.4, 2.7, 2.6, 2.2, 4.1,      3.4, 3.7, 4, 4, 4.5, 4.5, 1.7, 2.9, 2.4, 2.9, 3.7, 3.3, 4, 1.7,      3.5, 1.9, 3.3, 4.1, 3.7, 4, 3.6, 2, 2.4, 3.7, 2, 2.8, 3, 1.8,      1.8, 3, 1.9, 3.1, 1.1, 2.2, 1.6, 1.2, 1.4, 1.3, 1.4, 2, 1.9,      1.4, 0.7, 2.1, 1.2, 1.8, 1.3, 1.2, 0.8, 2.5), new_dev_unc60 = c(2.7,      2.2, 2.6, 2.3, 2.7, 2.4, 3.1, 2.8, 2.3, 2.9, 1.8, 2.4, 2.6, 1.3,      2.7, 2.6, 2.2, 4.1, 3.8, 3.9, 4.5, 4, 4.2, 4.6, 2, 2.6, 2.6,      3, 3.4, 3.5, 3.9, 2, 3.4, 2.3, 3.3, 4.2, 3.7, 3.9, 3.5, 2.1,      2.4, 3.1, 2.1, 3, 3.2, 1.9, 1.9, 3, 2.1, 3.1, 1.3, 1.5, 1.6,      1.2, 1.8, 1.1, 1.2, 2.1, 2, 1.4, 0.8, 1.9, 1.1, 0.9, 1.4, 1.5,      0.8, 2), new_dev_cal30 = c(2.4, 2.5, 2.5, 2.8, 2.4, 2.3, 3.1,      2.6, 2, 3.5, 1.6, 2.5, 2.5, 1.7, 2.9, 2.5, 2.4, 4.5, 3.3, 4.3,      4.3, 4.3, 4.4, 4.5, 1.1, 1.9, 4, 3.5, 3.1, 3.8, 3.5, 2.2, 3.5,      2, 3.6, 3.7, 4, 3.5, 3.8, 2.4, 2.7, 3.3, 2.7, 2.3, 3.2, 2.3,      0.9, 3.2, 2.9, 3.1, 0.9, 2.3, 1.3, 1.4, 1.5, 1.1, 1.3, 2.7, 1.7,      2.1, 0.7, 1.9, 1.2, 1, 1.6, 1.1, 1.4, 1.6), new_dev_cal60 = c(2.6,      2.5, 2.4, 2.3, 2.5, 2.3, 3.2, 2.8, 1.9, 3.1, 2.1, 2.7, 2.4, 1.8,      2.7, 2.6, 2.3, 4.3, 3.5, 4.3, 4.3, 4.3, 3.8, 4.4, 1.4, 2, 4.2,      3.3, 3.2, 3.3, 3.4, 2.1, 3.4, 2.5, 3.8, 3.8, 3.9, 3.6, 3.5, 2.4,      2.9, 3.1, 2.3, 2.6, 3.1, 2.2, 0.8, 3.2, 2.6, 3.1, 0.9, 1.7, 1.2,      1.5, 1.2, 1, 1.1, 2.6, 1.6, 1.9, 0.6, 2.3, 1.4, 1.2, 1.5, 1.2,      1.4, 1.4)), .Names = c("SiteNumber", "old_dev", "new_dev_unc30",      "new_dev_unc60", "new_dev_cal30", "new_dev_cal60"), class = "data.frame", row.names = c(NA,      -68L))      Here is what I did:               str(meter)     attach(meter)     boxplot(list(old_dev,new_dev_unc30,new_dev_unc60,new_dev_cal30,new_dev_cal60))     summary(firsttest <- aov(old_dev ~ new_dev_unc30 + new_dev_unc60 + new_dev_cal30 + new_dev_cal60 ))     plot(firsttest)     TukeyHSD(??????)


title: SVM and cross validation with a minimum finding algorithm

text: Just a simple question on parameter selection for SVMs. If I use a minimum finding algorithm to find the optimal parameters for a set of data, how do I "average" the parameters over a set of cross validation runs to come up with the best parameters for my test set run? When doing X-fold cross validation, I can see simply averaging the fit of the model for each parameters set, but not sure how to do it with a minimum finding algo. Thanks


title: Use of standard error of bootstrap distribution

text: _(ignore the R code if needed, as my main question is language-independent)_ If I want to look at the variability of a simple statistic (ex: mean), I know I can do it via theory like:               x = rnorm(50)          # Estimate standard error from theory     summary(lm(x~1))     # same as...     sd(x) / sqrt(length(x))      or with the bootstrap like:               library(boot)          # Estimate standard error from bootstrap     (x.bs = boot(x, function(x, inds) mean(x[inds]), 1000))     # which is simply the standard *deviation* of the bootstrap distribution...     sd(x.bs$t)      * * * However, what I'm wondering is, **can it be useful/valid(?) to look to the standard _error_ of a bootstrap distribution in certain situations?** The situation I'm dealing with is a relatively noisy nonlinear function, such as:               # Simulate dataset     set.seed(12345)     n   = 100     x   = runif(n, 0, 20)     y   = SSasymp(x, 5, 1, -1) + rnorm(n, sd=2)     dat = data.frame(x, y)      Here the model doesn't even converge using the original data set,               > (fit = nls(y ~ SSasymp(x, Asym, R0, lrc), dat))     Error in numericDeriv(form[[3L]], names(ind), env) :        Missing value or an infinity produced when evaluating the model      so the statistics I'm interested in instead are more _stabilized_ estimates of these nls parameters - perhaps their means across a number of bootstrap replications.               # Obtain mean bootstrap nls parameter estimates     fit.bs = boot(dat, function(dat, inds)                   tryCatch(coef(nls(y ~ SSasymp(x, Asym, R0, lrc), dat[inds, ])),                            error=function(e) c(NA, NA, NA)), 100)     pars = colMeans(fit.bs$t, na.rm=T)      Here these are, indeed, in the ball park of what I used to simulate the original data:               > pars     [1]  5.606190  1.859591 -1.390816      A plotted version looks like:               # Plot     with(dat, plot(x, y))          newx = seq(min(x), max(x), len=100)     lines(newx, SSasymp(newx, pars[1], pars[2], pars[3]))          lines(newx, SSasymp(newx, 5, 1, -1), col='red')     legend('bottomright', c('Actual', 'Predicted'), bty='n', lty=1, col=2:1)      ![enter image description here](http://i.stack.imgur.com/xHyqE.png) Now, if I want the variability of these _stabilized_ parameter estimates, I think I can, assuming normality of this bootstrap distribution, just calculate their standard errors:               > apply(fit.bs$t, 2, function(x) sd(x, na.rm=T) / sqrt(length(na.omit(x))))     [1] 0.08369921 0.17230957 0.08386824      **Is this a sensible approach? Is there a better general approach to inference on the parameters of unstable nonlinear models like this?** (I suppose I could instead do a second layer of resampling here, instead of relying on theory for the last bit, but that might take a lot of time depending on the model. Even still, I'm not sure if these standard errors would be useful for anything, since they would approach 0 if I just increase the number of bootstrap replications.) Many thanks, and, by the way, I'm an engineer so please forgive me being a relative novice around here.


title: Determining the right method to compare observed versus measured data

text: I'm trying to solve this problem out loud. I'm not great at statistics but I try sometimes to understand where things are and where they're going. I'm trying to figure out the best method to compare the two following datasets. The first is the actual observation and the second is of what my sensor determined the value was. Basically an object travels across the sensor that with properties. Each object can have anywhere from 2 to 5 components. We know the exact value of each one of those components. The big table below shows the data results from running those objects on the sensor. Observations are independent so observation 5 has no affect on observation 26. The dataset has 31 objects in it, and each object was measure using the sensor. We know the exact properties of the object and I put it under "Actual" **Note:** Total = 1 + 2 + 3 + 4 + 5 All what I'm trying to figure out is if the sensor was adequate or does it need further calibration? The old way of doing things was to get an average difference less than 5% to be acceptable. I don't agree with that methodology as I believe it's prone to error and does not capture the essence of measurements. So therefore, I'm investigating better data analysis tools for the problem. **Hypotheses:**   1. "Is the sensor performing adequately?"   2. Is there an overall statistical difference between Actual and Sensor values?   3. Is there a statistical difference between each observed pair?   4. At which point is the significance broken? i.e. at which confidence interval do the two datasets become different? **Analysis:** The first thing I did was to compare the average differences for each column. I also did the mean of the absolute differences and here are the results:                  Difference between Actual and Sensor (%)                           Total        1        2      3         4        5        Mean of Diff        1.06     2.80    -4.00   7.12    -10.86    17.01        Mean of abs diff    7.31     6.96     7.56   8.52      0.70     1.10      These numbers don't tell me much really but they show somewhat of a relevant comparison. I then proceeded to do a two paired t-test using alpha = 0.05 with the following results:                       df      Mean Actual     Mean Sensor     Var Actual      Var Sensor      pearson     t stat      P 1 tail    t Crit 1 tail   P 2 tail    t Crit 2 tail     Total   30      37,858          89,829,589      89,829,589      77,766,451      0.9463      1.4871      0.0737      1.6973          0.1474      2.0423     1       30      5,411.61        5,261           254,000         407,118         0.7829      2.1076      0.0218      1.6973          0.0435      2.0423     2       30      15,246          15,719          4,762,169       3,878,279       0.7845      -1.9103     0.0328      1.6973          0.0657      2.0423     3       27      17,943          15,396          40,196,438      33,216,653      0.4268      2.0742      0.0239      1.7033          0.0477      2.0518     4       1       6,770.00        7,600           14,257,800      20,480,000      1.0000      -1.5661     0.1809      6.3138          0.3618      12.7062     5       1       5,910.00        4,900           64,800          0.00            N/A         5.6108      0.0561      6.3138          0.1123      12.7062      I basically concluded that the means are not significantly different for "Total", "", "4", and "5". Am I correct on the analysis? I then did an F test analysis for two means to compare those results. Note that I wasn't sure if I needed to do the F-test but I did it to compare the variance. I concluded that according to the F test, it appears that my variance is high on all accounts that the results are statistically significant.                       F   P 1 tail    F Crit 1 tail     Total   1.1551  0.3477      1.8409     1       0.6239  0.1011      0.5432     2       1.2279  0.2887      1.8409     3       1.0299  0.4698      1.9048     4       0.6962  0.4427      0.0062     5       65,535  N/A         161.4476      **Data Table**                           Actual                                             Sensor     Observation   Total       1       2       3       4       5    Total      1       2       3       4       5               1   36820    5860   14840   16120                   36200    5500   15900   14800               2   42530    5160   12000    9840    9440    6090   43300    5200   13100    9200   10800    4900               3   48650    4770   18260   25620                   52000    4700   18300   29000               4   42290    5020   16080   21190                   48800    5600   18400   24800               5   40530    5350   15670   19510                   37800    4800   15300   17700               6   41610    5140   16000   20470                   41500    5700   17100   18700               7   35210    5310   15150   14750                   34000    5100   15600   13300               8   36150    5430   15520   15200                   33600    4900   15400   13300               9   17600    6770   10830                           19300    6600   12800              10   35920    5280   15540   15100                   28700    4700   13000   11000              11   37900    5740   15440   16720                   35100    5500   14700   14900              12   41530    5590   16560   19380                   39300    5400   16600   17200              13   33580    5330   14870   13380                   33400    5100   14800   13500              14   25550    4990   13450    7110                   24600    4800   13000    6800              15   45920    5420   16200   24300                   43300    5100   15800   22400              16   17750    5740   12010                           21700    6500   15200              17   29570    5230   14140   10200                   28800    4600   14600    9600              18   44870    5150   16990   22730                   47500    5600   18900   23000              19   35500    5780   14930   14790                   32200    5600   14200   12500              20   54070    5310   17170   31590                   50100    5000   16700   28300              21   35440    4590   13920   16930                   35000    3900   14400   16600              22   40680    5330   12840   12680    4100    5730   40000    5400   13400   11900    4400    4900              23   16540    7160    9380                           20200    6900   13300              24   47820    4960   17480   25380                   47100    4700   18100   24300              25   39410    5350   18070   15990                   42100    5900   20700   15500              26   40790    5470   16480   18840                   36400    5000   15700   15700              27   31760    5480   14570   11710                   33700    5900   16200   11600              28   33000    5090   14760   13150                   30900    4800   14600   11500              29   49390    5350   18110   25930                   45500    4900   18100   22400              30   54290    5110   17540   31640                   48800    4800   16700    27300   20900              31   40930    5500   17850   17580                   37300    4900   16700   15700


title: Examples of spatial generalized linear models

text: I've been reading some materials on Spatial data analysis, and I've a good background in GLMs. Right now I'm looking to find an example in spatial generalized linear models, but so far I've not found any. I'm aware of some packages in R such as GeoR. etc. but couldn't find a concrete example that includes codes, data, and analysis. I was wondering if anyone can help. Please share links or papers. Thanks!


title: Find the probability of having k red balls after d days

text: Its an interview question.. Initially there are n white balls. Each day you are allowed to take a ball. If the ball in hand is white, replace it by red ball. If it is red in color, put that ball into bag without doing anything. The question is to find the probability to have k red balls after d days.Give the general solution I started with making a binary tree.. initially on 1st day we have n white balls, when we pick 1 ball and apply the rule we get 1 red ball and (n-1) white balls on second day pick 1 ball, it could be white or red. if we pick white we are left with 2red and (n-2) white balls otherwise 1 red and (n-1) white balls. these two nodes again will each have two children each for the 3rd day.. 1 with 3 red balls and (n-3) white balls , other with 2 red balls and (n-2)white balls and so on But this is a direct recursive formula, is there a better solution? I think it can be solved with dynamic programming, but I am unable to connect dynamic programming with probability. any idea how to do this ? Also can someone help me with how the probability will be calculated in this question? Last, Can someone give me a good resource to study probability based programming questions?


title: Inverse hyperbolic sine transformation: estimation of theta

text: I'm trying to use an inverse hyperbolic sine transformation to reduce the effect of outliers in my target variable. Unfortunately, I don't appear to have access to the basic papers on it. I've found the formulation but am not sure how to estimate the theta parameter for it. Does anyone know? Thanks


title: Variability in k-means clusters results: setting set.seed() before?

text: Then I fit linear models to the `plot(n_clust, error)` aiming to identify the best combination of I'm trying to perform a k-means cluster on my data (matrix with 2000 cases and 10 variables). I don't know how many clusters should I choose. To solve this problem, I adopted a strategy in which different values of K are setted. Error results and the number of used clusters are stored in two vectors. Then I fitted linear models that gives me the highest sum R-squared. I used this strategy to choose de best K value. The point is that I have performed this process several times, noting that error vector varied between each run. This gave me different values of best K. An alternative to make k-means results "stable" is the use of `set.seed()` function prior `kmean()`. However I'm afraid that the result, despite fixed, have no consistency. Somebody could give me some "clues"? `set.seed()` will not just hind a variability? Is there another strategy for choosing the best K? Thanks!               n_clust=NULL     error=NULL     for(i in 1:200){     cl <- kmeans(scores, i, iter.max=100)     erro <- c(error,cl$tot.withinss)     n_clust <- c(n_clust,i)     }          r2=NULL     for(i in 3:197){     a <- lm(error[1:i] ~ n_clust[1:i])     b <- lm(error[i+1:200] ~ n_clust[i+1:200])     rsqd <- as.numeric(summary.lm(a)[8]) + as.numeric(summary.lm(b)[8])     r2 <- c(r2, rsqd)}          id_n <- 3 + which(r2==max(r2))


title: Interpreting the "coefficient" output of the lm function in R

text: I have created a linear model (which has multiple predictors) using the lm() function and I would like to interpret the "coefficients" that I get when I use the summary() function on the linear model.      Now I want to consider how the coefficients reflect on the predictors' influence in the model - am I right in thinking that a large value for the coefficient means that the corresponding predictor has a greater effect? I'm not sure what else I need to consider or if I'm even thinking along the right lines.      Also, am I correct in thinking these "coefficients" are in fact the Beta coefficients?


title: comparing R-squared and F-stat

text: I am doing multiple regression with Gas production (l/d) as response variable, and flow rate (l/h) and COD influent (g/l) as explanatory variables in R. I made two models.               mod1 <- gas ~ flow + COD           r^2 = 0.64          f stat= 56.62(2,64) with 0.0000.... p value.          mod2 <- gas ~ flow*COD          r^2 = 0.81          f stat = 89.14 (3,63) with also very low p value.      Now for both cases, what the r^2 says and the f stat says? I did the anova test of the both models. how to interpret them?               Model 1: Gas ~ Flow + CODin          Model 2: Gas ~ Flow * CODin          Res.Df   RSS Df Sum of Sq      F    Pr(>F)         1     64 24307                                       2     63 12835  1     11472 56.311 2.631e-10 ***


title: How to implement a hold-out validation in R

text: Let's say I'm using the `Sonar` data and I'd like to make a hold-out validation in R. I partitioned the data using the `createFolds` from `caret` package as `folds <\- createFolds(mydata$Class, k=5)`. I would like then to use exactly the fold `mydata[i]` as test data and train a classifier using `mydata[-i]` as train data. My first thought was to use the `train` function, but I couldn't find any support for hold-out validation. I mean, I don't want to partition my data and then test all the combinations for estimating the mean accuracy, I want to train it only once. Am I missing something here? Also, I'd like to be able to use exactly the pre-defined folds as parameter, instead of letting the function partition the data. Does anyone have any thoughts? Thanks in advance


title: Package ‘ﬁtdistrplus’

text: I am trying to use the package ‘ﬁtdistrplus’ in R to fit one non standard distribution to my data set. I am trying to copy the methods the package creators used in their tutorial for specifying a Gumbel distribution but unfortunately I am not able to do it. (The package documentation can be found in this pdf.) So my problem is that I have my CDF function and my density function which are not standard : $ \begin{align*} H_{G}(x)&=\frac{1}{P(1)}\left[ \frac{\Phi^{3}\left( \frac{x-\mu}{\sigma}\right)}{6}-\frac{\delta}{2} \Phi^{2}\left( \frac{x-\mu}{\sigma}\right)+\left( \frac{\delta^{2}}{2} + \beta\right) \Phi\left( \frac{x-\mu}{\sigma}\right)\right]\\ h_{G}(x)&= \frac{1}{\sigma P(1)}\left[ \frac{\left[ \Phi^{2}\left( \frac{x-\mu}{\sigma}\right)\right] }{2}-\delta \Phi\left( \frac{x-\mu}{\sigma}\right)+\left( \frac{\delta^{2}}{2} + \beta\right) \right] \phi\left( \frac{x-\mu}{\sigma}\right) \end{align*}$ with $ P(1)= \left[ \frac{1}{6}-\frac{\delta}{2}+\left( \frac{\delta^{2}}{2} + \beta\right) x\right] $ and \begin{equation*} \Phi(\dfrac{x-\mu}{\sigma}) =\frac{1}{\sqrt{2\pi}}\int^{x}_{-\infty}e^{-\frac{\left(t-\mu\right) ^{2}}{2\sigma^{2}}}dt \qquad \text{and} \qquad \phi(\dfrac{x-\mu}{\sigma}) =\frac{1}{\sqrt{2\pi}}e^{-\frac{\left(x-\mu\right) ^{2}}{2\sigma^{2}}}. \end{equation*} And here is some of my code:               ###### I define my density as in the gumbel ######          ds <- function(u,a,e,mu,sigma){       P = function(a,e){ (( (1/6)*(1^3) ) - ((a/2)*(1^2)) + (((((a)^2)/2) + e)*1)) }       D = function(u,mu,sigma){ (1/((sigma)*sqrt(2*pi)))*exp(-(((u-mu)^2)/(2*(sigma^2)))) }       K = function(u,a,e){ (((1/2)*(u^2))- (a*u) +(((a^2)/2)+e)) }       H = function(u,mu,sigma){ pnorm(u,mu,sigma) }       Fprim = function(u,a,e,mu,sigma){ (1/P(a,e))*(D(u,mu,sigma))*(K(H(u,mu,sigma),a,e)) }       return(Fprim)     }           ###### I define my CDF function as in the gumbel ######          ps <- function(u,a,e,mu,sigma) {       S   = function(u) (((1/6)*(u^3)) - ((a/2)*(u^2)) +                          ( ((((a)^2)/2)+e)*u))/(((1/6)*(1^3)) - ((a/2)*(1^2)) +                          (((((a)^2)/2)+e)*1))       cdf = S(pnorm(u, mu, sigma))       return(cdf)     }          ###### I define my quantile function as in the gumbel ######     qs <- function(u,a,e,mu,sigma) optimize(function(z) (ps(z)-u)^2, c(-10,10))$minimum               ###### The MLE estimation using fitdist and mledist ######     fgu <- fitdist(X, "s", start=list(a=0.035, e=0.005, mu=-0.52, sigma=1))          mledist(X, "s", start=list(a=0.035, e=0.005, mu=-0.52, sigma=1))


title: Trend test for two ordinal variables both containing three factors

text: I have two variables each of which have three levels (low, medium & high). I would like to be able to test whether there is a significant trend between the two variables i.e. as x goes up y goes up. The trend tests I've looked at all seem to have one variable with 2 levels and the other with n levels. Here is an example of the frequency table of my data:               mat <- structure(c(20, 19, 6, 10, 22, 4, 0, 12, 4), .Dim = c(3L, 3L), .Dimnames = list(         c("low_1", "medium_1", "high_1"), c("low_2", "medium_2", "high_2")))      If there was a function for it in R that would be useful.


title: R2 of validation sample

text: After removing 25% (21 observations) of the sample as a holdout, model selection on the original 75% of the a sample led to a six variable multiple linear regression with R2 of 54%. A simple regression of the holdout values on their predictions under the fitted model led to an R2 of 8% and and adjusted R2 of 2%. Regression slope parameter was non-significant. Can anyone think of an interpretation other than that there is potentially a serious overfitting issue here, enough that the original six-variable model is in question? Was there a better way to use the holdout sample to validate the model? Thanks


title: Population with known distribution (not normal) versus small sample

text: I have a complete population consisting of 6587 individuals (genes) and for each individual I have a distance value. These distance values are not normally distributed. I'm interested in determining if a sample of 348 individuals (which has a smaller mean distance value compared with the mean distance value from the complete population) is significantly different than expected. I have limited statistics background and would appreciate advice on what statistical test to use. Ideas include:   1. two sample t-test (probably not right, since its not a normal distribution)   2. wilcoxon rank sum test (buddy suggested this)   3. some other test that's easy to implement in R   4. custom test - one idea I had was to randomly pick out 348 individuals (with replacement) and determine their mean distance. Then repeat this 1000 times (I think this is called bootstrapping). I did this, and I get a normal distribution of values with mean 32000 (same as population mean ) and sd = 2000. The mean distance value for the actual sample of 348 individuals was 18000. So by 1-tail p-value from normal distribution of mean 32000 and sd=2000 this appears significant. Thanks in advance for any advice!


title: Cumulative Distribution Function

text: If I need to determine the value of a data point at each percentile of a distribution (such as lognormal, weibull, etc), can use its CDF and plug in %iles to get the value? For e.g. finding the value at 95%ile by plugging in 0.95 in the CDF formula. Basically, I need to determine that if the distribution graph were drawn what data points do I obtain at each %ile. Does CDF tell me the same thing? Please advise urgently Thanks,


title: bivariate probit with endogenous covariate testing

text: I am interested in learning more about testing for the bivariate probit model with an endogenous treatment regressor. I have figured some stuff out -- summary below, since I don't see much on this topic -- but other questions remain. Here's the setup. Suppose I have a binary outcome $y_1$, which depends on $x_1$, the error $\varepsilon_1$ and binary treatment indicator $y_2$, which itself depends on $x_2$, and $z$ and error $\varepsilon_2$. The dependence is of the form $y_1=\mathbb{1}(\beta_1'x_1+\alpha y_2+\varepsilon_1>0)$ and $y_2=\mathbb{1}(\beta_2'x_2+\gamma z+\varepsilon_2>0).$ The errors int the two equations are correlated with $\rho$. I am interesting in testing the assumptions of:   1. normality of the two errors terms   2. $\rho=0$   3. homoskedasticity of the errors   4. exogeneity and weak instruments for the treatment equation Normality can be tested with goodness-of-fit Rao/Murphy score test described by Chiburis et al. (2011) and Chiburis (2010), who provide Stata code do so (`scoregof`). This test embeds the bivariate normal distribution within a larger family of distributions by adding more parameters to the model and checks whether the additional parameters are all zeros using the score for the additional parameters at the biprobit estimate. It rejects when there is excess kurtosis or skewness in the error distributions. They do not recommend a variant of the Hosmer-Lemeshow test to do this based on simulations. The correlation between the errors can be tested using a likelihood ratio test based on the idea that if $\rho=0$, the log-likelihood for the bivariate probit will be equal to the sum of the log-likelihoods from the 2 univariate probits. If you calculate Huber/White sandwich errors, this becomes a Wald test. Some questions remain.   1. Should I worry about heteroskedasticity if I use robust errors? How can I test for this?   2. Can I differentiate heteroskedasticity from heterogeneous treatments effects?   3. Can I use linear IV diagnostics (like weak instruments tests) to check the non-linear probit model? Is there anything better? This seems strange to me since they estimate different treatment effects.   4. Are there other things I should check that I am not aware of?


title: How to scale new datas when a training set already exists

text: Here is what I have : A scaled training set, with labels. Segmented images, from which I extract new vectors to classify. My classifier is a KNN which would have obviously been trained using my training set. Now, I wonder how I should scale those new vectors I just got. Is this correct to scale them on their own, or should I do something else ? I wonder for example if an outlier would have an effect on the scaling and subsequent classification... [EDIT] adding an outlier (which I would like to detect using kNN algorithm) to the test datas does impact the scaling, so subsequent classification won't work properly. What should I do then ? [EDIT 2] This is how I scale my data : ![enter image description here](http://i.stack.imgur.com/9tpnF.png) Which in Scilab I translate to :               function dataout = scaledata(datain)          dataout = zeros(size(datain,1),size(datain,2));          for i=1:size(datain,2)         dataout(1:size(datain,1),i) = (datain(1:$,i) - min(datain(1:$,i))) / ...                                          (max(datain(1:$,i)) - min(datain(1:$,i)));     end          endfunction      Thank you