title: Use of multiple regression for data on 20 countries representing nearly the whole population

text: I have operational expenditure data for 20 countries (representing 90% of the whole population) which I want to compare with other variables from the countries (e.g. GDP, Population etc). Do I still need to worry about the P-value/t-stats or can I just focus on the coefficients and the R square value of my (linear) regression. How many different variables could I include given the number of observations? Thanks


title: WinBUGS sript: display correlation

text: Does anybody know if it is possible to display the correlation matrix using WinBUGS script? (with the GUI it can be found in th Inference menu -> correlations -> print) Thank you everybody!


title: Clustering coefficient for a clique

text: I would like to understand how to solve this exercise about clustering coefficient for a clique. As shown in the picture below if node pairs (a; b), (a; c), (a; d), (b; c), (b; d) are linked, then the pair (c; d) is likely to be linked in the future. To quantify the amount of "quad closure", we define a quad clustering coefficient" as   Q =Number of cliques of size 4 / ( Number of connected quadruples with 5 edges / k )   where k is some normalizing constant. But this definition is incomplete unless we specify the value of k to normalize Q. Which is the value of k such that the value of Q for a clique is exactly 1 ?   ![enter image description here](http://i.stack.imgur.com/W8x0N.png)


title: Influences of scored/weighted ratings

text: I am experimenting with a scheme for aggregating $N$ scored ratings. The gathered ratings themselves are integers, e.g. a chosen value between $1$ to $5$, and each rating $r_i$ is assigned with a score $s_i$. That is, a data point/record here is a $(rating, score)$ pair. However, $\sum_is_i\neq1$. Naturally, if we sort these $N$ ratings by their scores, we end up with a ranking of them. Now, I would like to aggregate all the ratings into a single measure, but face a potential problem: for example, given two records $(r_1=5, s_1=0.7)$ and $(r_2=4, s_2=0.8)$, if we simply multiply the rating with its associated score and use the resulting product for aggregation, then we will see a larger influence from $r_1$ than from $r_2$ on the result of aggregation (as $5\times0.7>4\times0.8$) even though $r_2$ scores _considerably higher_ than $r_1$. I wonder what are the appropriate/fair methods/techniques for aggregating in this context. Thanks! Edit: in response to whuber's comment below, here are some clarifications   * The score for a rating is separately calculated for each rating. We can think of the score as a measure of the 'soundness' of a given rating. As the scores may well be different for different ratings, we can then sort the ratings by their scores in descending order to produce a ranking of ratings.   * One score for one rating, so more precisely, the dataset consists of $N$ number of records each of which is a $(rating, score)$ pair, e.g. $(r_1, s_1)=(5, 0.7)$, $(r_2, s_2)=(4, 0.8)$ as the examples given above.   * By aggregation, I mean combining all given ratings (i.e. the whole dataset) into one single number which is still within the range of individual ratings. In other words, if any individual rating $r_i$ is $1\leq r_i \leq5$, then the aggregated $r^*$ should also be $1\leq r^* \leq5$


title: R randomForest - is my variable treated as quantitative or categorical?

text: I am using randomForest and obtain a model:               Shoppers.rf <- randomForest(repeater ~ . - id, data=trainData,                       importance=TRUE, ntree=1000, mtry=15, nodesize=50, maxnodes = 100)     varImpPlot(Shoppers.rf, type=2)      ![enter image description here](http://i.stack.imgur.com/EX63b.png) The diagram tells me that the 'brand_penetration' variable is the most interesting to predict the response. But if I look at the raw data, I get a feeling there the algorithm is treating 'brand_penetration' as a categorical variable instead of a quantitative one. Here is a barplot of the success ratio per brand_penetration: ![enter image description here](http://i.stack.imgur.com/RaFpa.png) How can brand_penetration possibly be useful to predict the response given the plot above, except if it is used as a categorical variable? But since I pass this variable as a 'int', I am puzzled. Should it not be used as a quantitative input by randomForest? I have tried converting to 'num', but get the same result.


title: Steady state probabilities for a continuous-time Markov chain

text: I have a finite state and time-homogeneous continuous-time Markov chain (CTMC) which is not irreducible. Will steady state probabilities exist for this CTMC? How to prove this?


title: How to model the relationship between geocoded and ungeocoded sales data?

text: I am trying to model sales data for stores at the Census block group level in order to predict sales at potential new restaurants. For example, I know that store 2, which has a giant flashing neon sign, has $2K in sales from block group 101, which is 2.5 miles away from the store 2 and where there are 600 households and 50 people living in college dorm. So far this is pretty standard. The fly in the ointment is that the average store has ~30% of the sales data that cannot be geocoded for some reason (new construction, college dormitories, military bases, lazy employees who take down the address in shorthand, and so on), so that I only know the store that handled the sales and not where those customers reside. My approach to modeling the ungeocoded data starts with aggregating demographic data from around the store's trade area and all the ungeocoded sales, so that even if I don't know where those customers are, I can at least take a stab at understanding the sales behavior based on what's around the store. For example, if my store is near a college campus or has lots of construction, I would it expect it to have higher ungeocoded sales, all else being equal. This works reasonably well, but the geocoded and ungeocoded sales models are not linked in any way, which is problematic. Essentially, my geocoded sales are measured with non-spherical error that is correlated with my explanatory variables. It's also the case that ungeocoded sales generally increase with geocoded sales. I tried to remedy the first by including the fraction of total sales that are ungeocoded in the geocoded sales model, and total geocoded sales in the ungeocoded model, but I don't know how to define those variables for potential sites whose sales I am interested in forecasting. I guess I can just set ungeocoded sales at 30%, then predict geocoded sales, and use that to forecast ungeocoded sales, but is there a better way to link the two models for the better estimation and forecasting?


title: Regression to remove consecutive effects

text: Suppose I want to know how my observable X depends on Y and on Z. But Y itself depends on Z (e.g. they are correlated). I want to know the effect of Y on X after removing the effect of Z on both X and Y. What is the right regression procedure for that? - Should I regress Z on X, take the residuals, and regress those on Y? Like,               X = a * Z + e1     e1 = b * Y + e2      I know I could do PCA, but I want the first order effect to be coming from Z exactly, the second order effect e2 to be orthogonal to e1.


title: Post-process the output of a Multinomial Naive Bayes text classifier

text: I have a multinomial text classification application where there are other features than the words in text which can be useful to do the classification e.g, contains email address, contains an URL, number of words, uppercase/lowercase. What I'm planning to do, is:   * Build a multinomial Naive Bayes classifier (with NLTK), which gives me the probability distribution over each category (percentage number)   * Extract other relevant features (has_email(0/1), has_url(0/1), ..)   * Use all that data (probability distribution plus other relevant features) as an input for a final step classifier. I have a few thousand documents for the training. The documents are short texts (from tens to few hundreds words). The multinomial classifier has 6 categories. I don't have much experience in NLP and machine learning. Is my approach the right one, or there's a better solution? Which would be the best type of classifier for the final step? Thanks in advance


title: Journals in statistical learning / machine learning

text: Can you please name some major and minor journals publishing articles in the field of statistical learning / machine learning. regards Anthony


title: How to measure the reliability of a consensus ranking (problem from Kemeny-Snell book)

text: Suppose that $k$ experts are each asked to rank a set of $n$ objects in order or preference. Let allow ties in the rankings. John Kemeny and Laurie Snell in their 1962 year book "Mathematical models in the Social Sciences" propose to solve next problem: > PROJECT $1$. **Develop a measure of the reliability of a consensus ranking > by $k$ experts.** For example, this may be based on the largest possible > change that can be brought about by changing the ranking of a single expert. > (Attention must be paid to the possibility of multiple consensus rankings.) > Prove some theorems concerning the most and least reliable consensuses > possible for a given $k$. The book gives notation for rankings and method for rankings aggregation (i.e. getting one "collective" ranking from many "individuals"). But no answer given for the problem above. First, I thought about Kendall's $W$ coefficient of concordance, but it looks like it doesn't suit. Any ideas are welcome!


title: Model checking and generalized multilinear regression

text: In my research I am analyzing the quality of game player softwares by making them play multiple times against a certain opponent and counting the number of victories. It seems to be very clear that the random variable of #victories should be modeled by a binomial distribution (as each game is either a victory or a defeat, and they play a fixed number of games). So, I have two different software types and I change a certain numeric factor for both of them. My hypothesis is that by increasing this factor the #victories will increase for both software types, however it will increase faster for one than for another. In order to test this hypothesis I am performing in R a generalized linear regression with many variables, assuming a binomial model. In R syntax, that would be:               glm(formula = cbind(Win, Lose) ~ Size + Group + Size:Group, family = binomial),       where `Size` is the numeric factor I was talking about and `Group` is a factor that defines the software type. I am looking at the coefficient of `Size:Group` to test if there is an interaction between these two factors (and therefore `Size` affects one `Group` type more than another). At first sight I get a good result, as the coefficient `Size:Group` is `-0.023752`, with `p = 0.000164`. However, R also prints the `Residual Deviance`, and in that case it is `34.027` on `10 degrees of freedom`. If I get the chi-squared value, I find: `1-pchisq(34.027,10) = 0.0001827632`. As far as I understand, this means that there is a statisticaly significant difference between each point and the prediction of the model. Therefore, seems that the binomial model is not adequate. Is that right? This is very puzzling to me, because by the nature of the `#victories` variable, it should be a binomial. If I just perform a simple multiple linear regression over the proportions for each `Size` (in R: `lm(formula = Win ~ Size + Group + Size:Group)`), I still find a negative coefficient for `Size:Group`, but a much higher p: `p = 0.079729`. Testing the model itself, however, gives a much better result (`Multiple R-squared: 0.8028, Adjusted R-squared: 0.7437, F-statistic: 13.57` on `3` and `10 DF`, `p-value: 0.0007372`). I am very confused, because the best model should be a binomial. Also, this linear regression is only considering the averages, there is no information about how many samples for each average, etc. Moreover, `p = 0.079` is a bit too high, and doesn't look like I could improve over that by running more samples, because the simple linear regression is actually not using any information about sample size. I also ran linear regressions for each `Group` type, testing only the effect of `Size` on `#victories`. I tried both a simple linear regression (`lm(formula = Win ~ Size)`) and also a generalized linear regression (`glm(formula = cbind(Win, Lose) ~ Size, family = binomial)`). The graphs of the simple linear regression actually "look better" than the ones for the generalized linear regression (I mean, the lines are closer to the data points). In other parts of my paper draft I perform point-to-point comparisons of the proportions (for a fixed Size, I compare the proportions of the two different `Group` types, and for a fixed group, the proportion of two different values for `Size`). For these comparisons I am assuming a binomial model (I am using `prop.test()` in R), due to the nature of the data. But maybe it is weird to assume a binomial model for these comparisons, but perform a simple linear regression instead of a generalized one? So, overall I am a bit lost about all this.


title: How to plot spectra of an AR(2) process

text: I am stuggling with this problem and was hoping to find some guidance to answer it. Let $y_t=\phi_1y_{t-1}+\phi_2y_{t-2}+\epsilon_t$, with $\epsilon_t\sim N(0,1)$. Now, I want to plot the spectra of $y_t$ in the following cases: Case 1: When the AR(2) characteristic polynomial has two real reciprocal roots given by $r_1=0.9$ and $r_2=-0.95.$ Case 2: When the AR(2) characteristic polynomial has a pair of complex reciprocal roots with modulus $r=0.95$ and frequency $2\pi/8$. Now, before plotting the spectra of $y_t$ in the following cases, I have tried to make use of the following important facts. The AR(2) process $y_t=\phi_1y_{t-1}+\phi_2y_{t-2}+\epsilon_t$ has the general linear process form $\psi(u)=1/(1-\phi_1u-\phi_2u^2)$ and hence $$f(\omega)=\frac{v}{2\pi}|(1-\phi_1e^{-i\omega}-\phi_2e^{-2i\omega})|^{-2}$$ This can be expanded to give $$f(\omega)=\frac{v}{2\pi[1+\phi^2_1+2\phi_2+\phi_2^2+2(\phi_1\phi_2-\phi_1)\cos(\omega)-4\phi_2\cos^2(\omega)]}$$ Now if the roots are real, then $f(\omega)$ has a mode at either zero or $\pi$; otherwise, the roots are complex conjugates and $f(\omega)$ is unimodal at $\omega=\arccos[-\phi_1(1-\phi_2)/4\phi_2]$ lying strictly between zero and $\pi$. So if anyone could help explain to me how I am supposed to relate the above facts with the two different case that would be very helpful. I guess what I am struggling with is what values to plug into $f(\omega)$.


title: How to calculate the average selling period from inventory list (uneven intervals)?

text: Will I be able to calculate the average sales time (number of days) for an item, simply from the following 6 data points? Or will I need to repeatedly collect data, such that I for each day (say in the next 14 days) have the following data?               PUT UP FOR SALE WITHIN    NUMBER OF ITEMS     The last day                       17,319     The last 2 days                    30,832     The last 3 days                    44,328     The last 7 days                    96,829     The last 14 days                  177,765     Total                           1,272,306      Once an item is sold it no longer appears in the list above. And all items are removed after approximately 2 months. The guess (of the top of my head) would be to estimate some sort of Sigmoid function, and then calculate the average sales time on that. But having only 6 observations -- and not knowing which initial parameters I should start with -- means that it doesn't converge to anything reasonable. I also remember learning something about the Possion distribution that related to time, but I've completely forgotten it.


title: Better to use cross-validation or training/holdout for predictive modeling?

text: I'm working with a small behavioral health care dataset (22,090 records) and have been asked to develop a predictive model that identifies patients at higher risk for hospitalization & health costs in FY2013 based on information in FY2012. The final predictive model will eventually be used to flag high risk members in FY2015 based on FY2014 data. In order to compare the performance of different methodologies (CART, SVM, logistic regression, etc.) and avoid overfitting, I'm considering two options:   * Use 5 or 10 fold cross validation on my existing data FY2012-FY2013.   * Train competing models on FY2011-FY2012 data and compare their performance on the FY2012-FY2013 dataset. Which approach will help me find the best-fitting predictive model: cross- validation or training/holdout?


title: How to add confidence intervals to predicted data when the response variable is log transformed

text: I'm trying to predict the impact of avg grasshopper density(ghavg) on plant biomass. Both the response and predictor are continuous data. I had to log transform biomass (logmass) for a normal distribution, giving the following dataset:               logmass = c(8.032393925,7.439531107,7.307924891,7.036315375,6.679316231,6.545784839,6.414481385,6.39297518,6.209312602,6.209312602,5.698486978,4.862609605,4.367692388,3.608137836)          ghavg = c(30.4,30.4,7.7,124.8,7.7,7.7,123.2,30.4,21.1,21.1,21.1,123.2,47.9,124.8)      In R I ran a glm:               biomass<-glm(logmass~ghavg,data.txt).       To predict the effect of grasshoppers on plant biomass (up to 300 m2) I ran:               nd<-data.frame(ghavg=0:300)     pred_mass<-predict.glm(biomass,type="response",newdata=nd)                    Then back transformed the output with:               trans_pred<-(exp(1)^(pred_mass))                              I think this works well but I would now like to add confidence intervals to the prediciton. I tried:               pred_massSE<-predict.glm(biomass,type="response",se.fit=TRUE,newdata=nd)       but the standard errors are clearly off when back transformed. Any assistance on how to properly add CI's to predicted values from transformed data would be greatly appreciated


title: Stationarity in OLS time series and asymptotic properties

text: I think I lack somewhat deeper understanding of this topic, but I thought stationarity is required in order for OLS to have asymptotic properties. "But stationarity is not at all critical for OLS to have its standard asymptotic properties"(Wooldridge, 2012) I thought stationarity is needed or otherwise OLS would not be consistent, but I guess I'm wrong. Could someone tell me why stationarity is not critical for LLN? Thanks in advance!


title: What does it mean that random forest is patented?

text: Looks like random forest has been patented multiple times: http://www.faqs.org/patents/app/20120321174 http://www.google.com/patents/US6009199 Any idea what this means in practice? Is any company who uses this algorithm liable for a lawsuit? If so half of silicon valley will be affected.


title: Fixed effects analysis in Stata

text: I have a large data set on which I need to run a fixed effect regression. The two variables that I want to control for: year (2001–2008) and industries (20 types). But I apparently need to set the data for a panel setting. How do I do this? With the menu bar `Statistics`$\rightarrow$ `Longitudinal/panel data`$\rightarrow$ `Setup and utilities`$\rightarrow$declare data set to be panel data, I choose the following: Panel ID variable = Industry   Check Time variable and choose Year But I get an error "repeated time values within panel" r(451). Am I using the right steps? Help would be highly appreciated :)


title: How to model the factors extracted from factor analysis with multinomial response

text: I want to model the 5 factors obtained from factor analysis with multinomial response. I calculated the factor scores by taking the mean of the raw scores of correlated variables in a factor. When I applied multinomial regression on these factors as explanatory variables along with the multinomial response, all these factors are insignificant. Where I am doing something wrong? Has anyone faced this type of situation?


title: How to estimate a pdf of x under the model of y = x+n, when the pdf of y and the pdf of n are given

text: I guess I come up with a classic question, but I failed to find any useful solutions by far. My question is about the following model $$y=x+n$$ where $x$ is a hidden random variable that cannot be observed, $n$ is the random variable of white noise, namely $n$ follows a known Gaussian distribution of variance $\sigma^2$ ($f_n= {\cal{N}}(0,\sigma^2)$), and $y$ is a random variable that we can observe. Suppose we know the distribution of $y$ as $f_y$, and I wonder how to find the pdf of $x$. Theoretically, it seems to be equivalent of the sum of two dependent random variables i.e. $x = y-n$, but in this question how $y$ and $n$ are correlated is unknown. Is there any existing solution? Thanks


title: Interpreting logistic regression

text: I need to perform a logistic regression to to see if a group of variables which are found to be significantly associated with an outcome (by univariate tests) have significant impact on the outcome when they are put together. In this case, from my reading, I gathered that in SPSS I need to use the enter method. When I applied this method, my new model decreased the prediction percent (not sure if this is the correct term) in the classification table in comparison to the constant only model (90.5 in constant only to 89.8 in new model). I have excluded multiple collinearity earlier.               Omnibus Tests of Model Coefficients                          Chi-squaredf    Sig.     Step 1  Step    27.441  7   .000              Block  27.441  7   .000              Model  27.441  7   .000          Model Summary                Step    -2 Log likelihood   Cox & Snell R Square    Nagelkerke R Square       1               164.185a                .086                   .184          (a Estimation terminated at iteration number 6 because parameter estimates changed      by less than .001)                    Hosmer and Lemeshow Test                 Step     Chi-square    df     Sig.     1            6.317      7     .503      Given my objective of seeing if the independent variables when put together can have a significant impact on the outcome or not, can I use this model for that purpose to identify the variables that have and do not have independent effects on the outcome when compared together? In answer to the comments:   1. As for the "measure of predictive accuracy", could I please know it in layman's terms please? My statistical vocabulary is inefficient to answer precisely though I understand what is should mean.   2. I ran collinearity diagnostics to see there is an issue of multicollinearity between independent variables. I cannot paste the ROC curve itself, but the tables is given below, given the values I guess I may still use the model as it is rather than changing any other parameter. Opnion on this is welcome. Thank you everyone for suggestions.               Area Under the Curve                     Test Result Variable(s): Predicted probability     Area    Std. Error a   Asymptotic Sig. b    Asymptotic 95% Confidence Interval                                                 Lower Bound             Upper Bound     .756    .042           .000              .67                           .839          The test result variable(s): Predicted probability has at least one tie between the      positive actual state group and the negative actual state group. Statistics may be biased.           a Under the nonparametric assumption                     b Null hypothesis: true area = 0.5                Now I have an additional problem, or two rather 1\. I have another outcome variable, again my onjective is to identify independant predictors for the given outcome amongst predictors with a p<0.01 in univariate analysis. I tried different combinations of covariates manually using the enter method, and non of a combinations seem to give a statistically acceptable models in terms of correctly classified proportions in classification tables or in ROC curves. Individially taken, these variables would give acceptable values (>0.8) in ROC despite not improving the correctly classified rates. Can I take this as being evidence of these covariants not haveing independant effects on the outcome? (number of independant variables available 7) As the dependant variable in this is 30 day mortality, I think I can try Cox regression. However, the only time dependant variable I have is length of hospitalization, I am not too sure if I can use this at the "time" cage in SPSS. Thank you in advance. (I am not too sure if I need to post this as a new question)


title: Trouble in comprehending question: Simple linear regression - LS and MLE

text: I've been given the simple linear regression model: $y_i = β_0 + β_1x_i + ε_i$ Under the assumptions of a simple linear regression model, the question they ask is: > Assuming the usual model assumptions hold, show that the least squares and > maximum likelihood estimators of $(β_0 β_1)'$ coincide What I can't seem to comprehend is the $(β_0 β_1)'$ part. I've never seen this before. Can anyone explain what $(β_0 β_1)'$ means? Also, how would this proof of LS and MLE differ from the usual $β_0$ and $β_1$ calculations? Note. I'm not asking for the proofs, I simply want to understand what they are trying to ask.


title: Evaluating correlation: which test to choose

text: I want to evaluate the correlation between two samples of data (300 samples, 2 parameters). I've done some researches and I reached these conclusions (correct me if I'm wrong):   * Pearson's correlation is used when: the variables are normally distributed, their relationship is linear and there are few outliers;   * Spearman's correlation is used when: the variables are NOT normally distributed, their relationship is not linear but monotonic (crescent/decrescent);   * Kendall's correlation is similarly to the Spearman's one, but computed differently. Now let's say I have 300 samples of 2 parameters (Throughput and Latency) and I want to evaluate their correlation. The two parameters are not normally distributed, their relationship is linear but not monotonic (there are some spikes). ![enter image description here](http://i.stack.imgur.com/6ds89.png) Given these assumptions I don't know which test to use because none of these characteristics satisfies ALL the "requisites" of the three tests. Graphically it's easy to say they are heavily correlated: ![enter image description here](http://i.stack.imgur.com/MZlXv.png) Anyway I did some analysis with SAS JMP and obtained the following results (for "rho"):   * Pearson: -0,70   * Spearman: -0,54   * Kendall: -0,40 The last two rejects the null hypothesis that X and Y are indipendent (there is the asterisk). What are the results of this analysis (if it's correct)? Another example consists on evaluating the correlation between two RAM's parameters: free and cache. When the system is not saturated they are heavily correlated. Again their distribution is not normal and their relationship is monotonic, so I assume I have to use Spearman/Kendall instead of Pearson. Unfortunately their relationship is linear, so basically I don't know what to do. (maybe I misunderstood some of the basic concepts like "linear relationship". Correct me if I'm wrong)


title: How to analyze longitudinal count data: accounting for temporal autocorrelation in GLMM?

text: Hello statistical gurus and R programming wizards, I am interested in modeling animal captures as a function of environmental conditions and day of the year. As part of another study, I have counts of captures on ~160 days over three years. On each of these days I have temperature, rainfall, windspeed, relative humidity, etc. Because the data was collected repeatedly from the same 5 plots, I use plot as a random effect. My understanding is that nlme can easily account for temporal autocorrelation in the residuals but doesn't handle non-Gaussian link functions like lme4 (which can't handle the autocorrelation?). Currently, I think it might work to use the nlme package in R on log(count). So my solution right now would be to run something like:               m1 <- lme(lcount ~ AirT + I(AirT^2) + RainAmt24 + I(RainAmt24^2) + RHpct + windspeed +            sin(2*pi/360*DOY) + cos(2*pi/360*DOY), random = ~1|plot, correlation =           corARMA(p = 1, q = 1, form = ~DOY|plot), data = Data)      where DOY = Day of the Year. There may be more interactions in the final model, but this is my general idea. I could also potentially try to model the variance structure further with something like               weights = v1Pow      I'm not sure if there is a better way to do with with a Poisson mixed model regression or anything? I just found mathematical discussion in Chapter 4 of "Regression Models for Time Series Analysis" by Kedem and Fokianos. It was a bit beyond me at the moment, especially in application (coding it in R). I also saw a MCMC solution in Zuur et al. Mixed Effects Models book (Chp 23) in the BUGS language (using winBUGS or JAG). Is that my best option? Is there an easy MCMC package in R that would handle this? I'm not really familiar with GAMM or GEE techniques but would be willing to explore these possibilities if people thought they'd provide better insight. **My main objective is to create a model to predict animal captures given environmental conditions. Secondarily, I would like to explain what the animals a responding to in terms of their activity.** Any thoughts on the best way to proceed (philosophically), how to code this in R, or in BUGS would be appreciated. I'm fairly new to R and BUGS (winBUGS) but am learning. This is also the first time I've ever tried to address temporal autocorrelation. Thanks, Dan


What is the best way to strip weekly and seasonal noise from a time series data set? Any recommendations on different approaches and there relative benefits?


title: Comparing the relationship between two sets of data where the one of the sets has more values than the other

text: I am working on a questionnaire with which I am supposed to find the relationship between one parameter and the other parameters, but the stationery parameter has six values and all the others have different sets of values. I was trying to work with the correlation formula in Excel but that is not working. For example, the mean of the stationery parameter is 5.227642276 and the values of the stationery parameter are:               5.829268293   5.390243902   5.146341463   4.853658537   4.926829268   5.219512195      and the mean of X(with which I want to find the correlation) = 5.542682927 the values are               5.780487805 5.682926829 5.341463415 5.414634146 5.43902439         5.414634146 5.341463415 5.926829268


title: How to interpret a differenced series with non zero mean?

text: I have a non-stationary series and a stationary differenced series, the differences have a (small) non-zero mean. I am wondering how we interpret that in terms of the original series. As the differences will revert to and move around the mean (-0.12) the original series will tend to 'drift' downwards? Is this what is meant by a stochastic drift (as opposed to a deterministic drift)?


title: Independence of a linear and a quadratic form

text: How can I prove the following lemma? Let $\mathbf{X}^ \prime$ = $ \left[ X_1 , X_2 , \ldots, X_n \right]$ where $ X_1, X_2, \ldots X_n $ are observations of a random sample from a distribution which is $N \left ( 0,\sigma^2 \right)$. Then let $\mathbf{b}^\prime = \left[b_1,b_2,\ldots,b_n \right]$ be a real nonzero vector and let $\mathbf{A}$ be a real symmetric matrix of order $n$. Then $\mathbf{b ^\prime X}$ and $ \mathbf{X} ^\prime \mathbf{A} \mathbf{X} $ are independent **iff** $\mathbf{b} ^\prime \mathbf{A}=0.$ * * * I know that $\mathbf{b^ \prime X } \sim N(0, \sigma^2 \mathbf{ b^\prime b})$ but I do not see how I could proceed here. My main difficulty lies on the fact that these two variables are very different; had it been two quadratic forms instead, then Craig's theorem would be of use. Any advice? Thank you.


title: Strange definition of RSS for elastic net in scikit-learn (Python)

text: I've been reading the manual about GLM, specifically the one on the elastic net in this link, and they come up with the following function to minimize: $$\frac{1}{2n}||Xw-y||^2_{2}+\alpha \rho ||w||_1+\frac{\alpha(1-\rho)}{2}||w||_2^2,\ \ \ \ \ \ \ \ \ \ \ (1)$$ where $n$ is the number of samples and $w$ are the coefficients of the regression. However, obviating the division by $2n$ (which I find unnecessary), I'm confused by the second part. In Zou & Hastie (2005) the RSS (Residual Sum of Squares) of the elastic net is defined as: $$RSS=||Xw-y||^2_{2}+\lambda \left(\alpha ||w||_1+(1-\alpha)||w||_2^2\right),\ \ \ (2)$$ and I really don't see the similarity between criterions (1) and (2) (like, the relations between the $\rho$ in eq. (1) and $\alpha$ of eq. (2), and the relation between $\alpha$ in eq.(1) and $\lambda$ in eq. (2)): am I missing something here?


title: combining p-values for averaging technical protein quantification replicates in python

text: I have some data about protein expression levels in a cell. For each identified protein, there is an expression level and an associated p-value indicating the confidence that the protein was identified correctly. Two of the samples were technical replicates (i.e. a single sample was split into two parts and analysed separately). I now need to average the expression levels of the two technical replicates, and their corresponding p-values. I thought to use Fisher's method to combine them, which seems to me like the right thing to do. The problem is that I need to convert the result from a $\chi^2$ value into a p-value. Excel has a CHIDIST function which seems to do the trick, but, since I'm likely to be doing this sort of thing a lot in the near future, I thought to write a script to do it to plug into our analysis pipeline. I'm using python to write it, but can't find an equivalent function to CHIDIST. I realize that the process is basically finding the probability of getting the given $\chi^2$ score, so I'd like to check my thinking.   1. Am I doing the right thing combining the p-values?   2. Is Fisher's Method the appropriate process for combining the p-values correct?   3. How many degrees of freedom? I'm thinking 2, as per the wikipedia page.   4. The final p-value is one minus the integral from negative infinity to the $\chi^2$ score over the $\chi^2$ probability density function with one degree of freedom. Thanks


title: Outliers for boxplot

text: If I have a set of data and I want to display it as a boxplot, when there are some outliers in the boxplot, we indicate outliers with a `*`. If there are two outliers having the same value, how to put that in the boxplot? Is it still *?


title: Minimum population size for chi-square test?

text: I'm analyzing data from an experiment in which two independent groups were exposed to an experimental setup without and with treatment. I am testing whether treatment changed the second group's behaviour by performing a chi-square test that compares group 2 (the observed) vs group 1 (the expected). The result indicates there is a significant change in behaviour X² p-value < 0.00014. Now, I am trying to test subgroups to understand better the change, i.e looking at gender, age, and other self reported metrics. My question is, given that group 2 N=40 if I look at age for instance I find people in their 20s and their 60s show significant change but other age groups don't. However people in their 20s N=12 and people in their 60s N=5. Is there a heuristic/rule that says there is a minimum number of people needed to consider a result significant? for instance anything below N=5 cannot be considered significant or anything below N=20% of the population? EDIT: Just to clarify, I am doing a chi-square test of independence (between group 1&2) not a chi-square goodness of fit test. EDIT 2: With this edit I consider the question closed. None of the answers/comments gave me a definitive solution, which I believe says more about the question than the answers. I was hoping for a definitive answer along the lines you need at least 5 ppl or 20% of your sample. It seems the answer is less direct as it is sensitive to many factors. Thanks.


title: Difference Between SPSS Modeler and SPSS Statistics?

text: I guess we can perform most of the operations of modeling in SPSS Statistics software itself, so why there is need of two software packages doing the same thing? If I am wrong please correct my doubts and help me decide when to use what


title: Given a set of points in two dimensional space, how can one design decision function for SVM?

text: Can someone explain me how one goes about designing a SVM decision function? Or point me to resource that discusses a concrete example. **EDIT** For the below example, I can see that the equation $X_2 = 1.5$ separates the classes with maximum margin. But how do I adjust the weights and write equations for hyperplanes in the following form. $$\begin{array}{ll} H_1 : w_0+w_1x_1+w_2x_2 \ge 1 & \text{for}\; Y_i = +1 \\\ H_2 : w_0+w_1x_1+w_2x_2 \le -1 & \text{for}\; Y_i = -1.\end{array} $$ ![enter image description here](http://i.stack.imgur.com/0A0lL.png) I'm trying to get the underlying theory right in 2-D space (as it's easier to visualize) before I think about higher dimensions. **I have worked out solution for this Can someone please confirm if this is correct?** weight vector is (0,-2) and W_0 is 3 $$\begin{array}{ll} H_1 : 3+0x_1-2x_2 \ge 1 & \text{for}\; Y_i = +1 \\\ H_2 : 3+0x_1 -2x_2 \le -1 & \text{for}\; Y_i = -1.\end{array} $$


title: How to plot a mekko chart in Excel?

text: So a question you will have from here! Q: What is Mekko? A. Its like a histogram. The bar chart can show information in two axes. Normally in a bar chart the y axis shows the variation and x shows just categories/labels. Check Mekko Graphics So my question is if there are any open source ways of getting these into Excel :) Update: All the charts on Types of Mekko Charts but specifically the Bar-Mekko as shown here: ![enter image description here](http://i.stack.imgur.com/3yDYm.jpg)


title: What is crisp logic (in the area of classification)?

text: I have a brief question. I am having trouble finding a good reference that explains what crisp logic is. What I think: I have two classification models, a decision tree and a ruleset, which I think are crisp models. Since they say that an instance is either class A or not. I have another classification model, a logistic regression, that is not crisp, since it gives the probability that an instance belongs to this class. I want to mention this in a presentation that I have later today. I tried looking up this subject, but it would be great if somebody could just confirm how I interpreted it. I could also use a good reference to use in my paper. Thanks,


title: Generalized Linear Model

text: How do I identify accuracy data in each class in GLM via cross validation (in Matlab)? Is R in GLM the overall accuracy?


title: How should one define the sample variance for scalar input?

text: I was horrified to find recently that Matlab returns $0$ for the _sample_ variance of a scalar input:               >> var(randn(1),0)   %the '0' here tells var to give sample variance     ans =          0     >> var(randn(1),1)   %the '1' here tells var to give population variance     ans =          0      Somehow, the sample variance is not dividing by $0 = n-1$ in this case. R returns a NaN for a scalar:               > var(rnorm(1,1))     [1] NA      What do you think is a sensible way to define the ~~population~~ sample variance for a scalar? What consequences might there be for returning a zero instead of a NaN? **edit** : from the help for Matlab's `var`:               VAR normalizes Y by N-1 if N>1, where N is the sample size.  This is     an unbiased estimator of the variance of the population from which X is     drawn, as long as X consists of independent, identically distributed     samples. For N=1, Y is normalized by N.           Y = VAR(X,1) normalizes by N and produces the second moment of the     sample about its mean.  VAR(X,0) is the same as VAR(X).      a cryptic comment in the m code for `var states:               if w == 0 && n > 1         % The unbiased estimator: divide by (n-1).  Can't do this         % when n == 0 or 1.         denom = n - 1;     else         % The biased estimator: divide by n.         denom = n; % n==0 => return NaNs, n==1 => return zeros     end      _i.e._ they explicitly choose not to return a `NaN` even when the user requests a sample variance on a scalar. My question is why they should choose to do this, not how. **edit** : I see that I had erroneously asked about how one should define the population variance of a scalar (see strike through line above). This probably caused a lot of confusion.


title: Is there a name for this chart - sort of a cross between a pie chart and a mekko plot

text: Is there a name for this sort of chart below (sourced from New Zealand's Ministry of Business, Innovation and Employment, for whom I work but I was not involved in creating this plot)? It consists of rectangles where the area is proportional to a variable, and resembles a sort of cross between a pie chart, a mosaic plot and a mekko plot. It is perhaps closest to a mekko plot but has the complication that we are not working with columns but a more complex jigsaw. The original looks a little better as there are white borders between the rectangles for each region. Surprisingly, it actually strikes me as not too bad a statistical graphic, although it could be improved through better use of colour mapped to something meaningful. A powerful interactive version showing the US 2011 budget has been used by the New York Times. An interesting challenge is to think of an automatic algorithm to draw one and make it look reasonable too. The rectangles need to be allowed to have different aspect ratios, within an acceptable range. ![enter image description here](http://i.stack.imgur.com/mptL6.png)


title: How to generate Normal variables parts of which are correlated (in R)?

text: My second attempt to explain the question. Start with a vector of numbers V1 of length M. The elements of V1 form a Normal distribution. Take (any) N elements from this vector. Generate a replica of this subset with some random error. Call this V2. Add to V2 elements so it becomes of length M **and** so that all M elements in V2 form a Normal distribution. As a result, V1 and V2 must be two normal random variables with one part that is correlated and one part which is not correlated across the two variables. And the resulting variables must have (the same) predefined mean and standard deviation. My approach has been the following (in R): ` rho<-0.9 #set the correlation level x<-rnorm(1000,0,2) y<-x*rho+sqrt(1-rho^2)*rnorm(1000,0,2) #this generates the correlated parts of X and Y x[1001:1500]<-rnorm(500,0,2) #this is the uncorrelated part of X y[1001:1500]<-rnorm(500,0,2) # this is the uncorrelated part of Y ` The remaining problem is that for now the correlated elements of `x` and `y` are normally distributed which is admissible but should not be necessary. In other words, how to attach elements to any collection of numbers so that the resulting set has a normal distribution? If that would clarify things, in reality the data I am trying to simulate comes in the following form: people assign a fixed number of statements (based on agreement with the statement) into bins which follow a quasi-normal distribution. It is suspected that people score (assign) a subset (but not all) of these statements in a very similar way. I want to simulate such data in order to explore the power of different statistical tools to detect such data structure.


title: Does computing this confidence interval make sense here?

text: I work in a medical company and we often analyze the days patients stay in hospital longer than needed. For example we have these data:               patient      days in hospital     days not necessary           1                  20                      2           2                   5                      0           3                  13                      1           4                   9                      3           5                  22                      0     ..      we now compute the average days per patient that were not necessary for treatment in hospital. In this example:               (2/20 + 0/5 + 1/13 + 3/9 + 0/22) / 5      where 5 is the number of patients in this example. Then we compute a 95% confidence interval for this value (in SPSS). The result is usually an interval of the kind               - 0.9 < \mu < 3.6      Because a negative value makes no sense we change this interval to               0 < \mu < 3.6      For me the whole process feels kind of wrong but I cannot specify it exactly. What would be a better way to do this, or does this make no sense at all? Edit: Some more insight on the data: We are authorized by health insurance companies and get data about patients directly from the hospitals. Our experts, usually doctors, look at each patient's case and make a report about possible unnecessary days of treatment. Reasons can be that an outpatient treatment would be enough or various others reasons. Sometimes patients also want to stay a day or two longer because they feel it is necessary. Of course health insurances have a interest in getting people out of hospital as soon as possible because that's a huge cost factor. Anyway, this is not relevant here. In this report we make a CI like the above example. Of course I cannot provide any real data nor anything that is close to it. The sample size is usually between 40 and 140 patients. There can be no negative days. I think we get negative values for the interval because it's so close to zero and the variance of the assumed normal distribution makes it so. days not necessary is always a fraction of the total days. So for patient 1 in the above example 2 out of 20 days were not necessary. The CI should say in which area the expected value of not-necessary-days/total-days per patient is.


title: Formatting graphs: when is it appropriate to use a fill under a line graph?

text: This is a data visualisation question - I hope it's okay to ask this here. When is it appropriate to use a fill under a line graph for a time series, like the graph below? (which is showing ping times over a day) ![Line graph with fill underneath](http://i.stack.imgur.com/55cTA.png) I guess it's more common to use a plain line, with no fill beneath, but is it OK to use a fill for visual variety? I'd be particularly interested in knowing about any perceptual research on the topic, or any style guides.


title: Linear regression where attribute vectors can be reordered

text: Is there an existing method to solve this kind of regression problem? Given a set of $n$ training data $\\{(x_1,y_1)\ldots(x_n,y_n)\\}$, find $\beta$ and a set of functions $\\{f_1,f_2,\ldots,f_n\\}$ to minimize $\sum_{i}\left(y_i-f_i(x_i)^T\beta\right)^2$ where each $x_i$ is a length $m$ vector, $\beta$ is a length $m$ weight vector, and $f_i$ is a function that permutes the elements of $x_i$. For example, if $x_i=(a,b,c)^T$, then $f_i(x_i)$ is one of the six possible permutations $(a,b,c)^T,(a,c,b)^T,(b,a,c)^T,(b,c,a)^T,(c,a,b)^T,(c,b,a)^T$. Note that without the permutation function $f_i$, the objective function is just the residual sum of squares, $\sum_{i}\left(y_i-x_i^T\beta\right)^2$, and this becomes a typical linear regression problem. This seems like a hard problem. It would be too computationally expensive to solve it by brute force, e.g. perform a normal linear regression on every possible reordering of every possible data point. The only saving grace is that I know a lot of the reorderings will be the same, e.g. the first 200 data points will have the same ordering, the second 200 data points will have the same ordering but possibly different from the first, etc. This kind of problem comes up if there isn't any apparent way to order the attributes in your data. For example, each data point may represent the incomes of a group of 3 people (i.e. $x_i = (income_1,income_2,income_3)$), and each data point comes from a different group. If you don't know a reasonable property by which to associate members across groups, you can't order the incomes in a meaningful way.


title: Brier versus AIC

text: There are four models I'm looking at and trying to evaluate which will give me the best predictions going forward. So I loop through each model and look at the error rate, the Brier Error and the AIC. In some cases the AIC and Brier Error are not consistent (results for model 3 and 4 below). If they give contrasting indications, which one should be paid attention to? Why? Or am I looking at this the wrong way? For those interested, my R code and results are below:               line<- seq(71,79, by=.5)          models<- list("target~total+tot_eit_h_h1+tot_both_h_h1","target~total", "target~total*spread*as.factor(loc)", "target~total+spread*as.factor(loc)")              for(j in 1:length(models)){            for(i in 1:length(line)){              data$target <- ifelse(data$result>line[i], 1,                                ifelse(data$result<line[i], 0,                                       rbinom(dim(data)[1],1,.5)))              model<-glm(models[[j]], data=data, family="binomial")         data$prediction <- predict(model, data, type='response')               brier=with(data, mean((prediction-target)^2))         error_rate <- mean((data$prediction>0.5 & data$target==0) | (data$prediction<0.5 & data$target==1))         aics<-c(aics, model$aic)         briers<-c(briers, brier)         error_rates<-c(error_rates, error_rate)         }       print(paste(mean(error_rates), "error rates"))       print(paste(mean(aics), "AIC"))       print(paste(mean(briers), "BRIER"))       cat("\n")       error_rates<-c()       aics<-c()       briers<-c()          }      RESULTS               [1] "0.381418853255588 error rates"     [1] "1633.50190616553 AIC"     [1] "0.232600213678421 BRIER"          [1] "0.38491739552964 error rates"     [1] "1641.54699258567 AIC"     [1] "0.234847243715002 BRIER"          [1] "0.314091350826045 error rates"     [1] "1451.51883987259 AIC"     [1] "0.200957523053587 BRIER"          [1] "0.324509232264334 error rates"     [1] "1467.88024134523 AIC"     [1] "0.205168360087164 BRIER"


title: Comparing two models

text: I am interested in comparing two logistic regression models. The two models are nested: model 1 contains all predictors, and model 2 contains all predictors except 1. My goal is to test if removing that one variables makes the modelling significantly worse. I understand that I can use the likelihood ratio test (LRT) for the comparison. But I'm wondering if I can use the Akaike Information Criterion (AIC) to compare the two models. Use of AIC seems more elegant, and gives stronger results in the study that I am looking at. On the other hand, I see that LRT is specifically designed for nested models. Is it acceptable to use AIC to compare the two nested models, or is the LRT test more accurate than AIC when the models are nested? Thank you, Fijoy


title: Interpreting plot of PCA results (from 3 to 2 dimensions)

text: I'm having trouble understanding how to interpret/explain the end result of dimensionality reduction via PCA. Namely, I've attempted to code up a simple example in R but can't really say what happened. I started with three vectors (x,y,z) where x and y are highly correlated and z is veeery loosely related to x. After finding the top 2 principal components and multiplying back by the original data, I have what I believe to be the "final result"... but the plot of the new 2-D data appears entirely uncorrelated. What am I seeing? What are the axes? I expected to see something that resembled the data on the x-y plane. I have plots of the mean-subtracted data in 3d and the reduced plot but not enough points to post them.


title: In a linear regression, why is the mean of all $y_i$ equal to 0?

text: ![enter image description here](http://i.stack.imgur.com/JQCjR.png) I don't understand why the mean for all $y_i$s must equal $0$. What property/properties does the above assertion rely on? I had thought that the mean of $y_i=E(y_i)$, so: $$ E(y_1)=\alpha $$ $$ E(y_2)=\alpha+\beta $$ $$ E(y_3)=\alpha+\gamma $$ $$ E(y_4)=\alpha+\beta+\gamma $$ Since:   1. $\alpha$, $\beta$ and $\gamma$ are constants   2. The mean of $\epsilon_i$'s are 0, as provided by the question. Am I wrong?


title: Comparing effective password spaces in textual and graphic passwords

text: I've read a lot of theory about graphic passwords and how people remember them. I wanted to test the theory, so I created an application that allowed people to create textual and graphic passwords and recall them after some distraction tasks. The graphic passwords were a combination of choosing an image from a grid and then positioning the chosen images on a canvas. In both cases, the only restriction was that the passwords weren't blank. I've already had people run through the application so I have some data to work with. I want to test two hypotheses: **Hypothesis 1** : People will choose passwords with equal effective strengths in both methods **Hypothesis 2** : People will be able to recall their passwords equally well I'm looking for any high-level or algorithmic advice knowledgeable people can give on testing these hypotheses. My thoughts after talking to my academic supervisor and a mathematician friend are below: Hypothesis 1:   * We should follow an established procedure for scoring textual passwords and attempt to generate a strength score that follows the same distribution for the graphical passwords. We can then compare the distributions based on knowledge of the theoretical password spaces.   * I'm looking for a way to compare 2 sets of semi-scientific scores representing the effective password spaces that are broken down into groups of passwords generated per person, whilst having a more scientific understanding of the theoretical password spaces. Hypothesis 2:   * Although we have data about the Levenshtein distances of the passwords, the reality is that passwords are generally either right or wrong based on whether they are an exact match. With this in mind, we also have produced binary data for each password attempt. Unfortunately I don't know how to compare the binary data well other than taking a simple average although I've heard it's possible with binomial distribution.   * Since each person recalled 3 passwords 2 times, I'm looking for a statistical method that takes into account the the person, the amount of distraction they've had (i.e. 5 minutes or 15 minutes) and whether the entry was an exact match or not to give some indication of how likely people are to get a password wrong after 5 or 15 minutes of distraction. I will also need to know if this difference is significant.


title: Algorithm for minimization of sum of squares in regression packages

text: I'd like to know what technique or techniques are used by regression packages, (in particular the `lm` function of R) to minimize the sum of squares. Does it use gradient descent? Thanking in anticipation.


title: What can be used instead of confirmatory factor analysis?

text: Can you please share your thoughts of the best method to deal with the following issue. In my study I'm using a well-known tool (questionnaire) which was validated in several settings but developed in 80's, and not direcly in the area I'm studing.I would like to "update" some of the items and add questiones relevant to my field.Unfortunately, because of the sample size requirements I can't use CFA. Is there any other method of validating what would be a "updated" version of the instrument?


title: Why even my training data failed during the prediction of libsvm

text: Currently I'm using libsvm for my one class classification problem. I have 10 samples in my training set, 5 samples in my testing set, both of my training and testing set is scaled by svm_scale, then I use smv_train with the scaled training set for training, and svm_pridict with the scaled testing set for testing. But this problem occurs which makes me very confusing. In the result, why even my training sample fails....               Accuracy = 26.666666667% (4/15) (classification)


title: Method for simulation for random variables with a particular probability density function

text: Suppose one has an expression for the probability density of a random variable- how does one simulate for the particular random variable. I understand that this isn't an issue for probability density functions like Gaussian which have a MATLAB implementations(say wgn()). I have a probability density function which doesn't have a Matlab implementation for realization. What does one do ? Are there any particular numerical methods for say turning a rand() into a particular probability density function.


title: Sample divided in three equal categories of individuals in logit regression

text: I am analysing a dataset composed by three groups of individuals: unemployed, unstable workers and regular workers. These three groups are numerically equal (400 individuals in each category). I am running a logit analysis to test the effect of occupational status (regular worker, unemployed and unstable worker) on a bivariate dependent variable (dummy 0 or 1). I use as regressors (independent variables) a categorical variables 1-3 to indicate occupational status and others regressors (educational level, gender etc.etc.). Is it correct to use occupational status as a regressor even if my dataset is not a random sample, but is built "artificially" to have three equal groups of individuals? I have terrible doubts on the procedure! If anybody have suggestions I would be glad. Thanks,Lara


title: MLEs from glmer {lme4} in R

text: Numerically deriving the MLEs of GLMM is difficult and, in practice, I know, we should not use the brute force optimization (e.g., using `optim` in a simple way). But for my own educational purpose, I want to try it to make sure if I correctly understand the model (see the code below). I found that I always get inconsistent results from the result from `glmer`). In particular, even if I use the MLEs from `glmer` as initial values, according to the likelihood function I wrote (`negloglik`), they are not MLEs (`opt1$value` is smaller than `opt2`). I think two potential reasons are: (1) `negloglik` is not written well so that too much numerical error in it, and (2) the model specification is wrong. For the model specification, the intended model is \begin{equation} L=\prod_{i=1}^{n} \left(\int_{-\infty}^{\infty}f(y_i|N,a,b,r_{i})g(r_{i}|s)dr_{i}\right) \end{equation} where _f_ is binomial pmf and _g_ is normal pdf. I am trying to estimate _a_ , _b_ , and _s_. In particular, I want to know if the model specification is wrong, what the correct specification is.               p <- function(x,a,b) exp(a+b*x)/(1+exp(a+b*x))          a <- -4  # fixed effect (intercept)     b <- 1   # fixed effect (slope)     s <- 1.5 # random effect (intercept)     N <- 8     x <- rep(2:6,each=20)     n <- length(x)      id <- 1:n     r  <- rnorm(n,0,s)      y  <- rbinom(n, N, prob = p(x,a+r,b))               negloglik <- function(p,x,y,N){       a <- p[1]       b <- p[2]       s <- p[3]            Q <- 100  # Inf does not work well       L_i <- function(r,x,y){         dbinom(y,size=N,prob=p(x,a+r,b))*dnorm(r,0,s)       }            -sum(log(apply(cbind(y,x),1,function(x) integrate(L_i,lower=-Q,upper=Q,x=x[2],y=x[1],rel.tol=1e-14)$value)))     }          library(lme4)     (model <- glmer(cbind(y,N-y)~x+(1|id),family=binomial))          opt0 <- optim(c(fixef(model),sqrt(VarCorr(model)$id[1])),negloglik,x=x,y=y,N=N,control=list(reltol=1e-50,maxit=10000))      opt1 <- negloglik(c(fixef(model),sqrt(VarCorr(model)$id[1])),x=x,y=y,N=N)     opt0$value  # negative loglikelihood from optim     opt1        # negative loglikelihood using glmer generated parameters     -logLik(model)==opt1 # but these are substantially different...      **A simpler example** To reduce the possibility of having large numerical error, I created a simpler example.               y <- c(0,3)     N <- c(8,8)     id <- 1:length(y)          negloglik <- function(p,y,N){       a <- p[1]       s <- p[2]       Q <- 100  # Inf does not work well       L_i <- function(r,y){         dbinom(y,size=N,prob=exp(a+r)/(1+exp(a+r)))*dnorm(r,0,s)       }       -sum(log(sapply(y,function(x) integrate(L_i,lower=-Q,upper=Q,y=x,rel.tol=1e-14)$value)))     }          library(lme4)     (model <- glmer(cbind(y,N-y)~1+(1|id),family=binomial))     MLE.glmer <- c(fixef(model),sqrt(VarCorr(model)$id[1]))     opt0 <- optim(MLE.glmer,negloglik,y=y,N=N,control=list(reltol=1e-50,maxit=10000))      MLE.optim <- opt0$par     MLE.glmer # MLEs from glmer     MLE.optim # MLEs from optim          L_i <- function(r,y,N,a,s) dbinom(y,size=N,prob=exp(a+r)/(1+exp(a+r)))*dnorm(r,0,s)          L1 <- integrate(L_i,lower=-100,upper=100,y=y[1],N=N[1],a=MLE.glmer[1],s=MLE.glmer[2],rel.tol=1e-10)$value     L2 <- integrate(L_i,lower=-100,upper=100,y=y[2],N=N[2],a=MLE.glmer[1],s=MLE.glmer[2],rel.tol=1e-10)$value          (log(L1)+log(L2)) # loglikelihood (manual computation)     logLik(model)     # loglikelihood from glmer


title: Data mining classification competition

text: I'm currently taking a data mining class, and for one our projects we're required to predict the class label for an unknown data set by first building a classifier on a training data set which already provides the class label. We're only required to get an accuracy of 80% to get a full mark on the assignment. I have already achieved this using the J48 Decision Tree algorithm (acc=84.08%). There is also an ongoing competition on who can get the highest accuracy (determined by a Judge system we can't see). I have two questions:   1. How can I use an ensemble method with to do this   2. Is there a way to optimize the parameters for each classifier? * * *                import java.io.*;      import weka.core.Instances;      import weka.filters.Filter;      import weka.filters.unsupervised.attribute.*;      import weka.classifiers.trees.*;      import weka.classifiers.Evaluation;          public class CompClassifier {     public static FileOutputStream Output;     public static PrintStream file;          public static void main(String[] args) throws Exception {         // load training data         weka.core.Instances training_data = new weka.core.Instances(new            java.io.FileReader("/Users//Weka/training.arff"));              //load test data         weka.core.Instances test_data = new  weka.core.Instances(new      java.io.FileReader("/Users//Weka/unknown.arff"));              //Clean up training data         ReplaceMissingValues replace = new ReplaceMissingValues();         replace.setInputFormat(training_data);         Instances training_data_filter1 = Filter.useFilter(training_data, replace);               //Normalize training data         Normalize norm = new Normalize();         norm.setInputFormat(training_data_filter1);     Instances processed_training_data = Filter.useFilter(training_data_filter1, norm);              //Set class attribute for pre-processed training data      processed_training_data.setClassIndex(processed_training_data.numAttributes() - 1);                  //output to file     Output = new  FileOutputStream("/Users//Desktop/CLASSIFICATION/test.txt");         file = new PrintStream(Output);              //build classifier         J48 tree = new J48();         tree.buildClassifier(processed_training_data);              //Clean up test data         replace.setInputFormat(test_data);         Instances test_data_filter1 = Filter.useFilter(test_data, replace);               //Normalize test data         norm.setInputFormat(training_data_filter1);         Instances processed_test_data = Filter.useFilter(test_data_filter1, norm);              //Set class attribute for pre-processed training data         processed_test_data.setClassIndex(processed_test_data.numAttributes() - 1);              //int num_correct=0;         for (int i = 0; i < processed_test_data.numInstances(); i++) {             weka.core.Instance currentInst = processed_test_data.instance(i);             int predictedClass = (int) tree.classifyInstance(currentInst);             System.out.println(predictedClass);             file.println("O"+ predictedClass);         }               }


title: Log-normal distribution versus chi-square distribution for comparing RMSE of nonlinear fits

text: I have parameter estimates fitting a particular nonlinear model for thousands of experimental cycles. My goal is to find a nice way to tag those cycles which didn't fit the model very well. Currently I am comparing RMSE values for this. I have a huge list of RMSE values and am looking for a cutoff RMSE past which to flag it. I plotted the histogram of RMSEs and found that they fit a log-normal distribution. However, I was informed that the proper distribution to use in this case is $\chi^2$. For my purposes the distribution I use doesn't seem to matter as long as it fits well, but since I was curious, I fitted a $\gamma$ distribution (which wikipedia informed me $\chi^2$ was a special case) to the RMSE histogram, but the parameters didn't seem to fit quite as well. Is the $\gamma$ distribution really that more meaningful? Also, if you have any alternative suggestions for approaching my problem I would love to hear them.


title: Non-parametric test with several confounding factors

text: I am pretty new to statistical testing with confounding variables, so please excuse if there will be some gross misunderstandings or basic questions. I want to test for systematic differences in species abundance ($y$) between several types of environments ($x$). The way I would usually approach is with a Kruskal-Wallis test. However, these different environments are not all in the same region. So I would like to use the region factor as block ($b$). The differences between samples in $y$ and $b$ are subtle, i.e. no obvious clustering can be seen corresponding to $x$ or $b$. For this I was using now the package coin for R, that is               independence_test(y ~ x | b, ...)      From here I have now a few theoretical questions:   1. How do I visualize this? I would prefer a boxplot, but how to include the blocking factor?   2. I would also like to test correlations (spearman) using a block design, what package would you suggest? Again, how does one visualize this (e.g. scatter plot, but with block factor)?    3. I would like to control for a continuous variable ($z$) at later stages. As I understand, this can not be done with "block"ing factor in `coin`? How would I do this correctly?    4. The last stage would be to control for **both** blocking factor b and continous variable z in a Kruskal-Wallis test. How can this be done? I would be grateful for any help on this topic.


title: Classification datasets with categorical features and a lot of classes

text: As stated in the title I would like to find some classification datasets with $c \ge 5$ classes and one categorical feature or more. It is even fine if there are some continuous features but I would like at least a categorical one. I have some genetic datasets where features take values in $\\{ G, T, A, C\\}$ but unfortunately the class is binary. I have a method that should improve the accuracy performance on those dataset. I used some of UCI datasets but I need more of them to show a significant $p$-value that states the improvement of performances.


title: Using sentiment lexicons or all words processing for sentiment analysis?

text: I am learning sentiment analysis to apply it to twitter real time data to predict user's mood. I ponder about using which alternative way to do that data mining job.   1. Use all words to process and create a model by SVM or such classifier   2. Get vast lexicon list that is obtained form a outsource and use it as a vocabulary and just use these lexicons as a features on classifier.   3. Just use lexicons and instead of using a model emerged by a classifier, take the real time value and look at frequencies of lexicons at the tweets and rank. (this is not the way for me but I want to see the public wisdom about)   4. Or use any other method. As a sub question, If predicting mood will be complicated for me I will just take care of the binary classification as positive or negative feelings. Which method will work best in that case?


> **Possible Duplicate:**   >  Confusion regarding a probability problem A certain auditorium has 30 rows of seats. Row 1 has 11 seats, while Row 2 has 12 seats, Row 3 has 13 seats, and so on to the back of the auditorium where Row 30 has 40 seats. A door prize is to be given away by randomly selecting a row (with equal probability of selecting any of the 30 rows) and then randomly selecting a seat within that row (with each seat in the row equally likely to be selected). Now, 1) Find the probability that Seat 15 was selected given that row 20 was selected ? 2) Find the probability that Row 20 was selected given that Seat 15 was selected ? To answer the first, given that Row 20 was selected, there are 30 possible seats in Row 20 that are equally likely to be selected. Hence Pr(Seat 15 | Row 20) = 1/30. The same kind of argument can be given to answer the second : given that Seat 15 was selected, there are 30 possible rows that are equally likely to be selected. Hence Pr(Row 20 | Seat 15) = 1/30. Now, it turns out that the first answer is correct whereas the second answer is incorrect. My question is where am I making mistakes in computing the second answer ?


title: What is the best way to discretize a 1D continuous random variable?

text: Say I have a 1-dimensional continuous random variable $X$, with PDF $f(X)$, CDF $F(X)$ and inverse CDF $F^{-1}$. What is the best way to discretize $X$? To keep things clear, let $Y$ denote the discretized version of $X$, and let $g(Y)$, $G(Y)$ and $G^{-1}$ refer to the PMF, CDF and inverse CDF of $Y$. Ideally, I would like for $E[X] = E[Y]$ and $Var(X) = Var(Y)$. However, I understand if these properties cannot hold for finite $N$ due to discretization error. Some thoughts:   * Should I just discretize $X$ into $N$ equally spaced points? Or should I discretize more at regions of high probability, by first discretizing a Uniform $[0,1]$ into equidistant points $u_1,u_2...u_N$ and then generating discrete values $Y_1 = F^{-1}(u_1), Y_2 = F^{-1}(u_2)$ and $Y_N = F^{-1}(u_N)$   * Say I have discretized values $Y_1, Y_2... Y_N$. What is the best way to create $g(Y)$, $G(Y)$ and $G^{-1}$ from these values? My thoughts here are that each $Y_k$ corresponds to an interval of $X$ with width $w_k$, and that the $Y_k$ should be the "midpoints" of these intervals. That is, $g(Y_k) = p(Y_k - w_k < X < Y_k + w_k) = F(Y_k + w_k) - F(Y_k - w_k)$.


title: Are mediation analyses inherently causal?

text: I am interested in testing a simple mediation model with one IV, one DV, and one mediator. The indirect effect is significant as tested by the Preacher and Hayes SPSS macro, which suggests the mediator does serve to statistically mediate the relationship. When reading about mediation I have read things such as "Note that a mediational model is a causal model." - David Kenny. I can certainly appreciate the use of mediation models as causal models, and indeed, if a model is theoretically sound, I can see this as very useful. In my model, however, the mediator (a trait considered to be a diathesis for anxiety disorders) is not caused by the independent variable (symptoms of an anxiety disorder). Rather, the mediator and independent variables are related, and I believe the association between the independent variable and the dependent variable can be explained largely by variance between the IV- mediator-DV. In essence I am trying to demonstrate that previous reports of the IV-DV relationship can be explained by a related mediator that is not caused by the IV. Mediation is useful in this case because it explains how the IV-DV relationship can be statistically explained by the IV-Mediator-DV relationship. My problem is the question of causation. Could a review come back and tell us that the mediation is not appropriate because the IV does not in fact cause the mediator (which I would have never argued in the first place)? Does this make sense? Any feedback on this matter would be greatly appreciated! **Edit** : What I mean to say is that X is correlated with Y not because it causes Y, but because Z causes Y (partially) and because X and Z are highly correlated. A bit confusing, but that is it. The causal relationships in this instance are not really in question and this manuscript is not so much about causation. I simply seek to demonstrate that variance between X and Y can be explained by variance between Z and Y. So basically, that X is correlated indirectly to Y through Z (the "mediator" in this case).


title: Distribution that describes the difference between negative binomial distributed variables?

text: A Skellam Distribution describes the difference between two variables that have Poisson distributions. Is there a similar distribution that describes the difference between variables that follow negative binomial distributions? My data is produced by a Poisson process, but includes a fair amount of noise, leading to overdispersion in the distribution. Thus, modeling the data with a negative binomial (NB) distribution works well. If I want to model the difference between two of these NB data sets, what are my options? If it helps, assume similar means and variance for the two sets.


title: Distributing values according to a given continuous probability distribution

text: I have a series of continuous probability distributions and the intervals for each distribution of possible values. Now I want to use the distributions to distribute values of a random variable among agents. I have thought of one way to do this, and I'm curious if this is an acceptable method. This is how I would do it:               Let P(A<X<B) be a given pdf, where A,B are min and max values for random      variable X.     Divide the spread between A and B into an arbitrarily large amount of 'chunks'      or ranges.     For each chunk, represented as a range, say between a and b, make a table      where P(a<X<b) relates to the range a,b.     For each agent receiving a value of random variable X, use a random number      generator, look up the result in the table, give the agent a value in the      range (a,b).      Assuming that this is an acceptable technique, I'm still curious as to what an acceptable amount of granularity for the 'chunks' would be. I hope this wan't too confusing and any feedback would be much appreciated. edit: All of my functions take the form of f(x) = a/(b*x)+c over some interval. They are all easily integrable.


title: Hausman's test for all $\beta$s – comparing FE vs RE models

text: I fit several two level models in SAS using `PROC MIXED`: an empty model with multilevel structure (null), a model with a level 2 covariate (partial model), and a model with level 1 and level 2 covariates (full model). SAS estimated the random intercepts for all models, and the random effect covariate in the partial model, but in the full model the estimates for the random effect were all 0. My colleague suggested this because the covariate is a poor predictor in the presence of the level 1 covariates, that within group variation may be larger than between group variation, and that I should do a Hausman's test to see whether the multilevel structure is needed. As I understand it, Hausman's test is $W = \frac{(\hat{\beta}_{FE}^\star - \hat{\beta}_{RE}^\star)}{Var(\hat{\beta}_{FE}^\star) - Var(\hat{\beta}_{RE}^\star)}$ which approximates a $\chi^2$ distribution with 1 degree of freedom, where $\hat{\beta}_{FE}^\star$ is the estimated beta coefficient from the fixed effects model, $\hat{\beta}_{RE}^\star$ is the estimated beta coefficient from the random effects model, and $Var(\hat{\beta}_{FE}^\star)$ and $Var(\hat{\beta}_{RE}^\star)$ are their variances. The relationship to the chi-square distribution is intuitive and it is easy to see how to implement this for a single coefficient. My question is whether I should calculate $W$ for each $\beta$ in my model (i.e. level 1 covariates and the level 2 covariate) or just the random effects covariate.


title: What are the general strategies in creating a Probabilistic Graphical Model?

text: While there is lot of theory and probability in the background to understand, I wanted to know if there are any resources/quick pointers on what to consider while modeling a problem using Bayesian Networks. Specifically, I had the following questions: - what kind of priors to use? - what do we learn from the data? and how to do prediction? Is there a blog/document that explains these from the ground-up with simple use-cases? (eg: modeling user-click behavior, weather forecasting, diagnosis etc). I am looking for examples which are a little detailed beyond toy- examples (a bit close to real life situations) used in courses. Any help is appreciated.


title: How can I quantify the phase/amplitude relationship between two signals?

text: I have two signals `A` and `B`. I want to show that high amplitude events in `B` are phase locked to oscillations in signal `A`. I have already identified candidate events signal `B`. I estimate `phase_of_A` using the `angle` of the `hilbert` transform of `A`, and I estimate the `envelope_of_B` using the `abs` of the `hilbert` transform of `B`. A scatter plot of the joint distribution shows a clear relationship between the two variables. Additionally if I randomly re-order the `phase_of_A` variable the structure is disrupted. My current thinking is use Monte-Carlo methods to demonstrate that the real data is more structured than the shuffled data, I'm just not sure what parameter/statistic to compute on the joint distribution. ![](http://i.imgur.com/lS7LT.png)


title: Bootstrapping - with population data (dependent variable is not normally distributed)

text: Thank you in advance for help! I am conducting a study using General Linear Modeling on the distribution of financial aid at a college. I am not looking to project my findings onto a larger population. I only want to make statements/conclusions about this one school. The dataset has over 1000 students, each of whom were offered an scholarship. My dependent variable (the scholarship award) is not normally distributed. The residuals of the model show non normal distribution and I have tried to take a log of my dependent variable... still not normally distributed. Bootstrapping was suggested to me as a possible solution. Can someone help me understand if bootstrapping might change the nature of my data? Can I interpret the output of the GLM parameter estimates the same way? Thank you very much.


title: (Nominal) raters with no gold standard

text: A friend of mine took a document and broke it up into parts, then asked 5 subject matter experts to classify each part into nominal category A, B, C, D, or E. (I'm not sure yet, but D may be "All of the Above" and E may be "None of the Above".) So we have around 200 sections of the document, each with five opinions as to what category the section falls into. The friend also classified each section herself, so we also have that sixth set of data. She believes that there are particular patterns to the sections, which has not been generally noticed by others in her field. Questions:   1. I've argued that we should not include her classification in any analysis, because it's her hypothesis -- the patterns she sees -- that we're looking to confirm, so we should only use the five other experts. Does this seem correct?   2. The test was originally given to seven experts, but two of them did not follow the instructions to pick a single category for each section. They gave answers like "A or B". Between the two of them, they chose multiple categories for 30% of the sections, so we chose not to use their answers. (My friend had also decided before giving the questionnaires to the subjects that any with more than 10% errors would not be used.) Does this seem like the correct choice?   3. Given the first two decisions, we decided that a section's category would be "strongly agreed to" if 4/5 or 5/5 of the experts put it into the same category. Does this seem reasonable, or should some kind of calculation be used to determine a level for strong agreement for each section?   4. Assuming that we decide we have enough sections with strongly-agreed-to categories to look for patterns, does anyone have pointers on the technique to use? For example, we might hypothesize that there is a pattern of sections of A followed by B followed by C (with perhaps E's or maybe D's thrown in). Is there a field or technique for this, or does it boil down to pattern matching and combinatorics?   5. I've created various contingency tables of experts/categories/sections/counts, done Cohen's Kappa, CA, MCA, etc, and it's all interesting, but perhaps the bottom line is that only about 1/3 of the sections have strong agreement (4/5 or 5/5 experts agree on the category), and I'm not sure that's enough coverage of the document to actually proceed to the next step and look for patterns. Any thoughts or suggestions about fields/techniques that might be helpful? (For example, I've briefly looked at IRT but it appears to not be useful without a gold standard and with nominal categorical data rather than ordinal.)


title: Expected value of modified geometric distribution

text: I am trying to find the expected value of $X$, where $X$ is the number of orders a customer will make in a lifetime. Assuming that there is a $p=.1$ chance of the customer placing an initial order, and then (given that the customer places that initial order) a $p=.9$ chance that the customer places an order after that (and each additional time the customer places an order after that as well. NB it is dependent — i.e. if the customer doesn't place an order, then he won't place any additional orders). So I have it set up like this: * * * $X = \\#$ of orders a customer places   $O = \\#$ of orders a customer places after the first one $E[X] = 0.1(E[O] + 1) + 0.9(0) = 0.1 + 0.1 \cdot E[O]$ $P(O = 0)$ = $0.1$   $P(O = n)$ = $0.1 \cdot 0.9^n$ for $n > 0$ So, $$ E[O] = \sum_{k=0}^\infty k \cdot P(O=k) = 0.1 \sum_{k=0}^\infty k \cdot 0.9^k = 9 \>, $$ leaving us with $E[X] = 0.1 \cdot (9 + 1) = 1$. * * * EDITED: using a shifted geometric distribution with $p=0.1$ for $O$ Is that correct? If not, where am I going wrong?


title: What are some software packages for doing panel aggregated logistic regression with no random or fixed effect intercept (not SAS)?

text: I am interested in doing Panel Aggregated Logistic Regression (aka Panel Binomial Logistic regression) without random intercept and I was wondering if anyone could suggest some existing packages. There are many packages that does aggregate logistic regression without random intercept but they are not panel-ready. Just to be clear, without panel, what I want to do can be done in R with GLM as such `model=glm(cbind(success,trial-success)~ x1+x2-1,family=binomial('logit'),data=d)` I am not certain but i think `xtmelogit` cannot do without random effect. `xtlogit` cannot handle binomial logit. and unblock + `xtlogit` is just not the same. I also know SAS can do this but do not own a copy. I know how to do this with some already existing packages (xtgee) but wanted to compare all the other existing packages across many different softwares. Thanks! Edit: Turns out -init- option on gllamm can do this as well


title: What is the correct ANOVA for repeated measures with multiple measurements for each within-subject factor?

text: In Andy Field's example for one-way ANOVA with repeated measures, he measures "retching time" for eight different celebrities eating four different gross foods. He tests things like sphericity, runs one-way ANOVA and post hoc tests. I have a similar set up but with multiple measures for each subject-factor combination. Imagine in Andy's example that for each celebrity-food pair I'd measured 50 different retching times. I've seen it suggested, on this site even, to just take the average across these 50 trials. Is this "legal". Shouldn't I have to verify that the variance for each celebrity-food pair is small? If so how? Should these trials go into one large analysis or do I just need to preprocess my data? As a last cry for help. How is the best accomplished with a program like SPSS?


title: Fuzzy C-Means / Latent Dirichlet Allocation

text: I have to compare the results of document clustering via latent Dirichlet allocation and fuzzy _c_ -means. How can I do this? I've got the option to compare the probabilities the documents are assigned to the cluster/topic. Are there other ways I could compare?


title: Differences between different Geometric Brownian Motion

text: I'm currently studying brownian motion and came across two kinds of geometric brownian motion. http://homepage.univie.ac.at/kujtim.avdiu/dateien/BrownianMotion.pdf (page 14) http://www.math.unl.edu/~sdunbar1/MathematicalFinance/Lessons/StochasticCalculus/GeometricBrownianMotion/geometricbrownian.pdf (page 2) I can't see the difference between the two formulas. I prefer the 2nd one since it seems simpler but both give different kind of expectation and variance and this is making it confusing for me to understand. Can anyone explain the difference between the two?


title: Coupling time series information from sources with multiple spatial resolutions/scales

text: I have many satellite raster images available from different sensors. From these, the coarser ones have a very abundant temporal resolution. The medium resolution rasters tend to have less acquisition dates but still some degree of information is available. The finer resolution ones have a very low temporal resolution, spanning from 2 to 6 observed dates in under two years. I was wondering if anyone knows of any efforts to study this type of multi-scale time series in any way? I would be interested in predicting future values at the finer scales using the information available from the coarser ones. It makes sense to me that the data must be related (yes the images cover the same regions) but I have no idea how to start coupling this information in a predictive model.


title: Is Kernel PCA with linear kernel equivalent to standard PCA?

text: If in kernel PCA I choose a linear kernel $K(\mathbf{x},\mathbf{y}) = \mathbf x^\top \mathbf y$, is the result going to be different from the ordinary linear PCA? Are the solutions fundamentally different or does some well defined relation exist?


title: How to estimate parameter for Biased Regression Estimator

text: Al-Gereari 2012 introduces the Modified Marquardt Estimator (MME) as: $\hat{\beta^{(r)}(m,k,J)}=T_r[\Delta_r(k)]^{-1}T_r'+m T_{p-r}[\Delta_{p-r}(k)]^{-1}T_{p-r}')(X'Y+kJ)$ where $X'X=\begin{pmatrix} T_r, T_{p-r}\end{pmatrix} \begin{pmatrix} \Delta_r & 0 \\\ 0 & \Delta_{p-r} \end{pmatrix}\begin{pmatrix} T_r'\\\ T_{p-r}'\end{pmatrix} $ corresponding to the linear estimator: $Y=X\beta+e$ where $Y$ is n*1 vector of responses   $X$ is n*p pbserved matrix of regressor variables   $\beta$ is a p*1 vector of unknown parameters e is a n*1 vector distrubuted $N(0,\sigma^2 I)$ $k$ is the ridge parameter to deal with ill conditioning   $J$ represents the prior information that the model shrinks towards as in Crouse et al 1995 However Al-Gereari does not give a method for estimating $m$. This estimator is biased (unless m=1) and so like all biased estimators should be compared on a minium mean square error basis MMSE basis i.e.: $(\hat{\beta^{(r)}(m,k,J)}-\beta)(\hat{\beta^{(r)}(m,k,J)}-\beta)'=Var(\hat{\beta^{(r)}(m,k,J)})+Bias(\hat{\beta^{(r)}(m,k,J)})^2$ I was proposing to estimate m by deriving the expressions for Var() and Bias() and then take take the derivative w.r.t m to find the value of m that minimises the MMSE for this estimator. Is this valid approach to estimate $m$? I then intent to compare this model to other biased estimators on an MMSE basis. Do I introduce any bias because the estimation of $m$ is also based on an MMSE criterion? I don't think so but want to check.


title: Are the degrees of freedom for Welch's test always less than the DF of the pooled test?

text: I am teaching a course on basics statistics, and we are doing the t-test for two independent samples with unequal variances (Welch test). In the examples I have seen, the adjusted degrees of freedom used by the Welch test are always less than or equal to $n_1+n_2-2$. Is this always the case? Does the Welch test always reduce (or leave unchanged) the degrees of freedom of the pooled (equal variances) t-test? And on the same subject, if the sample standard deviations are equal, do the DF's of the Welch test reduce to $n_1+n_2-2$? I looked at the formula, but the algebra got messy.


title: How to symmetrically adjust data by removing certain value

text: Say I have a data set in the following format: {a:200, b: 100; c:400} The sum of a, b and c equals 700. Now, I have another value which is supposed to be the actual total say 500. How do I adjust this difference of 200 in the initial data set (the values for a, b and c must be proportionately decreased summing 200 - something like a weighted difference)?


title: Add new distribution to JAGS or OpenBUGS?

text: I would like model some variables using a distribution that is not neither JAGS' nor OpenBUGS core distributions. Do you know how can I implement it? Thanks.


title: Correlation between two sets of time series data

text: I have an assignment to forecast snow depth at a certain location at a specific date. My approach to forecast this is to basically take previous years snowfall data and compare it to this years data so far. I need to figure out which snowfall pattern over the last 50 years is most correlated with this years pattern so far. I have matlab at my disposal and I have all the time series data but I don't know the best way to measure correlation between a certain years data and the current years data statistically. I have thought about doing an absolute value difference between this years daily data and a given previous years, then averaging it to find the average difference. Whichever one has the smallest average difference will be the year I choose to base the prediction off of. I don't think this is the most accurate method though. I am not very knowledgeable in statistics and want to be able to show that I used a statistically accurate method to measure correlation when I turn in my report.


title: Are these statistics sufficient?

text: Question (Casella and Berger 6.5): Let $X_1 \ldots X_n$ be independent random variables with pdfs: $f(x_i|\theta)= \begin{cases} \frac{1}{2i\theta}, & -i(\theta - 1)<x_i<i(\theta+1) \\\ 0, & \textrm{otherwise} \end{cases}$ Noting the inequality can be rewritten: $-\theta<\frac{x_i}{i} - 1<\theta$ Then we can use indicator variables to write the likelihood function in one line: $L(\theta| D )=\Pi_{i=1}^n f(x_i|\theta)= \left[ \left( \frac{1}{2\theta} \right)^n \Pi_{i=1}^n I_{[\frac{\min(x_i)}{i} - 1, \infty)}(x_i)I_{(-\infty, \frac{\max(x_i)}{i} -1]}(x_i) \right] \left[ \Pi_{i=1}^n \frac{1}{i} \right] $ This latter expression being of the form required for the factorisation theorem: $L(\theta|D)=g(T(D)|\theta)h(D)$ With $h(D)=\Pi_{i=1}^n \frac{1}{i}$ and $g(T(D)|\theta)$ the rest. The sufficient statistics being $\left( \frac{\min(x_i)}{i} , \frac{\max(x_i)}{i} \right)$ * * * I guess I'm wondering if I've missed any requirements to invoke the factorisation theorem. Thoughts appreciated. * * * UPDATE: $L(\theta| D )=\Pi_{i=1}^n f(x_i|\theta)= \left[ \left( \frac{1}{2\theta} \right)^n \right] \left[ \Pi_{i=1}^n \frac{1}{i} I_{[\min(\frac{x_i}{i}), \max(\frac{x_i}{i})]}(\frac{x_i}{i}) \right] $ I think this corrects the errors in set up of the indicator variable. Also, I realise it was in the wrong group, it should be part of $h(D)$ and not $g(T(D)|\theta)$ as it involves the individual $x_i$ and not $\theta$.


title: Faster alternative to multivariate LOESS?

text: I want to make predictions by creating a smooth response surface to 2 variables. I get good results using R's loess() function, but with 10 million observations, it is far too slow. Are there any other smoothing algorithms that would provide similar results at a lower computational cost?


title: LDA topic models for various and unknown domains

text: We have a large collection (1-2M) of documents of various domains (politics, design, programming, etc.). And lets assume that we don't know the exact number of domains. And our goal is to build relevant topic models for these domains. For example I want to find _N_ topics in each domain. For politics it could be: foreign policy, government, economic policy, etc. For design domain: web design, architecture, etc. And so on.      If I understand correctly, when using a LDA it is necessary to form a training set for specified domain (e.g. politics). But problem arises if there is no information about domain. How can I form a training set in that case? I was thinking about building index, extracting documents, that contain significant words (defined by TF-IDF). But it seems that it is not the best solution.      If instead of LDA we can use HDP and handle training set compiled with random documents. But there arises another problem - a huge vocabulary, which may not fit in the memory (we need much more words to describe several domains, than one).      So the problem is: what is the best way to build topic models (or one topic model?) of various (and unknown) domains? And them apply this models to new (unseen) documents.


title: Problems with mixed model simulation

text: So i was just having a nice evening trying to simulate some hierarchical models by myself and estimate their parameters. First model I tried to simulate and estimate was $$y_{ij}=\beta_0+\varepsilon_{ij},$$ where $b_0=\gamma_{00}+u_{0j}$. So actually what I was estimating was $y_{ij}=\gamma_{00}+u_{0j}+\varepsilon_{ij}$. I have generated some data and estimated parameters $\gamma_{00}$ and variances of error terms and everything was fine:               set.seed(1)     N <- 100     nj <- 100     g00 <- 10     e <- rnorm(N*nj)     j <- c(sapply(1:N, function(x) rep(x, nj)))     uj <- c(sapply(1:N, function(x)rep(rnorm(1), nj)))     d <- data.frame(j, y=g00+uj+e)     library(nlme)     lme(y~1, data=d, random=~1|j)          Linear mixed-effects model fit by REML       Data: d        Log-restricted-likelihood: -14520.94       Fixed: y ~ 1      (Intercept)         10.00215           Random effects:      Formula: ~1 | j             (Intercept) Residual     StdDev:   0.7752422 1.012683          Number of Observations: 10000     Number of Groups: 100       Then I tried different model: $$y_{ij}=\beta_0+\beta_1 x_{ij}+\varepsilon_{ij},$$ where $\beta_0=\gamma_{00}+u_{0j}$ and $\beta_1=\gamma_{10}+u_{1j}$. So I had to estimate the equation $$y_{ij}=\gamma_{00}+u_{0j}+(\gamma_{10}+u_{1j}) \cdot x_{ij}+\varepsilon_{ij}=\gamma_{00}+\gamma_{10}x_ij+u_{0j}+u_{1j}x_{ij}+\varepsilon_{ij}.$$ I did the same thing that I had before, but this time `lme` did not converge. I tried this, but it didn't seem to work.               g10 <- 10     u0j <- uj     u1j <- c(sapply(1:N, function(x)rep(rnorm(1), nj)))     x1 <- rnorm(N*nj)     d1 <- data.frame(j, y=g00+u0j+(g10+u1j)*x1, x1)     lme(y~1+x1, data=d1, random=~1+x1|j)      Here is what last call to `lme` spit out:               Error in lme.formula(y ~ 1 + x1, data = d1, random = ~1 + x1 | j) :        nlminb problem, convergence error code = 1       message = false convergence (8)      What can you suggest me to do? Maybe the problem is with my model specification, redundant parameters, singular matrix or something else?


title: Finding events related to a specific feature in a noisy dataset

text: I am working on a tool that has to discover the most likely events related to the presence of a feature. I have a two dataset containing events. The first contains the observations presenting the feature I am interested into and the other contains the observation that should not be related to the feature. You can think at my data as a single set of events with attached a FEATURE/NON FEATURE label. As I mentioned in the comment, an event may appear multiple times in the dataset. Each event represents is discrete and represents a categorical value. Not all the events labeled with FEATURE actually relates to the feature. They are simply noise due to the fact that I cannot observe the feature in isolation. Is there any method to find which are the events that are most likely to be specifically related to the feature I am interested into? I was thinking at some sort of correlation analysis or at a Bayesian analysis but I am not really into statistics and I'd like to have your opinion.