title: How to deal with a mix of binary and continuous inputs in neural networks?

text: I'm using the nnet package in R to attempt to build an ANN to predict real estate prices for condos (personal project). I am new to this and don't have a math background so please bare with me. I have input variables that are both binary and continuous. For example some binary variables which were originally yes/no were converted to 1/0 for the neural net. Other variables are continuous like `Sqft`. Sample of input data I have normalized all values to be on a 0-1 scale. Maybe `Bedrooms` and `Bathrooms` shouldn't be normalized since their range is only 0-4? Do these mixed inputs present a problem for the ANN? I've gotten okay results, but upon closer examination the weights the ANN has chosen for certain variables don't seem to make sense. My code is below, any suggestions?               ANN <- nnet(Price ~ Sqft + Bedrooms + Bathrooms + Parking2 + Elevator +                  Central.AC + Terrace + Washer.Dryer + Doorman + Exercise.Room +                  New.York.View,data[1:700,], size=3, maxit=5000, linout=TRUE, decay=.0001)      * * * **UPDATE:** Based on the comments below regarding breaking out the binary inputs into separate fields for each value class, my code now looks like:               ANN <- nnet(Price ~ Sqft + Studio + X1BR + X2BR + X3BR + X4BR + X1Bath             + X2Bath + X3Bath + X4bath + Parking.Yes + Parking.No + Elevator.Yes + Elevator.No              + Central.AC.Yes + Central.AC.No + Terrace.Yes + Terrace.No + Washer.Dryer.Yes              + Washer.Dryer.No + Doorman.Yes + Doorman.No + Exercise.Room.Yes + Exercise.Room.No              + New.York.View.Yes + New.York.View.No + Healtch.Club.Yes + Health.Club.No,         data[1:700,], size=12, maxit=50000, decay=.0001)      The hidden nodes in the above code are 12, but I've tried a range of hidden nodes from 3 to 25 and all give worse results than the original parameters I had above in the original code posted. I've also tried it with linear output = true/false. My guess is that I need to feed the data to nnet in a different way because it's not interpreting the binary input properly. Either that, or I need to give it different parameters. Any ideas?


title: Time series, x-y coordinates, regression, R

text: I have data in the form of these columns:               date, x coordinate, y coordinate, value A, value B, value C, value D, etc.       (I don't see the possibility to copy an Excel table into this text so I show just the headers of the columns) Now I have read and Googled a lot already but I can't find the code to make a regression model in R predict value A (numerical) based on the date, coordinates and the rest of the values (both numerical and text, no negative or complex numbers). Does anyone knows which code in R should work best? I already get nice regressions with the simple `lm()` function if I stick to years as a label instead of dates but then I of course lose the time- connection. For background its maybe handy to know that it is real estate data.


title: Aggregating all measurements per x-value in least-square fitting

text: We want to select one out of a given set of (continuous) functions that best matches a set of observations $\{(s_i,c_i) \mid 1 \leq i \leq N\} \subseteq\mathbb{N} \times \mathbb{N}$ (input size and counter) of unknown distribution. We have many observed counters per size¹. So far, we have used least-square fitting on the _average_ values per size. This obviously neglects any asymmetry in the data distribution (for a fixed size). Consider this example: ![enter image description here](http://i.stack.imgur.com/E4y2u.png) The blue blobs are the actually observed counters (one desaturated fleck per measurement, creating something like a violin plot), the red dots are the averages (per size) and the blue line is our estimation. Apparently, the actual measurements have some asymmetric distribution around the average. Now we think that we can improve the result by fitting against _all_ observations. This would increase the amount of data to be stored, transmitted and computed with a lot, so we would like to get around this. So this is our question: > Is there a statistic that can replace all observations per size without > changing the result of least-square fitting, i.e. is there an $S : > \mathbb{N}^\mathbb{N} \to \mathbb{N}$ such that > > $\qquad \displaystyle \arg\min \sum_{i=0}^N (c_i - f(s_i))^2 \ =\ \arg\min > \sum_{i=0}^{\max_j s_j} (S(C_i) - f(i))^2$ > > with $C_i = \{c_j \mid s_j = i\}$ for all (reasonable) functions $f$? We would also be fine with using another criterion but least-square on the right hand side, as long as it is (about) equivalent (given $S$). Also, approximate equality is fine as long as we have equality in the limit (for $N$ and/or $\max_j s_j$ to $\infty$). * * *   1. We count how often a given algorithm hits a given basic block when executed on inputs of certain sizes. For an average case estimation, we run the algorithm on random sample of inputs per size.


title: Why is semipartial correlation cited so seldom?

text: In OLS regression, I find the semipartial correlation (a.k.a. part correlation) to be a very useful indicator. When squared, it shows each predictor's unique contribution to explained variance in the outcome. If one tries several models, entering predictors in various orders, the semipartial correlation will always turn out the same by the final model. This indicator can also be a useful basis for creating Venn diagrams showing the relative importance of different predictors. Given all this, why is it cited so seldom? For example, after 6,300 questions on this site, it has never come up, nor do I see it cited in journal publications.


title: Likelihood ratio tests

text: Is likelihood ratio test ($F$-test) of significance of difference of two linear models the same as chi-square test of difference of $-2\log L$? SAS `PROC GLM` produces $F$-statistics and `PROC MIXED` $-2\log L$.


title: How to calculate a confidence interval for Spearman's rank correlation?

text: Wikipedia has a Fisher transform of the Spearman rank correlation to an approximate z-score. Perhaps that z-score is the difference from null hypothesis (rank correlation 0)? This page has the following example:               4, 10, 3, 1, 9, 2, 6, 7, 8, 5     5, 8, 6, 2, 10, 3, 9, 4, 7, 1     rank correlation 0.684848     "95% CI for rho (Fisher's z transformed)= 0.097085 to 0.918443"      How do they use the Fisher transform to get the 95% confidence interval?


title: How to interpret coefficients from a GARCH model?

text: Recently I have opened a question here to understand the output of a GARCH model. My goal is to understand if the series I'm checking is heteroscedastic or not. I'm using the `garch()` function from the tseries package. First I built a linear regression like this:               mod <- lm(a ~ b)      Then I need to check if the residuals of this linear regression present heteroscedasticity. I did:               g <- garch(resid(mod), order(c(1,1)))      and then               summary(g)      I get the follow output:               > summary(g)           Call:      garch(x = lm(A ~ B)$resi, order = order(c(1, 1)))           Model:      GARCH(1,2)           Residuals:          Min      1Q  Median      3Q    Max      -4.2058 -1.0262  0.1404  1.1069  3.6553           Coefficient(s):          Estimate  Std. Error  t value Pr(>|t|)         a0 3.361e-04  9.352e-05    3.594 0.000326 ***      a1 3.045e-01  4.486e-02    6.787 1.14e-11 ***      a2 1.209e-06  8.855e-02    0.000 0.999989         b1 4.938e-01  1.060e-01    4.660 3.17e-06 ***      ---      Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1           Diagnostic Tests:          Jarque Bera Test           data:  Residuals      X-squared = 18.84, df = 2, p-value = 8.108e-05                    Box-Ljung test           data:  Squared.Residuals      X-squared = 49.7251, df = 1, p-value = 1.769e-12       Thanks to the user who answered my question, I now understand that the `ao` is the intercept and the other `a1` and `b1` are the coefficients I need to check to understand if this time series is heteroscedastic or not. The problem (doubt) is that now I also see `a2` in the regression table: What does it stand for? Is it correct to say that if all coefficients have a $p$-value above 0.05 (...)


title: How do I calculate projected figures for the next year based on past performance?

text: I'm a SQL/C++ developer who recently has been asked to generate a report from our database to predict some future performance based on historical data; the problem is that I don't have much experience of this sort of data modelling. I initially thought I could take an average of each month's results and use that but after reading articles on the web about statistical models, estimation, forecasting etc. I felt that might not be sufficient. Plus, most of the formulas are over my head as I've never learnt statistics. What process would you recommend as the best way (for me) to calculate these predictions? Preferable something I can translate in to SQL or use in a spreadsheet (Excel)? **Update to link in comment** The link to the document "Statistical flaws in Excel" by Hans Pottel (www.coventry.ac.uk/ec/~nhunt/pottel.pdf) jthetzel has in the comment below no longer exists. I managed to find the document here, which can be downloaded by right-clicking and selecting "save as". _Please note the button to download the file took me to a page blocked by the college as possibly dangerous, hence I downloaded the file using the right- click method._


title: Conditions for Poisson approximation of the superposition of non-Poisson processes

text: It is well known that the superposition of $N$ Poisson processes is itself a Poisson process with an intensity given by $\sum_{n=1}^{N} \lambda _{n}$. Conversely a superposition including any non-Poisson component processes is not Poisson. However, the superposition of sparse processes converges to a Poisson process as $N\rightarrow \infty $ and as the sparseness increases. This is referred to on Wikipedia as the Palm-Khintchine theorem: a recent paper with discussion of the original work is given by Schuhmacher 2005. If I have a collection of sampled processes, the superposition of which is well modeled as a Poisson process, under what conditions of $N$ & $\hat{\lambda_{n}}$ can I conclude that all component processes are Poisson? P.S. This problem is known and has been studied in queuing theory. Newell 1984 provides a condition that if the component arrival processes satisfy $n(1-\rho)^2 {\gg } 1$, where $\rho$ is the utilisation, then the queue will approximate a M/M/1 queue (Poisson arrival process). How does this formula apply to the more general problem outlined above?


title: Simple question about multivariate/multiclass classification

text: From this link Text Classiﬁcation using Naive Bayes, there are two models described for classification, Naive and Bernoulli. My question is if i want to make this classifiers for multiclass (multivariate) classification i should use techniques like: one vs. one/rest and etc?


title: Can I select IV from their t value and confidence intervals in LMM in lmer?

text: As the issue of p value is still inconclusive and many papers point out that the confidence intervals are more informative in selecting factors in fixed effects, I would like to ask about interpretation of linear mixed model in R. Is it possible to select the important factors from the t value comparing to the range of confidence intervals? For example, in the information of fixed effects after running lmer, if the t value equals to -1.77 while the range of 95% CI is between -19.06 to -2.73. It means that this factor is significant to our data because it falls outside CI. What do you think? Cheers, Pat


title: How do I generate numbers according to a Robust Soliton distribution?

text: I'm working on project in Matlab which aim is to demonstrate how Luby Transform codes work. I need to generate generation matrix and I need to get some values from Robust Soliton Distribution, can someone explain me how to generate numbers according to this distribution, I need to generate one value for each column in my matrix, I saw Robust Solition formula here http://en.wikipedia.org/wiki/Soliton_distribution#Robust_distribution but I don't understand what each variable means. Thanks in advance Here is the formula form Wiki: ![enter image description here](http://i.stack.imgur.com/p2bog.png) ![enter image description here](http://i.stack.imgur.com/Hm9xh.png) ![enter image description here](http://i.stack.imgur.com/iDQIx.png)


title: Two-mode weighted network edge weight normalization

text: I have a two-mode network which looks into relationship between forum threads. Two forum threads are connected if the same user has replied to both. The nodes are sized based on number of people who have only replied to the given forum thread, and no other ones (it would be the users who are unique to the given thread) and edges are the number of people who were active in both of the threads. However, there are multiple problems with not doing anything with the edge weight: 1.) If there are two threads that are very active and have 500 users replying to them, and they have an edge of 100 between each other (100 common users), it shouldn't be as significant as if two less active forums share 100 common users; 2.) How do I solve the problem of a very active and a less active forum thread being connected, if at all, without using parallel edges? What would be a good way to calculate significance of edge weight?


title: which neural network is better

text: This is a homework question where I have not been able to reach any conclusion . I have an exam tomorrow . Please help me out . > We have a set of data from patients who have visited a hospital . A set of > features ( eg temperature , height ) have been also extracted for each > patient . Our goal is to decide whether a new visiting patient some $n$ > number of diseases or not . > > This problem is to be solved using neural network . We have two choices : > either to train a separate neural network for each of the diseases or to > train a single neural network with one output neuron for each disease, but > with a shared hidden layer . Which method do you prefer ? Justify your > answer . Thanks in advance . My opinion : Both methods are almost same because we have to train the same number of weights in both the case . But in one case we can train them all at once and in another we will train them individually so for training it will take more time .


title: Constant signs of correlation in the active set in least angle regression

text: I am trying to comprehend the proof of the Least Angle Regression algorithm and I am stuck at certain points. I would appreciate any help that I can get. Let me set the stage: I am following the paper in the following link On Page 8, in the paper, they make the following claim (restated in my own words). > The sign of the correlation of every predictor in the active set stays a > constant. Though I intuitively understand it I am trying to formally prove this claim. I have some ideas, but none give me conclusive proof. For example if I end up proving that the step taken in the equi-angular direction is less than the least squares step then I am done. But I am not able to prove that either. Can anyone help me?


title: Probability of a random variable to be the largest among others

text: Let us have $N$ random variables generated by uniform distribution. That is, $$u_i \sim \mathcal{U}(0,1),\quad i=1,\ldots,N$$. What is the probability of $u_N$ being the largest? I.e., how can I compute $$p\Big(u_N\geq \max(u_1,u_2,\ldots,u_{N-1})\Big)$$


title: Geometric Interpretation of Maximum Likelihood in a location model

text: Can you please check my work in the problem below? * * * Consider a location model: $$X_i =\theta +e_i, \quad i=1,\ldots,n $$ where $e_1, e_2, \ldots e_n $ are iid with pdf $f(z)$. Define $\mathbf{\mu}=\theta \mathbf{1}$ where $\mathbf{1}$ is a vector with all its components equal to $1$. Let $V$ be the subspace of vectors of the form $\mathbf{\mu}$. Then $$\mathbf{X}=\mathbf{\mu}+\mathbf{e} \quad \mathbf{\mu} \in V $$ It makes sense therefore to estimate $\mathbf{\mu}$ by a vector in $V$ which is "closest" to $\mathbf{X} $, i.e. $$\mathbf{\hat{\mu}}=\mathrm{Argmin} || \mathbf{X}-\mathbf{v} ||,\quad \mathbf{v} \in V $$ If the error PDF is the standard Laplace, show that the minimization of the above is equivalent to maximizing the likelihood when the norm is the $l1$ norm given by $$||v||_1=\sum_{i=1}^n | v_i | $$ * * * Of all the components of the location model, the error is the random one so it would make sense to make the one-to-one transformation $e_i=x_i-\theta$ in the PDF of $e_i$. The likelihood then becomes: $$L \left( \theta \right) =\frac{1}{2^n} exp \{ -\sum |x_i -\theta | \} $$ and the log likelihood $$ l \left( \theta \right) =-n \log2 -\sum_{i=1}^n |x_i -\theta | $$ Now it seems very clear that the maximization of $ -\sum_{i=1}^n |x_i -\theta |$ is equivalent to the minimization of $\sum_{i=1}^n |X_i-\theta| $ as desired. Does that suffice or am I required to do more? Thanks in advance.


title: Predicting the difference between two groups in R

text: I'm working on practice exercises with R. From the data, I now have to determine the difference in predicted male compared to female spending while keeping all other predictors constant. I replaced the sex `0,1` with `male, female`:               > spending$sex=factor(spending$sex)     > levels(spending$sex)=c("male", "female")     > tail(spending)         sex status income verbal gamble     42 male     61  15.00      9   69.7     43 male     75   3.00      8   13.3     44 male     66   3.25      9    0.6     45 male     62   4.94      6   38.0     46 male     71   1.50      7   14.4     47 male     71   2.50      9   19.2     > summary(spending)          sex         status          income           verbal           male  :28   Min.   :18.00   Min.   : 0.600   Min.   : 1.00        female:19   1st Qu.:28.00   1st Qu.: 2.000   1st Qu.: 6.00                    Median :43.00   Median : 3.250   Median : 7.00                    Mean   :45.23   Mean   : 4.642   Mean   : 6.66                    3rd Qu.:61.50   3rd Qu.: 6.210   3rd Qu.: 8.00                    Max.   :75.00   Max.   :15.000   Max.   :10.00            gamble           Min.   :  0.0        1st Qu.:  1.1        Median :  6.0        Mean   : 19.3        3rd Qu.: 19.4        Max.   :156.0        Now, this is where I got stumped, how do I set this up with constants and only use the `sex` variable?


title: How to estimate an upperbound for logistic regression by only 5 to 7 data points?

text: I have data that is of the form $y = \frac{\beta_1}{1 + \exp(\beta_2 + \beta_3 * x)}$. For the estimation of $\beta_1$ to $\beta_3$ I use the formulas of this paper: John Fox - Nonlinear Regression and Nonlinear Least Squares In this paper, $\beta_1$ is estimated by looking at the data. If I do that it works fine, even if I have only three points. From that I can calculate the two others. I tested my parameters with nls() in R and LevenbergMarquardt in C#. The models returned by them are satisfying. The problem is I don't want to look at the data to get a good estimator for $\beta_1$, I want my program to calculate it. For some time I used values that were a bit higher than the maximum of my values (something between $\max * 1.1$ and $\max * 1.5$. This worked fine as long as the points covered most of the function. It also worked fine if the data points were somewhere from the "top" of the curve, but when they all came from the area "below" the inflection point, this estimator was definitely lower than it should be and I couldn't fit the model. If I use something that is definitely higher than the maximum point (by multiplying it with ridiculously high values) the model doesn't fit in any useful way. Measurments could look like this: x = (40, 50, 60, 70), y = (1000, 950, 400, 200) -> easy to estimate x = (40, 50, 60, 70), y = (1000, 950, 800, 100) -> easy to estimate x = (40, 50, 60, 70), y = (500, 200, 100, 50) -> not so easy to estimate I guess I could find out where I am in the function (at the "bottom", at the "top", in the slope) by calculating the deltas in the given points and calculate an upper bound depending on that. **Does anyone have a hint for a better solution?** Additional info: If it can't be done, it's more important for me that measurements that can be fit are as good as possible and I accept that some measurements can't be fit at all. (although I want an implementation in C# I posted it here, I don't think the problem is language dependent) **update (Application of this):** x are temperature values und y the corresponding measurements. By nature it is given that it should look like a logistic curve with higher y values at lower temperatures and vice versa. The melting point is equal to the inflection point of the curve, which changes a lot with small changes of the model parameters. **update (some made up data with 7 data points and known inflection point at 60):**               //first I made up some data without any noise and enough (13) points     double[] x17 = { 30, 35, 40, 45, 50, 55, 60, 65, 70, 75, 80, 85, 90 };     double[] y17 = { 2000, 1920, 1820, 1700, 1500, 1280, 1000, 720, 500, 300, 180, 80, 0 };     //calculated inflection point: 59.642....          //Then I took three different parts of that data      (to simulate how much data I will have in the real world)     double[] x18 = { 30, 35, 40, 45, 50, 55, 60 };     double[] y18 = { 2000, 1920, 1820, 1700, 1500, 1280, 1000 };     //calculated inflection point: 59.263.... is ok!          double[] x19 = { 60, 65, 70, 75, 80, 85, 90 };     double[] y19 = { 1000, 720, 500, 300, 180, 80, 0 };     //calculated inflection point: 53.447.... to small!          double[] x20 = { 45, 50, 55, 60, 65, 70, 75 };     double[] y20 = { 1700, 1500, 1280, 1000, 720, 500, 300 };     //calculated inflection point: 59.868... almost perfect!      My current approach for an upper bound estimation is y0 * y0 / y1. This way I hope to take into account that I might not always have a value close to the maximum. **update:** The x and y values will never be negative. x will always be { 40, 45, 50, 55, 60, 65, 70 } unless a data point is lost. **update:** I did many, many tests with simulated data that should be easy to fit (I choose data points that were exactly on the curve = without any noise) and I see that it works fine, unless the first or second data point (where x=40 or x=45) is missing. I guess I will have to discard such measurements and the user will have to live with that.


title: Measures of internal consistency of survey items

text: I have a question on a survey that asks students how useful a feature of the software was that they used for learning. I then have three questions that assess in what way the feature was useful. The measures use a 7-point Likert scale from Very Strongly Disagree to Very Strongly Agree:               The [[feature]]     ...was useful for learning overall                               ...helped me to mentally organize the structure of the content.                              ...enabled me to control my learning pace                                ...made it easy for me to find the content I needed      The question I'm trying to answer is "If the feature was useful, in what way was it useful?" I'm wondering what kind of statistical tools can be used to do this. Could a covariance matrix determine this sufficiently? Could I compute Cronbach's alpha for each question paired with the first question? I have also come across Confirmatory Factor Analysis, but I don't yet understand it well enough to know if it is an applicable tool. I'm new to stats so I'm just looking for a pointer in the right direction. Edit: If it makes a difference, it is a within-subjects experiment where the student has used the software both with and without the feature.


title: How can I change the scores to make interaction effect in two way ANOVA significant?

text: My question is how to make an interaction effect become significant in two-way ANOVA. If my main effect A has an F value of F(1,80)=25.5 for example, and my main effect B has an F value of F(3,80)=15, but my interaction effect has a value of F(3,80)= 2.06 which is not significant, how do I change the score in my SPSS data file to make the interaction effect become significant so that I can run further interaction contrasts. Any help would be appreciated.


title: Find probability and expectation

text: 18 boys and 2 girls are made to stand in a line in a random order.Let $X$ be the number of boys standing in between the girls .Find $P[X=5]$ and $E[X]$. I proceed in this way: Note that $$P[X=0]=\frac{19.18!.2!}{20!}=\frac{19-0}{20\choose2}$$ $$P[X=1]=\frac{18.18!.2!}{20!}=\frac{19-1}{20\choose2}$$$$P[X=2]=\frac{17.18!.2!}{20!}=\frac{19-2}{20\choose2}$$So, $$P[X=x]=\frac{19-x}{20\choose2}$$,$x=0(1)18.$ so.$$P[X=5]=\frac{19-5}{20\choose2}=0.0737$$ and $$E[X]=\sum_{x=0}^{18}x.\frac{19-x}{20\choose2}=6$$


title: How to tell how extreme an outlier is?

text: I am analyzing some data and want to look at one particular point and see how "extreme" it is. Do I exclude this outlier from the data, calculate the dataset's standard deviation and average, then compare my outlier to THAT, or do I calculate the standard deviation and average WITH the outlier included and THEN analyze my outlier with respect to those metrics?


title: What to do when weak classifiers are almost identical in AdaBoost?

text: I have written some code which uses AdaBoost to generate a set of weak classifiers. I'm finding, however, that when I use the resulting strong classifier to classify examples, it seems that almost every weak classifier classifies an example the same way. Either they all think it's negative, or they all think it's positive. As a result, instead of having a nice confidence measure at the end of the classification, I get an essentially binary value (either positive and > 1 or negative and < -1), and nothing really in between. I'm using a weighted linear least squares weak classifier. When I examine the solved coefficients, I notice that they're very, very close in value among all of the weak classifiers. This seems strange to me, because I would expect that the AdaBoost method of weighting the difficult examples would cause the weighted linear least squares results to vary dramatically to 'fix' the difficult examples. What is this telling me about my training set? Is it saying that the data are not linearly separable? Could there be a problem with my weighting scheme? Note also that there may be outliers in the positive and negative example set. Could this be a root cause? I'm hoping somebody could provide some insight into how it could be that the weighted least squares solutions all seem to be almost identical. Thanks!


title: Computing a p-value using bootstrap

text: I have a relatively complicated model-fitting procedure, the output of which is a set of coefficients $\beta$ (one for each feature in the training set) and a set of hyperparameters $\lambda$ (which control the amount of regularization). I would like to compute p-values for each feature. At the moment, my procedure is this:   1. Determine good values for the hyperparameters $\lambda=\lambda^*$ using cross-validation.   2. Fit the model on the whole data set with these hyperparameters, obtaining $\beta=\beta^*$   3. For each feature $i$      * For $j=1$ to $n$        * Shuffle the values of feature $i$ to produce a new data set $X^j$       * Fit the model to $X^j$ with the hyperparameters $\lambda$ to obtain $\beta^j$     * Count the proportion of times that $|\beta^j_i| > |\beta^*_i|$ and call this the p-value for feature $i$ This matches my intuitive understanding of what a p-value is (i.e. each of the reshufflings produces a data set where the null hypotheses $\beta_i=0$ holds, since I have randomized away all the dependence) but I wondered how 'canonical' it is. I can think of a couple of other ways to do it, for example:   1. Refit the hyperparameters as well as the coefficients on each inner loop.   2. Use sampling with replacement to generate the new data set, instead of shuffling. Are either of those better? And, more importantly, can the numbers I'm generating be interpreted as p-values with any degree of reliability?


title: How can I create nice graphs automatically?

text: Eg. such as those on this page http://store.steampowered.com/hwsurvey Is there any ready-made software that can do this? Alternatively, any recommendations for other software that does something similar? I know this is not really a statistics question, but I feel very strongly that data should be presented in a neat and appealing way in order to be effective so I think this question would be uesful to some people. **Update (29/12/11)** : Thanks for all the responses to this question, I really appreciate all your suggestions. By automatically, I mean that I key in the data & the graph gets automatically updated. The purpose of my project is to collect 6-10 sets of data daily (or perhaps 2x daily) over an indefinite period of time, and I want to find a way to represent the data in a website (similar to the Steam link that I provided above), while the user backend is simple enough for a non-technical user to handle. If you have any more suggestions, please add to the answers! Thanks again!


title: Is there a standard process for finding linear relationships with matching?

text: Most of my reference on matching has been Rosenbaum's book "Observational Studies". From what I understand, the usual process of matching involves administering a treatment that is represented by a binary value of either 0 or 1, depending on whether the patient got the treatment or not. In the case that I am looking at, the treatment is represented as a continuous value with no particular meaning given to "0" treatment, making it difficult/non-applicable to assign a control. In my hypothesis there is a linear relationship between the amount of treatment and the outcome and I want to test for that. During these tests I need to apply all the corrections for overt biases included in the Matching procedure. What is the best method for testing for a linear relationship? What I have been doing, perhaps naively, is sampling at 2 different levels of treatment, $x_1$ and $x_2$, and with their mean outcomes, $u_1$ and $u_2$, I've been estimating the relationship as the line given by the two-point equation passing through $(x_1, u_1)$ and $(x_2, u_2)$. This turns out to be too inaccurate. Would it be better for me to sample at several levels of treatment $x_1, x_2, x_3, \dots,x_n$ and do a linear regression on that data or is there another standard way of doing it? Thanks


title: SAS for regression with categorical and quantitative explanatory variables

text: I am analyzing growth over time for 5 different cultivated forms (cultivars) of maize. Graphing the data reveals a clear linear pattern for all the cultivars in the time interval I am interested in. Multiple linear regression with categorical (5 cultivars) and continuous (7 time points) explanatory variables appears to be one way to approach this problem, but I am having trouble with the coding in SAS 9.3. I'd like to compare slopes and intercepts between all the genotypes. I should add that new (independent) organ samples were taken at each time point (e.i. I did not measure the same plant organ at each time point). So far, I've read the `proc reg` manual, and some very useful information at http://www.ats.ucla.edu/stat/sas/webbooks/reg/Chapter7/sasreg7.htm. Nonetheless, I haven't found a straightforward example on how to approach this problem. Any suggestions?


title: Linear regression for feature selection

text: Imagine we regress `y` on `x1`...`x4`. Now, we want to find out if `x5` is a stronger predictor than `x6` (given the other variables). Note that all variables are scaled. Would it be okay to use the residuals to see which one would be a stronger predictor?               y <- scale(rnorm(1000))     x <- scale(replicate(6, rnorm(1000)))          # Method 1:     res = lm(y ~ x[,1:4])$residuals      lm(res ~ x[,5] - 1)     lm(res ~ x[,6] - 1)      The goal here is to identify which variable is a stronger predictor (taking into account the other variables). As far as I can see, this indeed delivers different results from simply correlating `x5` and `x6` with `y` ( _method 2_ ) in turn. The benefit of doing it this way is that it would be less computationally expensive (with high amount of predictors) to compute rather than computing the whole equation. Also, the results still differ a bit from when we would compute them all at once, that is `lm(y ~ x[,1:5])` and `lm(y ~ x[,c(1:4,6)])` separately ( _method 3_ ).               results                   x5           x6     explain residuals    -0.003126777 -0.008349196     cor(x[,5:6], y)      -0.003499607 -0.006773532     explain at once      -0.003137124 -0.008407007      So: is there any kind of a shortcut that could produce the latter model without having to compute the large model? What would be the advice for feature selection? Is explaining the residuals a good approximation of how good the model would be including `x5` or `x6` from the start? Added some benchmark results (`10000x1002` matrix):                             x1001      x1002     time taken     method1     -0.01515   -0.00967       16s       method2     -0.01690   -0.01170    0.001s     method3     -0.01689   -0.01068       32s      This might actually suggest that `cor()` might be good enough, or does this have to do with the fact that here all `x`'s are independent of each other, while in reality this is most likely not the case?


title: Is it possible to conduct a regression if all variables are ordinal?

text: I would like to produce a regression analysis model. I have ordinal categorical data. I can use SPSS. I do not know what analysis to perform or what assumptions to check. What are the statistical tests that can be performed? How can I test the accuracy of the model?


title: What is the proper name for a "river plot" visualisation

text: In a famous plot, Charles Minard visualised the losses of the French Army in the Russian campaign of Napoleon: ![enter image description here](http://i.stack.imgur.com/A7z2Q.png) (another nice example is this xkcd plot) Is there a canonical name for this type of visualisation? I'm actually looking for an R package to create such plots, but I don't even know how to look for it. **EDIT:** As I could not find a good package in R do do this type of plots, I have created my own, called "riverplot" -- you can download it from CRAN. Here is a simplified version of the above diagram: ![enter image description here](http://i.stack.imgur.com/DZQYK.jpg) And an example of what other diagrams can be created with the package: ![enter image description here](http://i.stack.imgur.com/DV80v.jpg)


title: Significant regression for the dichotomised predictor variable and not the continuous version of a predictor variable?

text: I am interested in examining the association between cardiorespiratory fitness and PCA derived dietary pattern scores in 264 adolescents. My cardiorespiratory fitness predictor is vo2max, and my outcomes are a "Treat Foods" pattern and a "Fruit and Vegetable" pattern. Both vo2max and my pattern scores are continuous variables. I have also dichotomised vo2max based on an adequate health cutoff (the FITNESSGRAM Healthy Zones). The sample is relatively fit, so the categories are unbalanced with 88% categorised as having adequate cardiorespiratory fitness. I am using multivariate regression analysis in STATA 12IC, and my models appear to be linear. What I have found: 1) There is an association between vo2max and the "Fruit and Vegetable" score, but not between the FITNESSGRAM variable and the "Fruit and Vegetable" score. My assumption is that this may be due to a loss of power with the dichotomising. 2) While there is no association between vo2max and the "Treat Foods" pattern, there is an association between the FITNESSGRAM variable and the "Treat Foods" pattern. Can anybody suggest a statistical reason for this discrepancy, or a better way of handling the data, for the second result in particular?


title: Relation between Scoring rule and Loss function in Parameter estimation and model selection?

text: Initially, I had only heard of MLE and use it for almost everything, e.g. point estimate and model selection (with some penalty). Then, MSE appeared, which seems to play the same role as MLE does. I learned that both MSE and MLE are results from loss functions (entropy vs quadratic). Choose which one to use is completely up to one's belief. Now, I came across stuff called scoring rule. There are new terms like Logarithmic scoring function and Brier scoring function. However, I feel they are just another names of MLE and MSE methodology. Minimize Brier is essentially minimizing MSE, and Maximize Logarithmic is exactly the same as MLE. Am I right here? Why do we have different names for the exact same things between scoring rule and loss function?


title: Finding the MSE of the MLE of N?

text: If X is taken from a uniform discrete distribution, how can I find the MSE of the MLE of N? (I got disturbing and contradictory answers). I am having some trouble with calculating the mean squared error of the MLE and it seems that I am getting two answers depending on which definition I am using. I found that the MLE of N is just $N=X$, then I used the formula for mean- squared error to get: $MSE(N) = Bias^2(N) + Var(N)$ and $Bias(N) = E(N) -N = \frac{N+1}{2}-N = \frac{1-N}{2}$ and $Var(N) = \frac{N^2-1}{12}$ I ended up getting $MSE(N) = \frac{2N^2-3N+1}{6}$. However, using the definition of bias with regards to the expectation of the squared error I got: $E(N-X)^2 = \frac{(2 N^3-3 N^2+N)}{6}$ Why did I end up with an extra factor of N? Thanks!


title: Variance of sample median

text: Given a sample of n units out of a population of N, population median can be estimated by the sample median. How can we get the variance of this estimator?


title: Hypothesis testing with the geometric distribution

text: > A single observation of a random variable having a geometric distribution is > used to test the null hypothesis $\theta=\theta_0$ against the alternative > hypothesis $\theta=\theta_1 > \theta_0$. > > If the null hypothesis is rejected if and only if the observed value of the > random variable is greater than or equal to the positive integer $k$ , find > expressions for the probabilities of type I and type II errors. My idea: Type I error equals $\alpha$, which is the probability that we reject $H_0$ while it is true. Probability of type II error is $\beta$, which is the probability that we accept $H_0$ while $H_0$ is false. Let's say the Null Hypothetis is rejected, this means that $X\ge k$ ? I know that that $P(X\ge k ; \theta_0)= 1- \theta_0(1-\theta_0)^{k-1}$. Is this value equal to $\alpha$? How can I calculate $\beta$? $P(X\ge k ; \theta_1)= \theta_1(1-\theta_1)^{k-1}$. Who can help me making this clear ? :-) Thank you for your trouble :-)


title: Steps after cross-validation in Linear Regression

text: After making some variable transformations, I chose this model because I works well with the assumptions and gives the largest R-Squared and the smallest Overall MS in cross-validation between a set of linear models.               > gn<- lm(NA.~ I(PC^0.25) + I(((PI)^2)),data=DSET)     > summary(gn)          Call:     lm(formula = NA. ~ I(PC^0.25) + I(((PI)^2)), data = DSET)          Residuals:         Min      1Q  Median      3Q     Max      -425.22  -87.46   -2.30   79.11  396.14           Coefficients:                   Estimate Std. Error t value Pr(>|t|)         (Intercept) -2.047e+03  1.094e+02  -18.71  < 2e-16 ***     I(PC^0.25)   1.206e+03  4.231e+01   28.52  < 2e-16 ***     I(((PI)^2)) -5.242e-02  1.233e-02   -4.25 2.81e-05 ***     ---     Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1          Residual standard error: 127.2 on 319 degrees of freedom     Multiple R-squared:  0.7475,    Adjusted R-squared:  0.746      F-statistic: 472.3 on 2 and 319 DF,  p-value: < 2.2e-16      The output of the cross-validation is the following:               > a<-CVlm(df=DSET,form.lm = gn ,m=5)     Analysis of Variance Table          Response: NA.                  Df   Sum Sq  Mean Sq F value  Pr(>F)         I(PC^0.25)    1 14979750 14979750   926.5 < 2e-16 ***     I(((PI)^2))   1   292051   292051    18.1 2.8e-05 ***     Residuals   319  5157489    16168                         ---     Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1                    fold 1      Observations in test set: 64                Sum of squares = 1119289    Mean square = 17489    n = 64           fold 2      Observations in test set: 65           Sum of squares = 1072687    Mean square = 16503    n = 65           fold 3      Observations in test set: 65           Sum of squares = 865598    Mean square = 13317    n = 65           fold 4      Observations in test set: 64                Sum of squares = 1178545    Mean square = 18415    n = 64           fold 5      Observations in test set: 64           Sum of squares = 1042335    Mean square = 16286    n = 64           Overall (Sum over all 64 folds)         ms      16393                >       So, my question is, Which steps do you recommend to me after the cross- validation? Should I choose the training set applied to the fold 3(it gives the small mean square error)? Is there a way to improve the coefficients after the cross-validation?


title: Nonlinear optimization with constraints in R

text: I'm trying to calculate the optimal portfolio regarding Foster-Hart riskiness measure satisfying following equation: 1) E(log10(1+v/R))=0 where: R is the value of riskiness v is the vector of daily returns of the portfolio. v is composed as a sum of products of individual investment returns by portfolio weights: v = [w1....w10] %*% t(daily_change) My approach was to use the ALABAMA package (http://cran.r-project.org/web/packages/alabama/alabama.pdf) and using a vector of portfolio weights and R value as an input ([w1,....,w10,R]) and minimizing a function giving a value of R under the constraints of equation 1), sum(w1:w10)=1 and [w1:w10]>=rep(0,10). Unfortunately, the algorithm doesn't handle the problem well - it works if a constraint on R close to the minimum is set, otherwise it either gives an error or doesn't change the constrained value of R. Do you have any ideas how to approach this problem? The R value for a given portfolio can be easily calculated with the Newton method, but I couldn't figure out how to connect it with picking the optimal portfolio weights. I'll be extremely happy with any help, I'm in despair... I will gladly provide the R codes.


title: Using mutual information to estimate correlation between a continuous variable and a categorical variable

text: As for the title, the idea is to use mutual information, here and after MI, to estimate "correlation" (defined as "how much I know about A when I know B") between a continuous variable and a categorical variable. I will tell you my thoughts on the matter in a moment, but before I advice you to read this other question/answer on CrossValidated as it contains some useful information. Now, because we cannot integrate over a categorical variable we need to discretize the continuous one. This can be done quite easily in R, which is the language I have done most of my analyses with. I preferred to use the `cut` function, since it also alias the values, but other options are available as well. The point is, one has to decide _a priori_ the number of "bins" (discrete states) before any discretization can be done. The main problem, however, is another one: MI ranges from 0 to ∞, as it is an unstandardised measure which unit is the bit. That makes very difficult to use it as a correlation coefficient. This can be partly solved using **global correlation coefficient** , here and after GCC, which is a standardized version of MI; GCC is defined as follow: ![enter image description here](http://i.stack.imgur.com/0U3Q8.png) _Reference: the formula is from Mutual Information as a Nonlinear Tool for Analyzing Stock Market Globalization by Andreia Dionísio, Rui Menezes & Diana Mendes, 2010._ GCC ranges from 0 to 1, and therefore can easily be used to estimate the correlation between two variables. Problem solved, right? Well, kind of. Because all this process depends heavily on the number of 'bins' we decided to use during the discretization. Here the results of my experiments: ![enter image description here](http://i.stack.imgur.com/22WmE.png) On the y-axis you have GCC and on the x-axis you have the number of 'bins' I decided to use for discretization. The two lines refers to two different analyses I conducted on two different (although very similar) datasets. It is possible to find the complete procedure plus the code used here. It seems to me that the usage of MI in general and GCC in particular is still controversial. Yet, this confusion may be the result of a mistake from my side. Either the case, I'd love to hear your opinion on the matter (also, do you have alternative methods to estimate correlation between a categorical variable and a continuous one?).


title: Does RandomForest ignore spatial independence?

text: I have 5 variables for each countries of the world and I need to analyze their effect and interactions on an independent variable. Random Forest would be adequate for my scope as it deals with non-linear relationships and predicts variables importance. However, I'm wondering if spatial dependence may be an issue. I've never seen spatial dependence discussed in RF applications even if it has been widely used for spatial data.


title: Interpretation of Log-level difference-in-differences specification

text: When I run a standard Difference in differences specification with a log- transformed dependent variable, like: $log(Outcome_{it}) = \beta_1 + \beta_2*Treat_i +\beta_3*Post_t +\beta_4(Treat*Post)_{it} +\epsilon_{it}$ How do I interpret the coefficient $\beta4$ ? Normally in log-level models I would use the approximation $\%\Delta y= 100*\beta_4*\Delta x $ But this approximation is only valid for small changes in $x$ (and small $\beta$). In my case $x$ is a dummy variable and as such either $0$ or $1$ ( Treat*Post). Is this change considered to be small or do I have too use $\%\Delta y=e^( \beta-1)$ for interpreting the coefficient?


title: Suggestion for doing research in data mining and machine learning

text: I am a fresh graduate student. I want to do research in machine learning and data mining. There is no professor in our department doing this! I want to try to do this by myself, at least for a while. But I don't where should I start. What books or review papers should I read at the beginning?


title: Can lmer() use splines as random effects?

text: Say we're working on a random effects model of some count data over time, and we want to control for some trends. Normally, you'd do something like: `lmer(counts ~ dependent_variable + (1+t+I(t^2)|ID), family="poisson")` to include a quadratic shape for t. Is it possible to use some more sophisticated smoothing techniques like a LOESS smoother or splines to model that relationship?


title: What Statistical method is necessary to determine the correlation amongst a set of individual correlations?

text: I have a question regarding the statistical method necessary to determine a correlation between a set of individual correlations. Below is an example of a situational question I currently face. To clarify as requested. I have gathered data on over 400 hotel properties. I then determined if there is a correlation between two variables for each individual property. 400 correlations are then returned independent of each other. Now lets say if those 400 properties belong to a chain. I would like to see if there is an overall correlation across all of the previous correlations tested for the 400 properties. Is there a statistical method that describes how to approach this problem?


title: Continuous time markov chains simple explanation

text: I am looking into continuous time markov chains which I think are capable of answering questions like "what is the probability that after time $t$ (where $t$ is in discrete time steps in this case) a transition from current state $i$ will occur", or more precisely "what is the probability that a transition from current state $i$ to a state $j$ (where a connection exists with $i$) will occur after $t$"? I am looking for a simple step-by-step example tutorial if exists or a simple explanation here if possible. For example assume the following dummy observed sequences: $o1 = 11112223$, $o2 = 111112223$ and $o3 = 11113$ and compressed to avoid self-loops: $o1 = 123$, $o2 = 123$ and $o3 = 13$ A discrete markov chain would give the following transition matrix... $T = \begin{bmatrix} 0 & 2/3 & 1/3 \\ 0 & 0 & 1 \\ 0 & 0 & 1 \end{bmatrix} $ If I calculate the transition matrix for the uncompressed observations I gain nothing other than strong biases on self-loops... The way I would approach this would be something like: $T(t) = A(t) * T$ where $A(t) = \begin{bmatrix}a_i(t)\end{bmatrix}$, $a_i(t) = e^{-(mean\_duration_i - time\_elapsed\_since\_i) }$ Any help or directions would be most appreciated. Thanks.


title: Translog Returns to Scale Variation

text: I estimated the so called translog production function: $$ y_{it} = \beta_0 + \beta_1 k_{it} + \beta_2 l_{it} + \beta_3 k_{it}^2 + \beta_4 l_{it} \times k_{it} + \beta_5 l_{it}^2 + \epsilon_{it}, $$ where $y$ is the natural logarithm of gross output; $l$ is the log of labour; and $k$ is the log of capital. I also have a third input - intermediate inputs - but for this discussion I have omitted it for simplicity. I estimated this equation using clustered robust (to autocorrelation) standard errors (assumes non-constant variance and zero covariance) with fixed time and firm effects. I want to calculate the variance of the returns to scale. Returns to scale is defined to be $$ RTS = e_k + e_l $$ where $$ e_k = \beta_1 + 2 \beta_3 k_{it} + \beta_4 l_{it} $$ and $$ e_l = \beta_2 + \beta_4 k_{it} + 2 \beta_5 l_{it}. $$ That is, returns to scale is the sum of the output elasticities of capital and labour. The returns to scale for the translog function depends on the level of capital and labour. Assuming independence between capital, labour, and the estimated coefficients we have $$ var(e_k) = var(\beta_1) + 4var(\beta_3) var(k_{it}) + var(\beta_4) var(l_{it}). $$ The variance of the $\beta$'s is outputted in Stata, so this is easily obtainable. My question boils down to: what should I use for the variance of capital and labour? I could calculate using all my observations, which would return a constant level of variance. However, the level of capital a firm has probably depends on capital last year. Therefore I am not confident that this assumption makes sense. Could I use the error covariance matrix estimates in the original regression to somehow get variance on capital that is robust to autocorrelation? Or is there another way to get this? The independence assumption between capital, labour, and the estimated coefficients is a bit dodgy as well. What is a good assumption to place here? Perhaps that the coefficients are independent of capital and labour, but capital and labour are correlated. Your help is much appreciated. NB: This is copied from math.stackexchange. I have flagged the one on math.stackexchange to be deleted.


title: What is the minimum of two independent variables?

text: I am having trouble to understand the min {X,Y} where X and Y are independent random variables. From the online source, it says that min{X,Y} < x if and only if X>x, Y>y. First of all, what exactly do you mean by the term "minimum " of the two random variables? Does it refer to the domains that X and Y can take, and you the minimum is the smallest value in that domain or what?


title: Gretl - How to extract sub-series?

text: I have series with 62 observations. How do I create a series with just the first 26 observations so that I can run an estimation on it? Sorry, really really beginner with gretl. I tried the square bracket but it just gives a single element. (I know general purpose programming languages very well.)


title: Intuition behind cross-correlation function interpretation vs. correlation of lagged time series

text: Can someone please explain the difference behind **WHY** the cross correlation function `ccf()` chooses to keep the same denominator for all lags and chooses to ignore the reduction in observations? Here's an example of the two methods not matching:               x = c(1,2,3,4,5,6,7,8,9,10)     y = c(3,3,3,5,5,5,5,7,7,11)     round(cor(x,y),3)     [1] 0.896          # Think "Lag -1"       # x[-10] = 1,2,3,4,5,6,7,8,9     # y[-1] = 3,3,5,5,5,5,7,7,11     round(cor(x[-10],y[-1]),3)     [1] 0.894          # Think "Lag -2"      # x[-10:-9] = 1,2,3,4,5,6,7,8     # y[-1:-2] = 3,5,5,5,5,7,7,11     round(cor(x[-10:-9],y[-1:-2]),3)     [1] 0.878          print(ccf(x,y,lag.max=3))     Autocorrelations of series ‘X’, by lag              -3     -2     -1      0      1      2      3       0.197  0.466  0.699  0.896  0.436  0.221 -0.018       Notice how the Lag-0 cases matches the output of ccf(), but the negative "manual" lags do not. This is because (to my understanding) the cross correlation function will construct the "covariance" (numerator) by comparing the lagged items to the "full" 10-item mean(x) and mean(y); in addition, I believe the denominator will keep the "full" series as well. At the end of the day, I can prove why the above `Lag -1 of 0.894` does NOT match the `ccf() -1 of 0.699` but I'm struggling to understand **WHY** the `ccf()` functions chooses to do what it does? I'm guessing it has something to do with adjusting for some sort of bias...?


title: Why do we use the determinant of the covariance matrix when using the multivariate normal?

text: I am not well versed in statistics. I wanted to know why we use the determinant of the covariance matrix instead of having the covariance matrix itself when writing down the multivariate normal distribution. Why do we do this and what is the intuition behind it? I have noticed that for a basic Maximum Lkikelihood estimation with normal errors can be written as follows for the univariate case: $\frac{1}{\sqrt{2\pi \sigma^2}}\cdot \exp(-\frac{1}{2}(\epsilon/\sigma)^2)$ and as follows for the multivariate case: $\frac{1}{\sqrt{2\pi |\Sigma|}}\cdot \exp(-\frac{1}{2}(\frac{\epsilon}{|\Sigma|})^2)$ What I want to know is why do we use the determinant of the covariance matrix in the multivariate case. In the univariate case we have the variance sigma^2, but in the multivariate case we write the determinant of the variance.


title: Weighting on tables

text: I have a question about terminology. Suppose I have data categorized by three factors A, B, and C, with cell means $\bar{y}_{ijk}$ and cell frequencies $n_{ijk}$, where $i$, $j$, and $k$ index A, B, and C respectively. And suppose I want to summarize the results for factor A by computing some sort of weighted means $WM_i$, averaging over indices $j$ and $k$ with weights $w_{jk}$. Consider these four weighting schemes:   1. Use equal weights, $w_{jk} = 1$.   2. Use weights of $w_{jk} = n_{+jk}$ (where "+" denotes summing over that index).   3. Use weights $w_{jk} = n_{+j+}n_{++k} $ (outer product of the one-factor marginal frequencies).   4. Use weights $w_{ijk} = n_{ijk}$ (the only one where we use a different set of weights for each $i$). Scheme 1 yields the "unweighted" or "least-squares" means, and scheme 4 yields the ordinary means for A, ignoring B and C altogether. Scheme 3 yields weighted averages over $k$ of weighted averages over $j$ (or vice versa). **My question is what to call these schemes** ; _e.g.,_ as a character argument in an `R` function. Preliminarily, I am calling them "equal", "proportional", "outer", and "actual". The first one is pretty obvious. But maybe there is some existing standard terminology for some or all of the others that I am unaware of (or have forgotten). Any suggestions?


title: prop.test: correcting for large number of observations

text: This is not really a coding question but more of a statistical question. I'm doing a proportions test on multiple proportions for many subjects. For example, subject 1 will have multiple proportions (multiple "successes per total trials"), and subject 2 will have multiple proportions. And for each subject we're testing if all these proportions are the same. For each subject, there are multiple proportions where there is number of successes per total trials. The proportions could range from being 30 successes out of 60 to like 300 successes out of a 1000 (just to show the range of trials and successes for each subject). Furthermore, for each subject, there could be varying number of proportions. Subject 1 could have 50 proportions, whereas subject 2 could only have 2. The idea is that we're trying to test that for each subject that all of their proportions are the same, and then reject if they are different. However, I'm realizing that subjects that have many more proportions, will have more significant p-values than subjects that only have 2 proportions, when using the prop.test. I was wondering if there is a way to approach this problem in a different way. Any sort of correction I could do, or take into account the number of positions. Any suggestions would be helpfil.


title: AIC and BIC for Support Vector Machine inside e1071

text: After training a support using the e1071 package of R, how can I calculate an information criterion such as AIC or BIC?


title: Shewhart's method for working out "sigma" of a sample of lead times

text: I am not a statistics expert so please bear with me. It is about making predictions about how long something should take based on measurements of how long things took in the past. (Things like fixing software bugs) Up until now I have been using a normal standard deviation to be able to state things like "the average bug will be fixed within 40 days, and 98% of bugs will be fixed within 115 days" since 115 is 3 standard deviations away from the mean. Then we can do things like establish SLAs (service level agreements) and promise that bugs will be fixed within 115 days or we pay a fee. However I have been told that it is not appropriate to use standard deviation as lead times are not normally distributed, which makes sense, the normal distribution would enable us to say a small percentage of bugs will be fixed in -10 days for example. Lead times cannot be less than 0. I was told they might be Poisson distributed, which I am not sure about either, and that I could use Shewhart's method for calculating "the mean plus 2 sigma" and that would be a better value to use for lead times. I haven't been able to find any information about this Shewhart's method. I did find something that said the standard deviation of a Poisson distribution was simply the square root of the sample size, but this doesn't make sense either, as I can imagine two different teams for example that both have 300 samples, and an average of say 40 days, but one team has 95% of lead times between 30 and 50, and none over 60, whereas the other team has 95% of lead times between 30 and 100, and some over 110. Simply using the square root of 300 will give the same standard deviation for both which does not seem to make sense. Anyway I guess my questions are:   1. What is Shewhart's method for working out something like "the mean plus 2 sigma", and   2. If you don't know Shewhart's method, what would be a good formula to use to make statements like "x% of bugs are fixed within n days" based on measured lead time data. (Knowing that lead times cannot be less than 0, and are probably not normally distributed etc) Thanks a lot, I am really out of my depth here.


title: Should the mean precision and recall be weighted?

text: I have a bunch of experiments in which I am calculating precision and recall. I want to present a mean precision and recall for these experiments. Should these values be weighted by anything?


title: Mean when computing correlation between samples of unequal size

text: This question is in some way similar to this one, but about another nuance. I have two time series (update: stationary - with both mean and variance equal over time) with missing values in one of them, such as   * 1,2,3   * 1,Absent,3 (in reality I have many more observations). I want to compute their correlation, so I will not use the 2nd time-point in these data. But should I use it for computing the mean of the first series? I think that including this point in computing the mean will lead to a more precise estimate of the population mean for the first variable, but don't know whether it is correct though.


title: Fitting a special mixed model in R - alternatives to optim()

text: I would like to do something in R that SAS can do using SAS's proc mixed (there is some way to do in STATA es well), namely fitting the so called Bivariate model from Reitsma et al (2005). This model is a special mixed model where the variance depends on the study (see below). Googling and talking to some people familiar with the model did not yield a straightforward approach that is fast at the same time (i.e. a nice high level model fitting function). I am nevertheless sure, there is something fast in R that one can built on. In a nutshell one is faced with the following situation: Given pairs of proportions $(p_1,p_2)$ in $[0,1]^2$ one would like to fit a bivariate normal to the logit-transformed pairs. Since the proportions come from a 2x2 table (i.e. binomial data) each logit transformed observed proportion has a variance estimate that is to be included in the fitting process, say $(s_1, s_2)$. So one would like to fit a bivariate normal to the pairs, where the covariance matrix $\Sigma$ _depends_ on the observation, i.e. $(\text{logit}(p_1),\text{logit}(p_2)) \sim N((mu_1, mu_2), \Sigma + S)$ where S is the diagonal matrix with $(s_1, s_2)$ and depends entirely on the data but varies from observation to observation. mu and Sigma are the same for all observation though. Right now I am using a call to `optim()` (using BFGS) to estimate the five parameters ($\mu_1$, $\mu_2$, and three parameters for $\Sigma$). Nevertheless this is painfully slow, and especially unsuitable for simulation. Also one of my aims is to introduce regression coefficients for mu later, increasing the number of parameters. I tried speeding up fitting by supplying starting values and I also thought about computing gradients for the five parameters. Since the likelihood becomes quite complex due to the addition of $S$, I felt the risk of introducting errors this way was too big and did not attempt it yet, nor did I see a way to check my calculations. Is the calculation of the gradients typically worthwhile? How do you check them? I am aware of other optimizer besides `optim()`, i.e. `nlm()` and I also know about the CRAN Task view: Optimization. Which ones a are worth a try? What kind of tricks are there to speed up `optim()` besides reducing accuracy? I would be very grateful for any hints.


title: Comparing relative frequencies between two groups

text: If i have the relative proportion for two sample. For example : group 1:   * Male: 50%   * Female: 25%    * children: 25%  group 2:   * Male: 50%   * Female: 20%    * children: 30%  How can I measure the % of similarity between the two group?


title: Covariance of timeseries

text: If I have data on a time series with the following structure                Day   Value, p       1     1.2       2     2.1        3     4.0       4     3.2      How do I calculate the first order serial covariance of the value $p$ on day $t$ with the value from the previous day, i.e. $\text{Cov}(p_t,p_{t-1})$? (I am using Excel)


title: Get probability of no positive results in a population?

text: I have a population of say 100,000 records that can be either Positive or Negative (true probability of Positive is very small, but unknown - may even be zero). I want to take a sample and estimate how many Positive cases are in the population, to be able to say something along the lines of "I'm 95% confident there are less than 50 positive cases in the population" or even better show my confidence there are no positive cases eg. "I'm 87% confident there are no positive cases". How can I do the following:   * What tests do I want to be doing here? From doing some research it looks like I want to use a Proportion Hypothesis test such as on this page. Is this correct? How would I set it up for my case?   * I can adjust my sample size, how do I determine the smallest sample size I would need for such a test.


title: Re-implementing R pt function

text: I would like to re-implement (or even better find an analogue in GNU GSL which I can easily use) R's "pt" function, since it's not native to my application's language. My use case is to convert a t value, to a p value, along the lines of:               if ( t >= 0 ) {       2 * pt(t, df, lower, tail=FALSE)     }     else {       2 * pt(t,df)     }      Is there a GNU GSL function that produces similar output? If not, what is the underlying algorithm used by pt? The R documentation suggests:               For the central case of pt, a normal approximation in the tails, otherwise via pbeta.          For the non-central case of pt based on a C translation of          Lenth, R. V. (1989). Algorithm AS 243 — Cumulative distribution function of the non-central t distribution, Applied Statistics 38, 185–189.      Googling "normal approximation in the tails" is inconclusive and I don't have access to the applied stats text.


What are some applications of Chinese restaurant processes? I'm trying to learn a bit about non-parametric Bayesian methods, starting with Dirichlet processes and CRPs, but all the tutorials I've found are about theory, without describing any applications in depth. Names of papers would be great. I'm not really looking for state-of-the-art applications, but just some "canonical" examples (say, in natural language processing) of why Dirichlet processes and Chinese restaurant processes are useful and why I should care.


title: Time Period Predictive Modeling

text: I have been implementing some classification algorithms (Naive Bayes, SVM etc) recently on the iris data sets to get head start into the data science field. I enjoy working on machine learning problems and exploring more. I am currently stuck with this problem of predicting variables based on just the time data. I have never worked with a time/date ranges data before. All I have in my data set is 3 columns (start usage date/time, end usage date/time & ID's of cars). I am looking to build a predictive algorithm/model on this past vehicle usage to predict period of times where the vehicles are likely to not be used. Also I am addressing the issue of improving the efficiency here and so I am looking to see which cars are safe to turn off at a specific date, and for how long. This is because of the battery usage on the cars. My columns look like this:                 ****Start Usage       *******      End Usage   **********      ID         11/5/2012 14:33:24.0          11/5/2012 14:33:28.6        2515476       11/5/2012 17:45:22.0          11/5/2012 17:46:28.5        2515476         11/5/2012 19:13:46.3          11/6/2012 9:39:19.1         2515490      What I have done so far is I have grouped by the car ID's to get the different times at the which they have been used and so each car has like 30(approx) different usage times. A time series algorithm is what pops in my head when I think of this scenario. I am not able to think far enough due to my lack of experience dealing with time/date ranges. I would appreciate if any one could point me in the right direction to this predictive analytics problem in terms of algorithm usage. Is time series the only way to do it? Any suggestions on what would be the best approach?


title: Correlation with related (but not repeated, per se) measures

text: I have 20 subjects split between two conditions. Each participant rated two items in terms of two attributes, A and B. How do I calculate the correlation between A and B. How do I calculate the correlation (I'm performing analysis in R) between A and B given that there are two ratings per subject? I'm getting confused about this because the two items are not the same, so this doesn't seem like repeated measures, but the two item's ratings on each attribute are likely to be correlated.


title: Keyword probabilities

text: I'm not a huge stats buff and am wondering what the best approach to my problem is. Say I have a list of PPC* keywords and if a desired action was taken. Let's say I put this data into two word count list: keywords with purchase and all keywords. For example: Keywords with purchase word count:               keywords     with purchase     all keywords          buy          30                35     best         35                35     purchase     15                16     review       5                 30      Intuitively, I'd say I can divide keywords within the desired action bucket against all keywords to come up with a probability that the keyword will spur an action. For example:               buy  = 30/35 = .85     review = 5/30 = .16      Here I can see that the word ‘buy’ is more likely to trigger a purchase. So, I guess I'm looking for feedback on this approach. Is it the best way? I've been reading about Bayes theorem and haven't completely wrapped my head around. Should it be applied in this situation or should I be employing a different approach altogether? * * * *Pay Per Click (PPC). So, if you type in "Buy Shoes" to Google, you'll see organic and paid listings. In general, PPC refers to paid listing which you can get data. So, you can see if somebody comes to your website from a paid keyword and purchases/or doesn't purchase something.


title: Data input uncertainty + Monte Carlo simulation + forecasting

text: Consider a variable $Y$ (e.g., temperature). Suppose that we were able to estimate this variable each year for the past $N$ years using some type of model. This means we have access to annual estimated values for Y (denoted as $Y_1, \ldots, Y_N$) and associated standard errors $S_1,\ldots,S_N$. The goal is to produce a point forecast for the value of $Y$ at the time N+1 which incorporates the uncertainty present in the estimated annual values $Y_1,\ldots,Y_N$. One way to proceed is to use Monte Carlo simulation to create an ensemble of $B$=100,000 (or a large enough number) of synthetic time series obtained by shifting each of the original values $Y_1,\ldots,Y_N$ by a random z score (i.e., Gaussian white noise), scaled by the standard error $S_t$. (This approach assumes that $Y_1,\ldots,Y_N$ are independent.) Each synthetic series can then be used to produce (I) a point forecast of $Y$ at time $N+1$ and (II) an interval forecast of $Y$ at time $N+1$. My question is: **How do we summarize the information conveyed by the simulated point forecasts and interval forecasts to quantify the uncertainty present in $Y_1,\ldots,Y_N$ and its impact on the forecasting output?** For point forecasts, it makes sense to report the simulation distribution of the point forecasts. But what aspect of this distribution captures the uncertainty (e.g., spread)? For interval forecasts, it is not clear (at least not to me) how to proceed. Is there any way to quantify uncertainty in the forecasting input (i.e., $Y_1,\ldots,Y_N$) when it comes to width and/or coverage of these intervals? (Maybe by using some type of retrospective performance of the prediction procedure?)


title: Can the Burnham-Anderson book on multimodel inference be recommended?

text: As motivated by the recent change of the default model selection statistic in the R's forecast package from AIC to AICc, I am curious whether the latter is indeed applicable wherever the former is. I have a series of questions with this respect and here is the first one. I know that to replace AIC with AICc everywhere is what the well-known book in (1) by Burnham and Anderson (non-statisticians), as summarized here, recommends. The book is sometimes uncritically referred to by younger statisticians, see e.g. comments to this blog post by Rob Hyndman, but the statistician Brian Ripley advised in radically different way:               “Burnham and Anderson (2002) is a book I would recommend people NOT read until      they have read the primary literature. I see no evidence that the authors have      actually read Akaike’s papers." [quoted from [AIC MYTHS AND MISUNDERSTANDINGS][4] by     Burnham-Anderson]      It does follow from what Ripley writes on the AIC and related theory that the warning should be taken seriously. I have both a good collection of Akaike's own papers and the Burnham-Anderson book. I will eventually have my own opinion on the quality of the book, but it will also help to know what the community of statisticians, both young and old, think on that. In particular, are there professors of statistics (or other good students of statistics) who explicitly recommended the book as a useful summary of knowledge on using AIC for model selection? Reference: (1) Burnham, K. P. & Anderson, D. R. Model selection and multimodel inference: a practical information-theoretic approach Springer, 2002


title: Calibrating multiple binary SVM classifiers for one-vs-all multi-class classification

text: I'm classifying text using the one-vs-all approach. There are three classes. I've trained 3 different binary SVM classifiers using 10-fold cross- validation. The accuracy of the binary classifiers is fairly high (above 0.8 for all three classes). Where I'm stuck is utilizing the result of the 3 classifiers to output the most likely class for a sample. Since it's not meaningful to compare the classifier outputs directly, I'm using Platt's technique (PDF is available) to calibrate them. The graph below shows the calibration results: ![enter image description here](http://i.stack.imgur.com/Mj6mJ.png) The top row is the histogram of the classifier outputs. The bottom row is the posterior probability for each example. This graph is analogous to Figures 1 and 2 of Platt's paper. Since the sigmoid fits look good, I think I've implemented the fitting algorithm *described by the pseudo-code at the end of the paper) correctly. However, I don't get good classification results in the end. Specifically, there are too many "class 3" false positive. There are almost none "class 2" true positives. My questions are:   * Can anybody see what could be wrong?    * What should I try next? Should I scrap one-vs-all and do one-vs-one?   * Does one-vs-one require a similar type of calibration? Or can you just use the sign of the SVM output to determine the winner of a particular vote?


title: Chi-squared test for detecting trending terms

text: I'm trying to find bursty(trending) terms in a text stream. There are two frame in a stream; expected frame and observed frame. For each frame, I tokenize documents(tweets/blog posts etc.) up into terms(words), so I have two term-frequency lists: For Expected Frame:   * term a, 2500   * term b, 2000   * term c, 1500   * ...   * term m, 23   * term x, 11 For Observed Frame:   * term a, 2600   * term b, 1900   * term m, 1500   * ...   * term x, 15 Now, I need a method to find 'term m' as a bursty term. I think chi-squared test can tell which terms is bursty, in other words which term has changed significantly. As you can see, I'm not interested in detecting whether whole observed frame is changed or not, but each individual term. My questions is that does chi-squared test fit this problem and what other methods are there? If I use chi-squared, what is the critical value for 'term m', how can I calculate? `"chi-value= power( 1500-23 ,2) / 23"` is ok?


title: Topic Words Selection in Topic Modeling

text: I understand how generative model of topic modeling works; for each topic there is a distribution of words, and for each document there is a distribution of topics. Question is how words are determined for each individual topic ? Does this change algorithm to algorithm ? **Edit:** Such as: Topic-1 : internet, web, search etc. Topic-2 : literature, novel, etc I wonder how "internet, web, search" are chosen for topic-1. Why "internet" and "web" take place in same topic instead of separate topics? Is it because they frequently coexist in different documents?


title: How to estimate peak load when volume increases

text: My company has a software system, which, amongst other things runs payroll for its client companies. In the near future, the number of companies running payroll will jump from 100 to 3000. I'm trying to investigate whether the system will be able to cope with the new load. I can do a bit of basic analysis on the logs and see that typically on the busiest days of the month, there are about 5x as many payroll runs as on the average day. The peaks could be because there is a group of companies which are connected and all run payroll on the same day, or it could be because if the numbers come out wrong, they have to make corrections and run payroll again until it's right. When the new companies come online, obviously the average load will increase by 2900%, as almost all companies only run payroll once a month. What I'm trying to estimate is how big the peaks are likely to be. Clearly as volume increases I would expect to see a smoothing effect, and the load charts to look less 'peaky'. Is there a good method for estimating the peak size, or a confidence value for the peak size, or can I only guess? Even just a point in the right direction would be great.


title: Regression with missing values

text: Please help. I have 3 IVS and 1 DV in my model. The 3 IVs are rated by 150 Middle Managers while 50 Directors rate the DV. I need to conduct multiple regression. However for cases (rows in SPSS), where Middle Managers rate the 3 IVs, I have missing values for Directors' rating of the DV. How should I conduct regression? Do I need to restructure the data? Thanks.


title: Alternatives to Horvitz–Thompson estimator

text: I am trying to estimate a mean of a population in a sample = ${\{x_1,..,x_n\}}$. I currently use 2 estimators:   1. Sample mean = $\frac{\sum {x_i}}{n} $   2. Horvitz–Thompson estimator (HT) = $\frac{\sum {p_i*x_i}}{\sum {p_i}} $ Where $p_i$ is the probability that the value$(x_i)$ is part of the sample during the sampling Note that it is a sampling with replacement, so there could be more than one $x$ that have the same value, i.e. value($x_1$) = value($x_5$) = 7 The HT estimator performed better than sample mean in my experiments. But for some samples, $p_i$ is unknown. So my questions are   1. Are there alternatives to HT estimator? (such as other estimators that use weighting scheme based on sample only)   2. Let's say I just know that the sample is from a sampling with replacement but don't know exactly how a sample is taken (e.g. someone gives me the sample and ask me to estimate the mean of the population.) and I want to use HT estimator. Can I approximate $p_i$ based on a frequency that the value($x_i$) appears in the sample? For example $p_i$ = $\frac{\#\text{value($x_i$) in a sample}}{n} $ Thank you!


title: Two-Component distribution modelling

text: I am trying to plot a two-component mixture distribution $F_X(x):.7N(0,1)+.3N(3,1)$ with density $$f_X(x) = {0.7 \over {\sqrt{2\pi}}} e^{−x^2/2} + {0.3 \over {\sqrt{2\pi}}}e^{−(x−3)^2/2} $$. I have done the histogram for it but having trouble with plotting the true density,so far this is my code               x <- rep(na,2000)     mu1 <- 0     sd1 <- 1     mu2 <- 3     sd2 <- 1          for(i in 1:n){         u <- runif(1,0,1)         if (u < .7){             x[i] <- rnorm( 1, mu1, sd1)         }else             x[i] <- rnorm( 1, mu2, sd2)         }          hist(x,prob=T)


title: Intuition for the standard error of the difference of sample means

text: I read in Wilcox, 2003 p. 247 that the standard error of the difference between two sample means is (assuming the normality and homoskedasticity assumptions): $\sqrt{\frac{\sigma_1^2}{n_1} + \frac{\sigma_2^2}{n_2}}$ Rather than simply adding the two sample standard errors as in: $\frac{\sigma_1}{\sqrt{n_1}} + \frac{\sigma_2}{\sqrt{n_2}}$ What is the intuition behind taking the square of the sum of the two variances divided into their respective sample size, rather than the sum of the standard errors?


title: How to compute the standard errors of a logistic regression's coefficients

text: I am using Python's scikit-learn to train and test a logistic regression. scikit-learn returns the regression's coefficients of the independent variables, but it does not provide the coefficients' standard errors. I need these standard errors to compute a Wald statistic for each coefficient and, in turn, compare these coefficients to each other. I have found one description of how to compute standard errors for the coefficients of a logistic regression (here), but it is somewhat difficult to follow. If you happen to know of a simple, succint explanation of how to compute these standard errors and/or can provide me with one, I'd really appreciate it! I don't mean specific code (though please feel free to post any code that might be helpful), but rather an algorithmic explanation of the steps involved.


title: Are the slopes of two segmented lines important in recognition of a breakpoint?

text: I have found one breakpoint in my linear regression but the slope of two lines have the same decreasing slope, I was wondering if the change of the slope is necessary in recognition of breakpoints.


title: R2, P-value, or accuracy of model?

text: Which is more important when selecting a model? The R2 or the accuracy of the model? Here are 2 different models using the same data with different transformations for the Xs, the Ys are a percentage and the same in both. This is a classification problem in which correct classification has a higher growth rate associated with it. Therefore, what matters in reality is accuracy and growth. So for what matters in reality I'd choose B, but how do its not just better due to randomness? **Model A**               Residual standard error: 0.03033 on 338 degrees of freedom     Multiple R-squared: 0.03352,    Adjusted R-squared: 0.02205      F-statistic: 2.922 on 4 and 337 DF,  p-value: 0.02127      AIC: -1410.358          Accuracy:     PointEst  Lower     Upper     0.6012085 0.5305505 0.6678887          Growth: 6.116476      **Model B**               Residual standard error: 0.03035 on 338 degrees of freedom     Multiple R-squared: 0.03032,    Adjusted R-squared: 0.01884      F-statistic: 2.642 on 4 and 338 DF,  p-value: 0.03368      AIC: -1413.970          Accuracy:     PointEst  Lower     Upper     0.6174699 0.5471007 0.6832359          Growth: 8.032874      As you can see, both models have the same number variables. The only difference is how I transform/treat the input variable data before running the regression. From the time where I noticed the difference between the two models, I ran them with an additional 20 observations and both performed equally as well in classification.


title: Why is Standard Deviation used to evaluate the effect of a change in an independent variable on the dependent one?

text: I'm reading some economics papers about the relationship between inequality and growth and some of them have sentences like these: > an _increase of 0.07 (one standard deviation_ in the sample) in the income > share of the top 20 percent lowers the average annual growth rate just below > half a percentage point and > the estimated coefficients imply that an increase in, say, the land Gini > coefficient _by one standard deviation_ (an increase of 0.16 in the Gini > index) would lead to a reduction in growth of 0.8 percentage points per year Why is "one standard deviation" used? Why is it preferred to a unitary change? Thanks


title: Test relationship of firm investment to GDP individually or in a panel regression

text: I need to test how responsive are firms' investment to the level of the country's GDP. I have 90 firms, all in the same country; 10 annual observations for the GDP; and 10 annual observations for each firm's investment level. Do I have to test the relation for each firm separately; if yes is it possible to run a regression with only 10 observations? Alternatively can I use panel data regression, noting that the level of GDP is the same for all firms in a given year; it changes only over time not across firms.


title: Multinomial Logisitic Regression in Mlogit For Dummies - Building a Model and Extracting its Coefficients

text: I am trying to do a multinomial logistic regression on some data that I generated. I am using R and the package mlogit. My data looks like the following:               Class X1 X2  X3     V +0.0655197 +0.6418541 +1.8110291     V-0.6713268 -0.0262458 -0.3602958     V +0.2357610 -0.3602958 -0.6943458     M +0.3900129 +0.5583416 -1.7800082     M +0.5714871 -0.2767833 +0.5583416     M +1.0732807 -1.7800082 -0.3602958     S +0.9553640 -0.3602958 +0.6418541     S +0.1139899 +0.1356030 +0.3889280     S +0.4717283 -0.2852090 -1.1229880      My model is               Class = B1*X1 + B2*X2 + B3*X3      So far, I have:               library(mlogit);     allData <- read.table("Features/AllFeatures.dat", header=TRUE);     allData$Class<-as.factor(allData$Class);     mlData<-mlogit.data(allData, choice="Class");     myData<- mlogit(Class~1|X1 + X2 + x3, data = mlData);     print(summary(myData));      Which gives me:               Coefficients :                                Estimate Std. Error t-value  Pr(>|t|)         S:(intercept)                -0.6832392  0.0951834 -7.1781 7.068e-13 ***     V:(intercept)                -0.6696254  0.0943282 -7.0989 1.258e-12 ***     S:X1                         -0.1362492  0.1134039 -1.2015    0.2296         V:X1                         -0.0052649  0.1128722 -0.0466    0.9628         S:X2                         -0.0198451  0.0973608 -0.2038    0.8385         V:X2                          0.0183261  0.0974789  0.1880    0.8509         S:X3                          0.1728694  0.1110473  1.5567    0.1195         V:X3                          0.0230260  0.1101147  0.2091    0.8344      However in the following: http://cran.r-project.org/web/packages/mlogit/vignettes/Exercises.pdf The author gets:               Coefficients :     Estimate Std. Error t-value Pr(>|t|)     ic -0.00623187 0.00035277 -17.665 < 2.2e-16 ***     oc -0.00458008 0.00032216 -14.217 < 2.2e-16 ***      How can I change my function call to get the same? I want the actual coefficients for the model, not for the comparisons. Also: What is the difference between 'wide' and 'long'? What form is my data in? Also: Do you have some mathematical references for this type of multinomial logistic regression aimed at engineers? Thanks!


title: A practical example for MCMC

text: I was going through some lectures related to MCMC. However, I don't find a good example of how it is used. Can anyone give me a concrete example. All I can see that is they run a Markov chain and say that its stationary distribution is the desired distribution. I want a good example where the desired distribution is hard to sample from. So we create a Markov chain. I want to know how to select the transition matrix so that its stationary distribution, of the Markov chain is the target distribution Thanks


title: How to check the ARIMA model by simulating innovations?

text: I am working on a monthly river flow time series and somebody has asked me to check my model by simulating innovations (or something like that). Can somebody help me about how can we simulate and how can check my model with that? I appreciate if you suggest a article or document in this regard. Also, if my question is wrong and this can not be done let me know so I will delete it. Thanks


title: Fitting multiple distributions with a joint parameter

text: I would like to fit multiple distributions that share one of their parameters. As a simple example, let's say I have two different datasets which I know follow a Gaussian distribution. let's say I know $\mu_1$ and $\mu_2$, and now I want to estimate $\sigma_1$ and $\sigma_2$. However, I know that $\sigma_1=\sigma_2$. Now, I would like to find this joint parameter $\sigma$ that best fits both datasets. The maximum likelihood approach tells us how to find parameters by maximizing the likelihood of one distribution, but how do I maximize the likelihood of both? Perhaps I should maximize the sum of the likelihoods or the product of the likelihoods?


title: Conditional Expectations and Variances

text: > $\textbf{Background:}$ When $\mathbb EX$ is hard to calculate, it is a > common trick to use the following formula: $$\mathbb EX=\mathbb E[\mathbb > E(X|Y)].$$ And similarly, $\mathbb VX$ can be calculated using the following > formula: $$\mathbb VX=\mathbb E[\mathbb V(X|Y)]+\mathbb V[\mathbb E(X|Y)].$$ > > $\textbf{Question:}$ Recently, I saw similar formula for conditional > expectation and variance as follows. $$\mathbb E(Y|X)=\mathbb E[\mathbb > E(Y|X, Z)|X],$$ and $$\mathbb V(Y|X)=\mathbb E[\mathbb V(Y|X,Z)|X]+\mathbb > V[\mathbb E(Y|X, Z)|X].$$ I do not know how to prove them. Could anyone > provide some hint or reference, please? Thank you!


title: Compute areas of Venn diagram given covariance matrix?

text: Given a covariance matrix, it's simple to draw a venn diagram representing the independent and shared variance of two variables. How might one go about the same (at least, computing all pertinent areas; actual visualization not necessary) for the 3-variable case? Is this even possible?


title: Inverted SPSS results: Logistic regression command vs. Genlin?

text: I want to do a logistic regression in SPSS. However, since I analyse unemployment spells the subjects are sometimes repeated (violating the independence assumption of the regression). One way of removing the within subject variation is by applying a Genlin model with the repeated subject subcommand (in essence a GEE model). Thus, I tried out a Genlin model with binomal probability and the logit link, comparing it to a standard logistic regression. I used the exact same variables in the two procedures. However, the results that was delivered from the Genlin procedure was inverted relative to that of the logistic regression. For instance: Exp(B) for women (of the independent variable sex/gender) was just above 2.0 in logistic regression while being at 0.49 in Genlin. The same happened with every independent variable.   * Any suggestions to why the results of the Genlin procedure is inverted?    * Is there any way to get the Genlin results in accordance to the logistic regression?


title: Beginner's question—Dealing with ordinal data (scoring of a medical effect)

text: I'm a biology PhD student dealing with ordinal data (scoring of a medical effect: 0–3 points) for the first time and I'm looking for information regarding differences in treating ordinal numbers in   * calculating Mean, Median, SD, SEM etc.   * performing significance tests between two ordinal datasets   * calculating correlations For the last point, I stumbled across the Spearman's coefficient. So far so good, but can I treat ordinal datasets the same way as cardinal datasets in terms of Mean, Median, SD, SEM and significance tests? Any comments appreciated, I'm also happy with good literature, if it's not too hard to understand.


title: Sampling extremely large and diverse dataset?

text: General theories would be great. My specific problem is I'm trying to find a specific portfolio of all the stocks in the market. The possibilities are huge because I need the stock combination(itself massive) and weights for each stock(even larger than the quadrillion combination). To make matters more complicated, I'm taking daily data but am trying to build a long term model. Here's an example of what I'm trying to do(its impossible because of the size but just so it makes it more clear), I take all the possible portfolio that have a return of more than 3% today(lets pretend I need to study this on the portfolio as a whole and it cannot be derived from individual stocks). Then, I look at that portfolio and see if it did the same thing yesterday, day before that and so on..after going back X days, then I return the portfolio. This way is not feasible because the entire population is huge(it would take a few billion years to process all combinations) and the results that work everyday are huge as well. If I'm trying to build portfolios how can I intelligently draw samples from this? I've tried genetic algorithms to mutate successful ones but the problem is it may work but I might be in a local peak, so I tried to do random sampling as well but that didn't work very well. I'm trying to think of ways of not sampling randomly. Any ideas? (if its not clear here's another problem that is more general but also illustrates where I'm lost when it comes to sampling): say you wanted to know what percent of the entire world liked Democrats every year. The world is too large to poll everyone and its very diverse so doing it randomly may give you areas with high population(i.e. india/China but you'd miss out on smaller countries). Plus to complicate things, the answer this year will cause you to go back to look at the answer they gave last year and so on. I'm interested in the sampling techniques, do you somehow segment the world into smaller groups and equally take samples(then a small country would have just as many votes as a large one? or do you do it randomly? etc..). Thanks. p.s. I know this is highly dependent on my research set but I'm trying to learn about general ideas so I can use existing methods as a model of thinking about this problem. I'm have taken basic stats courses back in school but I'm not sure if this question is about sampling or if its about random sampling. Any suggestions are welcomed.


title: Can the ratio importance sampling estimate by made to be unbiased with resampling?

text: Consider approximating the following integral: $$ \mathcal{Z} = \int h(x) \pi(x) dx $$ Where $\pi$ is known only up to a normalizing constant, that is, $\pi(x) = \hat{\pi}(x)/\mathcal{Z}_\pi$. We can simulate $\{x^{(i)}\}_{i=1}^m$ from an appropriate proposal distribution $q$ and use the ratio importance sampling estimator: $$ \mathcal{Z} \approx \sum_{i=1}^m w(x^{(i)})h(x^{(i)}) / \sum_{i=1}^m w(x^{(i)}) $$ Where $w(x^{(i)}) = \hat{\pi}(x^{(i)})/q(x^{(i)})$. This is a consistant estimator and actually even has lower asymptotic variance than the standard version (which can be used when $\pi$ is normalized) but it cannot be said to be unbiased. In the event that one needs provably unbiased samples could this be achieved through the use of the resampling bootstrap? That is, we collect our weighted samples as usual: $\{w^{(i)},x^{(i)}\}_{i=1}^m$ then proceed to a resampling stage to generate a set of unweighted samples: $\{\hat{x}^{(i)}\}_{i=1}^m$ which are distributed according to $\pi$. At this point we can use the standard empirical average to yield an **unbiased** estimate: $$ \mathcal{Z} \approx m^{-1} \sum_{i=1}^m h(\hat{x}^{(i)}) $$ Of course, this is not likely to be a practical approach since the bias in the original estimate is not typically a problem, my question is my attempt to settle an argument over the fact that the bias doesn't matter, and even if it does, it can be worked around in this way.


title: $M$-estimation in multivariate linear regression model in R

text: Is there any function for $M$-estimation in multivariate linear regression model in `R`. I can estimate the $\beta$'s in my model by using the `rlm()` by rewriting the $y$-variables into one column but, I would like to use one function to get the $\beta$'s.


title: Variance of a sample of random variables

text: I have a sample of 100 items, each associated to a random variable for which I can compute expected value and variance: $X_1, X_2, ..., X_{100}$. From these, we can define the mean $\overline{X}=\frac{1}{100}\sum{X_i}$. I'd like to test the hypothesis $H_0:\mu=0$, (where $\mu$ is the true population mean from which the 100 items were sampled) but for that I need the variance of the sample of $X_i$'s. So the question is: how do I compute the variance of the sample, given the individual expectations and variances of the random variables? Thanks a lot! **EDIT: More info about the question** I have a set of 100 items, and there is a function that assigns a score to each of them. The problem is that computing that function is actually very expensive, so instead I have a process with which I can estimate the score: the more effort (e.g. money) I put into the process, the better the estimate (less variance). So initially the variance of each estimate is maximum, it decreases as I put more effort into the process, and eventually variance is 0 and expectation equals the actual score if I run the process completely (i.e. the original function). Those 100 items represent just a random sample of a wider population of items, so I'd like to test the hypothesis of the population score being different from (or larger than) zero.


title: Modeling trimmed mean

text: In OLS, the conditional mean $E(Y \mid X)$ is modeled as a function of some regressors $X$, i.e. $$ E(Y \mid X) = X \beta. $$ Is there a regression technique that allows to model the conditional trimmed mean of $Y$? (Least absolute deviations will lead to the extreme case of maximal trimming.) Is there a direct relation between this question and least trimmed squares estimator?


title: Calculate effect of independent variable

text: If $a\times b\times c=t$ and I change $a$, $b$, and $c$, how do I calculate the effect of each change. I believe the equation is: $(a+\Delta a)(b+\Delta b)(c+\Delta c)=(t+\Delta t)$, and then I solve for $\Delta a$, $\Delta b$, and $\Delta c$. Is this correct? And how do I describe what I'm solving for?