title: Which distribution has the larger variability? Is variance absolute, or relative to the mean?

text: I apologize in advance, I'm new to statistics. I have a large (millions) dataset (the US Census American Community Survey) with 286 attributes. I've calculated the mean, variance and standard deviation for each attribute, and I would like to order them by "variability" (roughly, that is - if the concept I'm driving at is not a statistical one, I'll settle for one that is just intuitively appealing.) Obviously I can rank them by sigma, highest to lowest (ranging from very high to very low in absolute terms, on the order of thousands down to <1.0). But is that meaningful? Does a large sigma (or variance) mean it's more "variable", or should I be looking for the largest _RATIO_ of sigma to mean? I can't find any reference to a concept like that in my textbooks, but it seems to me to convey more meaning than a variance/sigma on its own. (The point, if you haven't guessed, is to reduce the dimensionality to something more manageable by discarding the attributes with the least variability.) Thanks very much, and cheers - Ed


title: Using taxonomic levels as factors in random forests: does it make sense? Is it needed?

text: I want to test the effect of a set of predictors (ecological and morphological factors) on a categorical response variable (an animal behaviour). As far as I've read, random forests do not make assumptions about data independence. Therefore, can I use species in my analysis disregarding their relatedness, that is, not using phylogenetic signals (which measure tendency of related species to resemble each other more than species drawn at random from the same tree)??? Is it correct? Otherwise, I could find for example, that one behaviour is well explained by body size, not because body size has some kind of relationship with that behaviour, rather because most species that show that trait are indeed of the same taxonomic family, and that family is mainly composed by same size species. If this is the case, can I use a taxonomic level as explanatory variables (e.g. genus, family or order)??? Would it give me the real importance of relatedness on that behaviour? (e.g.if in my importance list a get that my taxonomic level is the most important factor, can I interpret that the behaviour is mostly dependent on phylogenetic reasons rather than eco- morphological factors???) What about classification trees on this regard? I haven't found anything on papers and books in this regard. Thank you in advance for the help!


title: Is mixed measures ANOVA the correct test for my data?

text: I'm fairly new to statistics and I'm still trying to figure out the best way to analyse the data I have. The experiment has 2 groups of participants who perform 2 repetitions of a task that consists of 5 stages. All participants completed both repetitions for all stages, but one group had 8 participants while the other group only had 6. I have a about 100 dependent variables that I wish to examine, so my data looks a bit like this:               ID   Group    Repetition    Stage   DV1    DV2     ...     1    A        1             1       212.9  179.9   ...     1    A        2             1       144.8  134.7   ...     2    B        1             1       146.3  156.8   ...     2    B        2             1       128.6  178.2   ...      Group is a between-subjects factor while Repetition and Stage are within- subjects factors. I would like to determine whether Group and Repetition have a significant effect on each dependent variable within each stage (I am not interested in the effect of stage itself). I'm doing the analysis in R so I have the following code:               options(contrasts=c("contr.sum","contr.poly"))     mydata            = read.csv("data.csv",header=TRUE)     mydata$Group      = factor(mydata$Group)     mydata$Repetition = factor(mydata$Repetition)     mydata$Stage      = factor(mydata$Stage)     # for each stage     mydata = mydata[mydata$Stage==1,]     for (i in 5:(ncol(mydata)))      {        fit = aov(formula=as.formula(paste(names(mydata)[i],                                      "~ Group * Repetition + Error(ID/Repetition)")),                   data=mydata)     }      My questions are:   1. Is mixed measures ANOVA a valid test for this data? What's the correct way to test whether my data fits the assumptions of ANOVA in R? If this is not a reliable test, what's a possible alternative?   2. Have I defined the mixed measures ANOVA in R correctly? The various tutorials I've read define it in different ways so I'm a bit confused.


title: R Library for computing cepstral coefficients from Auto-Correlation coefficients

text: My main goal is to compute cepstral coefficients in R from an ACF vector (Auto-Coorelation function vector with discrete time-step). Correct me if I am wrong in terminology or whatever as I am not a statistician/signal-processor. My first question will be   * Are the elements of an ACF vector i.e. the auto-correlation-coefficients (ACC) linear-predictive-coefficients (LPC)? I think they are since an AR(p) process is a linear predictive model. If my above hypothesis is true, then I have found MatLab code for converting LPC coefficients to cepstral coefficients (CC): MatLab Code. Then my question would reduce to:   * Is there an R library which will convert LPCs to CCs? Thanks for all the help, and please correct my terminology if I am wrong!


title: Building linear model from exact variable measurements for use with noisy variable measurements

text: I want to build a linear model to predict a scalar output from a vector of noisy scalar variable measurements. I have two separate training data sets. One has output data and corresponding exact variable measurements. The other has exact variable measurements and corresponding noisy variable measurements. The noisy measurements of some variables are noisier (higher error variance) than others, and the noisy measurements of some variables are biased. I do _not_ have a single data set with output data and corresponding noisy variable measurements. How should I build my linear model? Should I use the exact variable measurements and ignore the fact that when the model is applied/used, noisy measurements of the variables will be input to the model? Or is there some way I can make use of what I can figure out about the noisy measurements of each variable when building the model? Can R help me with this problem?


title: Error in lmer - how to fix

text: I am asking another question (as each one is a new issue and I am trying to get these models run). However I don't want to become a nuisance either. However I ran the following and wonder whether anyone could see what I have done wrong. I know I am comparing a lot of factors but really want to find out which ones really contribute and would like to be able to figure out what is wrong so that I will know what to do in the future. The code I used is as follows: Note I was going to use MuMNi to run all the model combinations:               ABmodel <- read.csv(file.choose(), header=T)      I typed out `as.factor()` for `Brood`, `Site`, `Age.class`, `Tide`, `tide.h.l`, `tide.inc.out` and `MF.vs.OF` then attached lme4:               local({pkg <- select.list(sort(.packages(all.available = TRUE)),graphics=TRUE)            if(nchar(pkg)) library(pkg, character.only=TRUE)})      Then I attempted run the following model (I know some factors seem to be related but they are different):               fullABmodel1 <- lmer(Feeding ~ Site + Age.class + tide.h.l + tide.inc.out +                           Tide + HDp + NHDp + Foraging.Eff + MF.vs.OF + Age +                           HDr + NHDr + AllDr + (1|Brood), data=ABmodel, REML=FALSE)      But I got the following Error               Error in mer_finalize(ans) : Downdated X'X is not positive definite, 7.       I wish that I could attach my CSV file so that you could see it (although it contains additional columns for other analyses). Not sure what this means. Please note my understanding of the intricacies of models is limited. Does it help to know when I ran it as an `lme` in `nlme` I got a message saying               Error in MEEM(object, conLin, control$niterEM) :          Singularity in backsolve at level 0, block 1


title: Is there a way to predict change based on data?

text: > **Possible Duplicate:**   >  Is there a way to predict future change based on data? I do not know anything about statistics. I am a software engineer. I had a question in my mind and I think here is the place to get my answer. Suppose that I have data for 3 months. The data shows my visitors on a daily basis. Is there a possible way to know how many visitors I will have in e.g. 6 months? If yes, I guess there are some factors to be considered. What are they? Also, the number of visitors in the 180 day will be predicted only approximately and must have some risk on it. How can the risk be defined? Any links or info about this are very welcome. Thank you community.


title: Effect of two demographic IVs on survey answers (Likert scale)

text: I have 4 questions in a survey with Likert-style answers (i.e. 1–5) – and the participants have also answered the demographic questions (country and study- years) – which test do I need to find if the demographics, i.e. country and study-years, have each affected answers to each of the questions significantly? Many thanks.


title: question on plotting error values for treatment effects from multiple regression output

text: Say I have some data, where a dependent variable, `dv`, is a function of some independent variable, `iv`, and a categorical predictor, `cat`. Here are some example data below generated in R:               a<-c(1:100)     err<-rnorm(100,sd=30)     b<-a + err     c<-a + err + 20     cat1<-rep(0,100)     cat2<-rep(1,100)     iv<-c(a,a)     dv<-c(b,c)     cat<-c(cat1,cat2)     data<- data.frame(dv,iv,cat)      I then model `dv` as a function of `iv` and `cat` with this code:               summary(lm(dv~iv + cat, data=data))      and get the following output                               Estimate Std. Error t value Pr(>|t|)         (Intercept) -5.49606    5.09626  -1.078    0.282         iv           1.18250    0.07848  15.067  < 2e-16 ***     cat         20.00000    4.53086   4.414 1.67e-05 ***      Now, I want to plot the effect of cat using a standard bar graph- means and errors. So, based on the model, I calculate what the value of `dv` should be when `cat` is 0 and when `cat` is 1, using a common value of `iv` of 50. For my particular data set, I would get `dv` values of 53.62894 and 73.62894, for `cat` levels 0 and 1, respectively. My question is: which term from the model should I use for the error bars? Should I just use the std error estimate of the `cat` predictor? Or something more complex that integrates the error values of the intercept and `iv` parameter as well?


title: On tests for correlation in case of bivariate observations

text: We have bivariate observations {($X_i, Y_i$), i=1,2,...n}. We are interested in testing if they are uncorrelated. One test is based on $r$ (correlation coefficient): $r\sqrt{\frac{n-2}{1-r^2}}\sim t_{n-2}$. Another is based on fitting a regression model $Y_i = a+b X_i+e_i$ and testing for b=0. Are both the tests identical? Which one should be preferred?


title: What is the training score used for?

text: This is a pretty basic question. However, I read that the training score isn't useful in machine learning. From what I can tell, the training score only tells us if we have an overfit model. [E.g. if we have a training score of 100, something probably went wrong]. But are there any other reasons to need a training score?


title: Why does taking the Discrete Fourier Transform of a histogram make sense?

text: I think I have a serious misunderstanding of the Fast Fourier Transform. Currently, I'm using ROOT to try to analyze some data, and the Fast Fourier Transform (FFT) function applies only to histograms. Why does this make sense? How can one get an accurate transform knowing only the values which were measured and not exactly _when_ they were measured? (The following is clarification as per Glen_b's comment: (iii) - I have some data of voltage/time that I'd like to take the Fourier Transform of, which I thought would be a discrete Fourier transform, since the data is (obviously) discrete. I was under the impression that the Fast Fourier Transform is just a method of taking a discrete Fourier Transform, and the discrete Fourier Transform is directly analogous to the continuous one. (ii) & (i) - That being said, I only have limited experience with the continuous transform (from class). I now have some voltage/time data (the type of data you'd see on an oscilloscope) that I'm trying to transform to isolate some noise frequencies. ROOT, the platform through which I am performing the transformation, only has a Fast Fourier Transform method for histograms. Histograms are just counts of numbers, independent of the time. I'm confused as to how one can take the Fourier Transform of counts of numbers and get anything useful from it without knowing the time at which they occurred.)


title: Has the parameter $\frac{1}{\sqrt{2(n-1)}}$ a name?

text: Has the parameter $\displaystyle \frac{1}{\sqrt{2(n-1)}}$ a name? $n=\text{number of observations}$


title: Correct Evaluation of Random Forest on fixed training/test set

text: I have to test the performance of Random Forest on the same dataset (text classification) with about 118.000 instances of which about 1/3 is used for training and 2/3 is used for testing. The division is given, so I cannot modify either the test or the training set. I need to conduct an investigation about which parameters and/or features are best for the classification problem given a certain metric. Now, for the randomness of random forest each run I obtain a different output predictor, that will generate different results, and so the metric will vary from run to run. So my initial idea was to simply run the random forest a certain number n of times (say n=4) and take the average value of the given metric as representative of that particular configuration, in order to compare with the other possibilities. Probably a better approach would make use of a statistical test in order to establish to which degree we can assure that one model is better than another. How do you suggest to perform this considering the particular case I've mentioned above?


title: Modelling a probability distribution on different feature sets

text: I have a binary classification problem, and I use method A and method B to extract features, F1 and F2, for this problem from dataset X. Now, I train two models, y1 and y2, separately on the two extracted feature sets. Concretely:               p(y1 | F1) = ...     p(y2 | F2) = ...      Now I chain them together:               Y = p(y1 | F1) * p(y2 | F2)      And if the result of the above multiplication is greater than 0.5, I decide that the test example belongs to classification set `1`. Is this approach correct?


title: Bayes' factor vs. Bayes' Discriminant Rule

text: When we are comparing two models against some data, will we obtain the same (set of) posterior odds for the models both when we use the Bayes' factor and when we use the discriminant rule? If not, which one is considered more "accurate" and why is the other being used? I have found that the Bayes' factor is described also as: $B_{i,j} = \frac{\pi_{data}(D|M_i)}{\pi_{data}(D|M_j)} = \frac{L_{max}^{(i)}}{L_{max}^{(j)}} \frac{W_i}{W_j}$ where $W_i$ is the Ochham's factor: $W_i = \int_{V_{\theta_{i}}} \pi_{prior}(\theta_i|M_i) \frac{L_{(i)}(\theta_i)}{L_{max}^{(i)}} d\theta_j$ Is $\frac{L_{max}^{(i)}}{L_{max}^{(j)}} \frac{W_i}{W_j}$ anyhow relevant to the ratio of 'allocations' at each model ($\frac{allocations-to- model_i}{allocations-to-model_j}$), made by the Bayes' Discriminant rule? Edit-: After the response of conjugateprior, I would like to add this: I am wondering whether a discriminant analysis (such as ML discriminant rule i.e. Bayesian discr. for priors 1/2 and 1/2) is an acceptable way of doing model selection when we are comparing two models from which we inferred posterior information about some parameter/hyperparameter via Bayes' theorem. I.e. although we assume that the data are generated by one model, if you do discriminant analysis, won't the percentage of the mixture reveal the model that the Bayesians' describe as 'more likely'?


title: Simple Question regarding t-tests

text: I have the mean from a previous study (the Brief Symptom Inventory) that I would like to compare the mean found in my study. Both samples have a population over 100. I do not have the entire data set for the previous study. I typically use SPSS for my stats. I was hoping to use a paired sample t-test as I am comparing the mean scores of the same scale. However, I was told by the program this was not possible with the values I entered. Is there another way that I may do this? Basically, looking to compare the means to determine if there is a statistical difference without having the entire population.


title: Which design is most appropriate?

text: I have data from an experiment testing a new treatment. Each subject has 4 data points, 2 of the new treatment and 2 of the control. It is like a paired analysis, but instead having 1 observation of treatment and 1 of control, I have 2 of each. Which design is it (complete block, ....). How should I analyze this data ? A mixed model (response is continuous) ? How bad will it be if I average the 2 of each treatment per each subject ?


title: Determine Number of Clusters Using Silhouette Method and Distance Matrix

text: I have a tree and its associated distance matrix. Now I want to cut the tree to give a desired number of clusters. Luckily there are a number of methods for determining the number of clusters as discussed many times on Cross Validated. I want to use the silhouette method to determine the number of clusters since I only have the distance matrix. But after doing some reading on this method, it seems that the method requires the number of clusters as input. And I want to get the number of clusters as output. **How do I use the silhouette method together with the distance matrix to determine the number of clusters?** I could guess you would maybe have to iterate the number of clusters and score with the silhouette method.


title: Interpreting output from the "dr" package in R

text: I'm dealing with a "large p small n" modelling problem and have decided to use the dr package to, as I understand it does, remove unnecessary / confounding parameters from my model, rather than do this manually. However I'm having a bit of trouble discerning as to what I am supposed to do, and would be very grateful for any help or guidance. Here is what I have done: I first used random forests in rattle to arrive at list of the 134 most important variables Then put those 134 variables into a formula object in r i.e               f1<-formula(dep_var ~ 0 + b1 +b2 + .... b134)      I have used the `dr` function to create the necessary inputs to `dr.compute`               reduction<-dr(f1,data)               computed.reduction<-dr.compute(x=reduction$x                                        ,y=reduction$y                                    ,method="sir"                                    ,weights=reduction$weights                                        ,M=reduction$M                                    ,numdir=5)      I have then passed the computed.reduction object to the dr.step function to get it to remove variables that are unnecessary               stepped.reduced<-dr.step(computed.reduction                            ,update=TRUE                            ,minsize=20                            )      However this is not doing anything, returning the message               No more variables to remove      I cannot find anything contained within the `stepped.reduced` object that would correspond to a shortened list of variables or anything explicit enough in the documentation. Any help greatly appreciated I would also be very grateful for any information about other packages/methods I can achieve automated removal of variables from a formula.


title: Is it possible to use the Breusch-Pagan Lagrange multiplier test (xttest0) in Stata for unbalanced data?

text: Is it possible to use xttest0 in Stata with unbalanced panel data? I want to test whether the I should use pooled OLS or random effects estimation. What does this test actually do?


title: How to test group difference between nominal and ratio data with nonparametric tests?

text: In my data the dependent variable ‘vocabulary uptake’ is measured as scores ranging from 0-10 and the independent variable is yielded from a ‘reading intervention’ with 4 different reading conditions (glossed, no glosses, etc.). There are about 80 participants in each of the four groups. The data is not normally distributed (I tried to transform the data). Given the fact that there is a natural limit (0), this is probably not surprising. I want to find out what the best way is to determine whether these groups are different or not. My questions:   1. I am sure that the independent is a ‘nominal’ variable, but I am in doubt about the dependent – I cannot decide whether to treat it as ‘ratio’ data or as ‘ordinal’. The latter would make sense, as the scores can be ranked, but ‘discrete ratio’ makes sense as the values are measured, have numerical value, the difference between, e.g. 1 and 2 is the same as between 3 and 4, and a score of 0 means that there is none of that variable. ‘Ordinal’ data would be easy, but if I categorize the data as ‘discrete ratio’, my problem is that – to my knowledge – the standard nonparametric tests only work with continuous data. What can I do?   2. What test + post-hoc test could I use to find out if there are significant differences between the groups? When working with a nominal dep. and a discrete ratio data indep. variable I would violate the assumptions for most nonparametric standard tests, which assume dichotomous independents and continuous/ordinal dep. variables (Mann-Whitney U, Wilcoxon) or categorical data for the independent variable and continuous/ordinal data for the independent variable (Kruskal-Wallis, Friedman’s test). Also there is a problem with tied values.  I am working with SPSS. Thanks for any suggestions! Your help would be GREATLY appreciated!


title: Questions about $R^2$, VIFs and very non normal input variables

text: I have been working with a small part of my dataset trying to eliminate variables and do some micro models. When analysing my micro set I initially found a few high correlations with inputs (0.95 etc) which so I simply removed the variables I then, fed all the variables left of interest into an OLS to which I got some quite high VIFs and standard errors in my estimates. I then, used a bit of PCA and my inflated VIFs were gone.   I didn't use PCA to reform the variables into factors as I wanted to maintain my variables. The $R^2$ of my model with no PCA or VIF consideration has a 0.15 better fit than my adjusted model but with the high standard errors. How do I determine which model is 'most right'? i.e. - Do i need to be concerned with VIFs? I have read the article here: http://www.statisticalhorizons.com/multicollinearity Another side question. Many of my original inputs were extremely heavily skewed due to top out real values. In fact this is where the peak lies in the distribution. Is it right to assume this isnt a problem as:               - my output variable is pretty much normal     - my model residuals by predicted are normal      i.e. - Do you only need to worry about input transformation when residuals aren't normal?


title: R Forward and backward Selection

text: I have a data set with large number of attributes some are not relevant and some are relevant for the regression model. My approach was to do forward and backward selection to identify a starting point, such as which attributes I should drop from the analysis. I did a forward and backward without any log transformation for the attributes and issue is that the best model provided by forward selection and best model provided by backward selection are different. I am confuse on what I should pick to go forward with. In addition when you have large number of attributes what is the best way to identify predictors for regression model.


title: Distribution of  total number of car accidents in a lifetime of the car

text: the number of car accidents every year has Poisson distribution with mean $\lambda$ and lifetime of the car has exponential distribution with mean $\frac{1}{\mu}$. how can find distribution total number of car accidents in a lifetime of the car?


title: How can I create a scatterplot in R using the plot function to control for covariates?

text: I would like to know if it is possible to create a scatterplot while controlling for covariates, such as in partial correlation. I am using R software and my code is below for the basic scatterplot. I am not interested in multiple lines of best fit or multiple scatters per graph. I am also not interested in creating a lattice of scatterplots with all of the variables. I can't seem to find any code that will allow me to parse the effect of the covariates from my x-y scatterplot. All three covariates are continuous. Any ideas would be helpful. Thank you               plot(pmc$reject, pmcp$LPA, main="r(Parenting, Left Amygdala)",         xlab="Parenting Age 2", ylab="Amygdala Reactivity Age 15", pch=19)     abline(lm(pmc$reject ~ pmcp$LPA), col="red") # regression line (y~x)


title: How to compare powers of two tests based on their p-values?

text: Suppose we are testing null $H$ versus $K$. Here $H$ and $K$ are viewed as two disjoint sets of distributions of sample $X$ Given a sample $X$, one test rule rejects null if and only if $T_1(X) \geq c_1$, and the other $T_2(X) \geq c_2$. I was wondering how to compare the powers of the two tests at a sample distribution $F \in K$, based on their p-values? One way I saw from an example (without explanation) is to first fix a level of significance $\alpha$, and then calculate the probability of each test's p-value no greater than $\alpha$, when the sample distribution is $F \in K$, i.e., $$ \mathrm P_{X \sim F} \{[ \sup_{F' \in H} \mathrm P_{X' \sim F'} (T(X') \geq T(X))] \leq \alpha \} $$ The bigger the value is, the more power the test has at $F \in K$, but I was wondering why this is true? Does it depend on $\alpha$?


title: Interpretation of the OLS regression coefficient of a proportion dependent and interaction independent variable with dummy

text: I have two regression models; I am not sure how to interpret the coefficients. In the first one the dependent variable is earnings management proxied by (abnormal accruals scaled by total assets), the independent variable is an interaction between a dummy variable and a ratio variable M/B : $A: \text{earnings mgt}= -0.057MB + -0.10OC + 0.057MB*OC$ where OC is a dummy variable. Should I say OC managers are likely to manage earnings by 5.69 cents per dollar in MB more than other managers. Or OC mangers are likely to manage earnings by 0.0569 percent for one percent increased in MB more than other managers The second regression has interaction between two dummy variables instead. MB is presented as dummy =1 if higher than industry average $B: \text{earnings mgt} = -0.056 MB + -0.004 OC + 0.078 MB*OC$ How can I interpret the coefficients in this case?


title: Statistical test to compare different treatments over time

text: In my experiment I have 4 treatments (the control plus the variables: A, B, C) and I am comparing the growth among them. I have 4 replicates for each treatment. Thus 16 samples in total. Every other day, a cell count was made. In the total 11 cells counts for each treatment was made. In the end, the mean for the replicates in each day were calculated and plotted ![Algae Density](http://i.stack.imgur.com/0cf9S.jpg). I have 11 density values (means) for the control + 11 density values for the A + 11 density values for the B + 11 density values for the C. Some data (number of cells in every other day of experiment, i.e. day 0, day 2, day 4,...):               Control Means: 111 / 263 / 370 / 460 / 627 / 715 / 802          100x Means:__  101 / 226 / 301 / 408 / 495 / 613 / 668          200x Means:__ 103 / 245 / 296 / 402 / 525 / 614 / 570          400x Means:__ 113 / 229 / 317 / 407 / 613 / 556 / 573      The growth rate was calculated and plotted as well ![Growth Rate](http://i.stack.imgur.com/7QfVF.jpg) Some data (growth rate means, µ , between days, i.e. µ mean between day 0 and 2, µ mean between day 2 and 4, µ mean between day 4 and 6,...):               Control Means: 0.42 / 0.17 / 0.10 / 0.15          100x Means:__ 0.39 / 0.14 / 0.15          200x Means:__ 0.42 / 0.09 / 0.14          400x Means:__ 0.35 / 0.16 / 0.12      I'd like to see if there is any difference among these treatments or if the algae has the same growth independent of the treatment I use. (p.s.: in the treatments I used waste water diluted in the same seawater I used to prepare the controls using dilutions of 100x, 200x and 400x to grow the algae)


title: Mutual information/pointwise mutual information for measuring prediction

text: I want to measure how well I predict a vector $Y$ (vector not a label) for observation $X$. Both $X$ and $Y$ have the same set of features ($1\times n$). For that, I thought of "scoring" the correlation between each pair of the $n$ features, based on previous data I have, i.e. create $n\times n$ matrix $S$, where $S(i,j)$ is the score how well feature $i$ goes with feature $j$. When predicting for observation $X$, a vector $Y$, I want to take the dominant features (say, above threshold $e$) from $X$ and $Y$, and look in $S$ for their score, if it is above threshold $t$ then I can say my prediction was ok. My questions are:   1. Is there any other way to do that?    2. What can I use as the scoring function? I thought of mutual information, or the frequencies of $i$ and $j$ appearing together in the data set.


title: Difference in difference model with 2 treatment groups

text: I have an experiment which has 2 treatment groups(effects) and a control group. Up until this point in my analysis, I've been carrying out a DID using a regression equation of the form: Y= γD_t + β1(TREAT1) + β2(TREAT2〗+ τ1(TREAT1*D) + τ2(〖TREAT2*D) + ε Basically I've been including both treatment effects in the same regression equation which I've been running on STATA. My question is should I be including both interaction terms in one equation, or should I be splitting my data and analysis so that I'm estimating the ATE for both groups separately? Thanks very much Rob


title: Comparing two paired data sets

text: I am trying to see compare tow paired data sets which do not have normal distributions, and I am not sure how to do that. Can I use t-test? I would appreciate your help very much. This what they look like: ![enter image description here](http://i.stack.imgur.com/3ayaK.png)


title: Empirical change of DAG for Bayesian Network

text: While working on some problem in bioinformatics I applied bayesian network algorithm for classification purposes. As predictors I took a window of sequence of aminoacids and for dependent variable I used some feature of interest of a central acid. For DAG estimation I used several different algorithms which gave me identical network structures, and all of them used BIC score during the DAG search. Obtained DAG structure consisted of two chains: X7-X5-X3-X1, and X6-X4-X2, while all Xs where connected with Y. Using this network I got my classification accuracy up to 65%. The problem is that motivation for my small research was to check if Bayesian Networks can beat Markov chain procedure developed by my collegues for same problem. Basically my collegues obtained 81% accuracy by applying simple Bayesian procedure for model in which sequence of amino-acids was a Markov chain (they also used some pretty tricky empirical methods to obtain such a great result). So, because my collegues showed that Markov Chain of second order gives best classification accuracy, I just modified my Bayesian Network by connecting two chains described above, so that each amino-acid was connected with two others (not one as before). I then used obtained structure and got classification accuracy 92% (which is pretty awesome for my task). So, my questions are:   1. How can I justify changes of DAG I made to improve my model. Of course I can (and I will) mention results of my collegues, but still, its more intuitive explanation. Maybe there is better reasoning then intuition here   2. What can be a possible reason for algorithms not finding optimal structure? (My current explanation is that BIC has just too strong penalty for model's complexity and algorithms just got into some local extremum, and weren't able to get out of it because of BIC's complexity penalty. DO you think it's worth to try and develop some more accurate score for my particular problem, to obtain correct DAG structure? What may be the other explanations?)


title: What is a Classifier?

text: I cannot find the general definition of what is a classifier? I understand how it can work, but I can't come to a definition.


title: arules R - How to ignore unknown item label in appearance list?

text: I am mining for association rules using arules. I often run new transaction sets with the same code. I dont change the list of restrictions on which items may appear in the rhs and lhs, but sometimes the transaction set does not have an item that is in my rhs or lhs list. Is there any way to make R ignore this?               > library('arules');     > txn = read.transactions(file="Query1.csv",rm.duplicates=FALSE,format="single",sep=",",cols=c(1,2));     > llist <-scan('gpis.txt', what="", sep="\n");     Read 5505 items     > rlist <-scan('hccs.txt', what="", sep="\n");     Read 127 items     > rules <- apriori(txn, parameter = list(confidence = 0.5), appearance = list(rhs = rlist, lhs = llist));     Error in asMethod(object) :        "blah" is an unknown item label, "blab" is an unknown item label.....


title: Based on groups A and B, how to classify new set of observations into the groups using regression?

text: Let's say there's two product types that already established, Search goods and Experience goods(Search goods refer to goods whose quality is easily evaluated before purchase. And Experience goods refer to goods whose quality can only be evaluated(experienced) after purchase. for more detail, please refer to here.) Based on the previous survey, I have a list of goods classified as Search goods or Experience goods. Between the two types of products, as many researchers pointed, there are different features of online reviews that cover those two types of products.(for example, online reviews for search goods are shorter than reviews for experience goods.) Also, I have large amount of reviews for products of the two types. Now, I want to classify new products into the product types, based on the features of reviews of those products. In this case, I think using regression with dummy coding will be nice. dummy coding will be:                                   Type1    Type2     Experience goods      1        0     Search goods          0        1     New product           0        0      Based on this coding, I can check the significance of Type1 and Type2 variable.(I guess, If Type1 is significant and Type2 is insignificant, the new product can be classified as Search goods, vice versa. Or, If both of the Types are insignificant, the new product cannot be classified in that grouping scheme.) Main point is, I'm not sure that it is persuasive to conduct this modeling in order to classify the new product. If it isn't please let me know better idea for classification. (But regression will be highly preferred.)


title: How to test for a home bias in a data set?

text: I'm trying to think how I can test a data set I have, for transactions between buyers and sellers, for a home bias. My hypothesis is that buyers would rather buy from a merchant from their home country (and that merchants don't necessarily care to whom they sell, as long as they make a sale). So, let's suppose there are: $m$ sellers, $n$ buyers (let's assume $m < n$ for convenience), and $x$ different countries. Each buyer buys one unit of the product, which is exactly the same for each seller. Also let's assume that a buyer with a home bias would rather buy from a seller from his own country, and if there is no home bias the buyer is indifferent between the sellers (uniform distribution). I thought about doing a $\chi^2\mathrm{-test}$, but the thing is that if there are many countries, and we assume that the difference would be for just one country, than the $\mathrm{d.f.}$ would rise and so the critical value, and we won't be able to find an effect. For example, there are a total of $61$ countries, we'd expect $100$ purchases from each. One country has $160$ purchases and the other $60$ have $99$, the $\mathrm{d.f.}$ is $60$, and the critical value at the $95\%$ confidence level is ~$79$. The $\chi^2$ value is $36.6$, so we won't find a statistically significant difference. I hope I'm clear about my problem. Would love to hear your thoughts. Thank you


title: Quadratic form of a bivariate normal

text: This is a homework problem. Let $(X,Y)\sim N(\mu_1,\mu_2,\sigma^2_1,\sigma^2_2,\rho)$. Show that if $\sigma_1,\sigma_2 >0,|\rho|<1$, then $$ \dfrac{1}{1-\rho^2}\left\{\dfrac{(X-\mu_1)^2}{\sigma^2_1}-2\rho\dfrac{(X-\mu_1)(Y-\mu_2)}{\sigma_1\sigma_2}+\dfrac{(Y-\mu_2)^2}{\sigma^2_2}\right\}$$ has a $\chi^2_2$ distribution. EDIT: There is a generalization of this: $(X-\mu)^T\Sigma^{-1}(X-\mu) \sim \chi^2_n$ where $X=(X_1,\ldots,X_n)^T$. The proof is theorem 7 here (http://www2.econ.iastate.edu/classes/econ671/Hallam/documents/QUAD_NORM.pdf) Can anyone think of a better way to show this?


title: Jarque-Bera normality test in R

text: Jarque-Bera normality test has significant p-values even when there is skewness and kurtosis. Does that mean test is infering data distribution is approximately normal?


title: Machine learning algorithms for context-based relevancy scoring

text: I'm building a Getting Things Done app. My principle complaint with most GTD / TODO apps is they show me too much. When I'm at work, I don't want to see tasks I'm planning to do on the weekend. When I'm at home, I don't want to see work tasks. When out running errands, I just want to see my errands tasks, etc. I've been thinking about this problem for awhile and with some help from Bret Victor's article "Magic Ink" I've designed a system I think will work quite well. At the heart will be a system that learns what projects/tasks/tags I want to see based on the context I'm in. So when I'm at work, the system should know to only show me work items. When I'm at home, work items should be hidden. What I'm trying to figure out right now is which machine learning algorithm would be "best" for my system. How the system will work is at each page load and every 10 minutes or so while the application is open, a context object will be created on the client based on the current time, day of week, location, and anything else I can think up. Based on the context, every project/task/tag will have a relevancy score calculated. The most relevant projects/tasks/tags will be shown on the home page of the app. Irrelevant items would be hid completely. The system will learn based on what the user does after opening the app. If she clicks on a project on the front page, that will strengthen the connection between the current context and that project and any tags on that project. If she has to change one of the global filters (say at 9pm she turns on the work tag and turns off all other tags as she has a late meeting with her team in India), the system should quickly learn that at around 9pm, items tagged with "work" are really relevant. It should also learn over time more subtle things like this meeting happens only on Mondays and Thursdays but not on the third week of the month and also never when she's not at home (e.g. she's traveling so the meeting is canceled). I've read a fair bit about machine learning and have some practical experience but I'd like some expert guidance about what would be the appropriate algorithm for this use case. My best guess right now is to use some variant of a neural network but I might be completely off-base. I'd prefer simpler algorithms that I can implement myself as I'd like to learn more about ML. Also that the algorithm is not a black box that I can't tweak.


title: About variance mixture models and probability distributions

text: I was wondering if anyone knows a good resource to learn about variance mixture models ? My interest is in particular the normal variance mean mixture. I know what they mean with their definition of $Y=\mu+\beta X+\sigma\sqrt{X}Z$, with $Y$ the variance mixed distribution, X the mixing distribution and Z the standard normal distribution, which are independent of each other. In an article on ArXiV the immediately suggest the generalized hyperbolic distribution as a solution. Now I was wondering how one would go from the prescription above for $Y$ to a formula for the distribution. Is it by guessing and hoping it works, or is there some kind of procedure ?


title: Time dependent coefficients in R - how to do it?

text: **Update** : Sorry for another update but I've found some possible solutions with fractional polynomials and the competing risk-package that I need some help with. * * * ## The problem I can't find an easy way to do a time dependent coefficient analysis is in R. I want to be able to take my variables coefficient and do it into a time dependent coefficient (not variable) and then plot the variation against time: $\beta_{my\_variable}=\beta_0+\beta_1*t+\beta_2*t^2...$ ## Possible solutions ### 1) Splitting the dataset I've looked at this example (Se part 2 of the lab session) but the creation of a separate dataset seems complicated, computationally costly and not very intuitive... ### 2) Reduced Rank models - The coxvc package The coxvc package provides an elegant way of dealing with the problem - here's a manual. The problem is that the author is no longer developing the package (last version is since 05/23/2007), after some e-mail conversation I've gotten the package to work but one run took 5 hours on my dataset (140 000 entries) and gives extreme estimates at the end of the period. You can find a slightly updated package here - I've mostly just updated the plot function. It might be just a question of tweaking but since the software doesn't easily provide confidence intervals and the process is so time consuming I'm looking right now at other solutions. ### 3) The timereg package The impressive timereg package also addresses the problem but I'm not certain of how to use it and it doesn't give me a smooth plot. ### 4) Fractional Polynomial Time (FPT) model I found Anika Buchholz' excellent dissertation on "Assessment of time–varying long–term effects of therapies and prognostic factors" that does an excellent job covering different models. She concludes that Sauerbrei et al's proposed FPT seems to be the most appropriate for time-dependent coefficients: > FPT is very good at detecting time-varying effects, while the Reduced Rank > approach results in far too complex models, as it does not include selection > of time-varying effects. The research seems very complete but it's slightly out of reach for me. I'm also a little wondering since she happens to work with Sauerbrei. It seems sound though and I guess the analysis could be done with the mfp package but I'm not sure how. ### 5) The cmprsk package I've been thinking of doing my competing risk analysis but the calculations have been to time-consuming so I switched to the regular cox regression. The crr has thoug an option for time dependent covariates:               ....     cov2        matrix of covariates that will be multiplied                  by functions of time; if used, often these                  covariates would also appear in cov1 to give                  a prop hazards effect plus a time interaction     ....      There is the quadratic example but I'm don't quite follow where the time actually appears and I'm not sure of how to display it. I've also looked at the test.R file but the example there is basically the same... ## My example code Here's an example that I use to test the different possibilities               library("survival")     library("timereg")     data(sTRACE)          # Basic cox regression         surv <- with(sTRACE, Surv(time/365,status==9))     fit1 <- coxph(surv~age+sex+diabetes+chf+vf, data=sTRACE)     check <- cox.zph(fit1)     print(check)     plot(check, resid=F)     # vf seems to be the most time varying          ######################################     # Do the analysis with the code from #     # the example that I've found        #     ######################################          # Split the dataset according to the splitSurv() from prof. Wesley O. Johnson     # http://anson.ucdavis.edu/~johnson/st222/lab8/splitSurv.ssc     new_split_dataset = splitSuv(sTRACE$time/365, sTRACE$status==9, sTRACE[, grep("(age|sex|diabetes|chf|vf)", names(sTRACE))])          surv2 <- with(new_split_dataset, Surv(start, stop, event))     fit2 <- coxph(surv2~age+sex+diabetes+chf+I(pspline(stop)*vf), data=new_split_dataset)     print(fit2)          ######################################     # Do the analysis by just straifying #     ######################################     fit3 <- coxph(surv~age+sex+diabetes+chf+strata(vf), data=sTRACE)     print(fit3)          # High computational cost!     # The price for 259 events     sum((sTRACE$status==9)*1)     # ~240 times larger dataset!     NROW(new_split_dataset)/NROW(sTRACE)          ########################################     # Do the analysis with the coxvc and   #     # the timecox from the timereg library #     ########################################     Ft_1 <- cbind(rep(1,nrow(sTRACE)),bs(sTRACE$time/365,df=3))     fit_coxvc1 <- coxvc(surv~vf+sex, Ft_1, rank=2, data=sTRACE)          fit_coxvc2 <- coxvc(surv~vf+sex, Ft_1, rank=1, data=sTRACE)          Ft_3 <- cbind(rep(1,nrow(sTRACE)),bs(sTRACE$time/365,df=5))     fit_coxvc3 <- coxvc(surv~vf+sex, Ft_3, rank=2, data=sTRACE)          layout(matrix(1:3, ncol=1))     my_plotcoxvc <- function(fit, fun="effects"){         plotcoxvc(fit,fun=fun,xlab='time in years', ylim=c(-1,1), legend_x=.010)         abline(0,0, lty=2, col=rgb(.5,.5,.5,.5))         title(paste("B-spline =", NCOL(fit$Ftime)-1, "df and rank =", fit$rank))     }     my_plotcoxvc(fit_coxvc1)     my_plotcoxvc(fit_coxvc2)     my_plotcoxvc(fit_coxvc3)          # Next group     my_plotcoxvc(fit_coxvc1)          fit_timecox1<-timecox(surv~sex + vf, data=sTRACE)     plot(fit_timecox1, xlab="time in years", specific.comps=c(2,3))      The code results in these graphs: Comparison of different settings for coxvc and of the coxvc and the timecox plots. I guess the results are ok but I don't think I'll be able to explain the timecox graph - it seems to complex... ## My (current) questions   * How do I do the FPT analysis in R?   * How do I use the time covariate in cmprsk?   * How do I plot the result (preferably with confidence intervals)?


title: Computational intelligence-based methods for time series forecasting?

text: I know that ARIMA is one approach widely used in forecasting of univariate time series. What are some _computational intelligence_ based methods that are reported to be successfully used in the same area? If possible, please recommend some articles/reports presenting such approaches. Thanks in advance.


title: Tuning the birthday paradox

text: There are _m_ disjoint sets _X_ _i_ and a function _f_ : [1.. _n_ ] → [1.. _m_ ] such that _Y_ _j_ = _X_ _f_ ( _j_ ). I am given oracles / samples from random variables _x_ _i_ and _y_ _j_ which uniformly select elements from _X_ _i_ and _Y_ _j_ respectively, but each sample from _x_ _i_ and _y_ _j_ has a cost (independent of _i_ or _j_ ). I would like to compute _f_ by sampling from _x_ and _y_ and looking for collisions, since _X_ _i_ and _Y_ _j_ are either equal or disjoint, but I'd like to keep the cost down, so using the fewest samples from _x_ and _y_. If _n_ and _m_ were approximately equal, I think I should just sample evenly from _x_ and _y_ , and the birthday paradox should allow me to identify at least one value of _f_ ( _j_ ) fairly quickly (taking about the square root of the size of the largest set in _Y_ ), and I think in general all of _f_ fairly quickly. However, in my situation _m_ is large compared to _n_ (say _m_ = 50000 different _X_ s, _n_ = 50 different _Y_ s). Each _X_ _i_ has about 106 elements. It seems in this case one should "focus" on _y_ more than _x_ since most _X_ _i_ will end up being irrelevant (but still expensive), but on the other hand, only sampling from _x_ once per _i_ seems to completely negate the birthday paradox, requiring an expected number of | _Y_ _j_ | samples from _y_ _j_ to find _f_ ( _j_ ). > How often should I sample from _y_ compared to _x_ in order to minimize the > total (expected) number of samples required to compute _f_? In other words, express _p_ as a function of _m_ , _n_ , _k_ = max(| _X_ _i_ |) where 0 ≤ _p_ ≤ 1 is the minimizing argument of the function which gives the expected number of samples of the process which samples uniformly from the _x_ _i_ with probability _p_ , and uniformly from the _y_ _j_ with probability 1− _p_ , and stops as soon as _f_ has been computed (that is, as soon as for every _j_ , one has sampled a value of _y_ _j_ equal to a value sampled from an _x_ _i_ ).


title: Minimum items for questionnaire

text: Is there a minimum item for questionnaire? Some researchers said that the total items for teenagers ranging from 80-110 items


title: Does the number of IV after stepwise model selection depend on the amount of data?

text: I have 3 different DV that I try to model with 3 distinct models (linear mixed models) using the same set of IV. I found that the DV that I have the least amount of data for also has the lowest number of IV's in the final model, while the DV that has the largest amount of data has a larger number of remaining IV (including interactions). I tend to believe that this might be an artifact of data-availability. Might that be the case? And if so: Is that a serious problem when reporting these results in peer reviewed journal? I know that many argue against stepwise reduction procedures, however my feeling is that (at least in my case) the selection mostly comes to the same final models irrespective of minor changes in the inital model. Are the constraints against model selection mostly valid in case of large numbers of IV or is it even in case of a relatively small number of initial IVs that there might be completely different results?


title: Using R's aov() for two-way repeated-measurement ANOVA

text: We analyze blood vessel samples of patients in an organ bath which gives us force values. Each subject (patient) donates two types of blood vessels (pulmonary vein and pulmonary artery). The organ bath measurements include dose-response curves with a pharmaceutical, i.e. different concentrations are applied. From what I understand, this should be analyzed using a two-way RM ANOVA test. There are two variables (vessel type and concentration), and all levels of both are applied to each subject. The data table looks like this:               subject  vessel  conc  force     20110818 PA      -12   0     20110818 PV      -12   0     20110818 PA      -11   0.09     20110818 PV      -11   0.15     ...      There are no missing values or anything like that to consider. Perusing several R tutorials, I came up with the following query:               aov.ex2=aov(force~(vessel*conc)+Error(subject/(vessel*conc)),data=data.ex2)      However, the results are somewhat unclear to me:               > summary(aov.ex2)          Error: subject               Df Sum Sq Mean Sq F value Pr(>F)     Residuals  1 15.370  15.370                         Error: subject:vessel            Df Sum Sq Mean Sq     vessel  1 166.72  166.72          Error: subject:concentration                   Df Sum Sq Mean Sq     concentration  1 3134.3  3134.3          Error: subject:vessel:concentration                          Df Sum Sq Mean Sq     vessel:concentration  1 148.32  148.32          Error: Within                Df Sum Sq Mean Sq F value Pr(>F)     Residuals 219   3031  13.840                     This is quite a lot of information, but in contrast to the nice examples on the web I can't seem to find the useful part of the information: p values and whether or not there are significant differences. Is there something wrong with my input, or does the result simply mean there is nothing interesting to see in these data?


title: Grouped income data in Stata or SAS

text: I have an income distribution dataset which I am attempting to manipulate in Stata (or SAS). I need to come up with a table listing the number of households earning between one income level and another. The income levels are regularly spaced, say at $10,000 intervals for example. In other words, for $i=30,000$, I want to know the number of households $H_i$ receiving between 20,000 and 30,000. And I want to know this for each i, at 10,000 dollar intervals, from 10,000 up to oh let's say 280,000. $H_i=\int_{i-10000}^{i}f(x)\: dx$ $(0 < i\leq280000)$ (Where f(x) is the underlying density function representing households plotted along the y-axis, and x is income level, plotted along the x-axis. The distribution will have the skewed (lognormal-like) shape typical of income distributions.) Does anyone know the Stata (or SAS) command that would give me this info? In Stata I can get the income at each percentile of the population using the "centile" command, but this doesn't quite give me what I want. Thanks


title: Evaluating and combining methods based on ROC and PR curves

text: I am evaluating and combining a few binary classification models. I am using the ROC and PR curves to evaluate their performance. The problem I am having is that as I try to improve the method, I am improving the AUC-ROC but the PR curve suffers. For example: ![ROC curve](http://i.stack.imgur.com/yH47K.png) ![PR curve](http://i.stack.imgur.com/4MuJj.png) As an aside, I am actually adding a weak learner to Method 1 to arrive at Method 2, and then adding another weak learner or two to arrive at Method 3. When I was only evaluating AUC-ROC, it looked fine, but when I saw the PR curve, it seems I have been degrading the performance. Now it seems that the weak learners are doing better at points lower in the ranked list. But this is only for one dataset training/test split. What would be a principled way to investigate what is going on and come up with a way to use the weak learners so as to improve both ROC and PR curves? **Update:** To visualize this, I am showing the weak learner that I am adding to Model 1 to arrive at Model 2 here: ![weak ROC](http://i.stack.imgur.com/tVUPu.png) ![weak PR](http://i.stack.imgur.com/7WulV.png)


title: Likelihood of a document-term matrix

text: I'm trying to deeply understand the Probabilistic Latent Semantic Indexing (PLSI) algorithm. The generative model underlying the algorithm is the folowing:   1. Select a document $d \in D$ with a probability $P(d)$   2. Select a topic $t \in T$ with probability $P(t|d)$   3. Generate a word $w \in W$ with a probability $P(w|t)$ where the sets $D$, $T$ and $W$ are finite and fixed and all the probability distributions are multinomial. Given an observation, which is a corpus $C$ of documents, the author computes the likelihood function as: $$L = \prod_{d\in D}\prod_{w \in W} n(d,w) p(d)\sum_{t \in T} p(w|t)p(t|d)$$ where $n(d,w)$ is the number of occurrences of the term $w$ within the document $d$. My question is: given that we are not considering the order of the words within a document, shouldn't we consider, for each document, also the _multinomial coefficient_ ? $$ C_d = \binom{\sum_{w}{n(d,w)}}{n(d,w_1), n(d,w_2), \dots , nd(w_W)}$$ Could you explain me why it has been omitted?


title: Classifier for only one class

text: In a simple classification, we have two classes: class-0 and class-1. In some data I only have values for class-1, so none for class-0. Now I am thinking about making a model to model the data for class-1. So, when new data come, this model is applied to the new data and finds a probability saying how probable that new data fit this model. Then comparing with a threshold, I can filter inappropriate data. My questions are:   * Is this a good way to work with such problems?   * Can a RandomForest classifier be used for this case? Do I need to add artificial data for class-0 which I hope the classifier regards as noise?    * Any other idea may help for this problem?


title: Confusion regarding dynamic programming in CRF paper

text: I was referring to this paper related to CRF. However, I didn't get the part related to dynamic programming. Further I have question related to feature function selection. In this paper they have given this feature function $$ F_j({{\bf{y}}},{{\bf{x}}}) = \sum_{i=1}^{n}f_j(y_{i-1},y_i,{\bf{x}},i) $$ I didn't get why we have to sum from $i=1:n$. Is it because it allows to have particular $y_i$ and $y_{i-1}$ values to be satisfied by any ordered pairs like $(y_2, y_1),(y_3,y_2)$ and so on. Now suppose I have a feature function such that I have to have this value $A, B$ at particular positions like $y=1$ and $y=2$, then I cannot sum from $i=1:n$ right. Or even if I do it will be all zeros right? Now my question related to dynamic programming. I didn't understand the expression given after the following in the paper it must be possible to eﬃciently compute the expectation of each feature function with respect to the CRF model distribution for every observation sequence $x_{(k)}$ in the training data given by... So can anyone please explain with some short examples. I mean considering training data with 1-2 short sentences and with few labels. I think it will help me to get clear. Also the paper says that if the output sequence has $n$ elements and there are possibly $y$ lables, then there are $n^y$ possible label sequences. But I have some confusion. Lets say I have two elements and the possible labels are $[A, B, C]$, then the possible sequences are:               A A     B B     C C      A B     B A     A C      C A     B C      C B      which is equal to $9$ instead of $2^3=8$. So can anyone please provide some insights?


title: Lucene-based text feature construction

text: When doing the feature construction for text mining, does Lucene has a better performance in terms of classification/clustering result than the traditional bag-of-word approach?


title: Logistic regression for multiclass

text: I got the model for the logistic regression for multiclass which is given by $$ P(Y=j|X^{(i)}) = \frac{\exp(\theta_j^TX^{(i)})}{1+ \sum_{m=1}^{k}\exp(\theta_m^T X^{(i)})} $$ where k is the number of classes theta is the parameter to be estimated j is the jth class Xi is the training data Well one thing I didn't get is how come the denominator part $$ 1+ \sum_{m=1}^{k}\exp(\theta_m^T X^{(i)}) $$ normalized the model. I mean it makes the probability stay between 0 and 1. I mean I am used to logistic regression being $$ P(Y=1|X^{(i)}) = 1/ (1 + \exp(-\theta^T X^{(i)})) $$ Actually, I am confused with the nomalization thing. In this case since it is a sigmoid function it never lets the value be less than 0 or greater than 1. But I am confused in the multi class case. Why is it so? This is my reference https://list.scms.waikato.ac.nz/pipermail/wekalist/2005-February/029738.html. I think it should have been to be normalizing $$ P(Y=j|X^{(i)}) = \frac{\exp(\theta_j^T X^{(i)})}{\sum_{m=1}^{k} \exp(\theta_m^T X^{(i)})} $$


title: Independent gamma and normal distribution

text: Let $Z = X + Y$ where $X \sim N\left(\mu, \sigma^2 \right)$ and $Y \sim \Gamma\left(k, \theta \right)$ using this parametrization of the Gamma distribution. Also assume $X$ and $Y$ are independent. Then what is the distribution (pdf) of $Z$? * * * I thought this question looked straightforward, but it doesn't seem so. For example, I have tried using the convolution formulas here, but can't seem to find a closed form expression for the integral. I have also tried multiplying the moment generating functions (mgfs) of $X$ and $Y$, but it does not seem to match up to any known mgfs. Is it possible to find a closed-form solution for the distribution of $Z$? If not, how about some sort of approximation?


title: PCA with on data with unequal length

text: I am wondering how to perform PCA on data with unequal lengths. In the simplest case you could split up the data matrix so that one block would be $T_{1}Xn$ and the rest is $T_{2}Xn$. The second block has effectively missing data in the first $T_{1}-T_{2}Xn$ block. Many of the applications of PCA on missing data are variants of the EM algorithm methodology (these typically will assume the distribution is multivariate Gaussian) that are more suited to applications where the data is missing at random. Oftentimes, these techniques do not particularly apply to data with uneven lengths. However, there are techniques to estimate the mean and covariance matrix in these cases and it is possible to apply PCA to the resulting covariance matrix. However, when the number of variables becomes large (absolutely or relative to the number of observations), then it may no longer be convenient to estimate the covariance matrix in this fashion.


title: Show that if $X\ge 0$ , $E(X)\le \sum_{n=0}^{\infty}P(X>n)$

text: If $X$ is a random variable and also let $X\ge 0$. I want to show $E(X)\le \sum_{n=0}^{\infty}P(X>n)$.


title: Regression for two unequal samples

text: I have been working on a research theme entitled,"Service Quality and Shareholders Value in select Commercial Banks of India". I have collected the data from the sampled 1200 customers from all the banks, under study and I am using time series data of 7 years for my dependent variable i.e. shareholders value. The problem I am confronted with is how to assess the link of service quality having 1200 responses and shareholders value with just 7 entries using regression, if any.


title: Procedure for fitting an ARMA/GARCH Model

text: I want to try fitting an ARMA/GARCH model but want a methodological approach rather than fitting different models and picking the best one. However, I'm not sure how to choose my AR and MA terms for my mean equation, same thing for my variance equation. Hence I have the following queries:   1. Are there any conventions/guidelines for choosing these by looking at the ACF/PACF of returns and squared returns?   2. How are the coefficients estimated?Are the coefficients of the ARMA estimated first, then a GARCH fitted to the residuals?or are they estimate in parallel?   3. Also I read on wikipedia that the Ljung-Box test should be applied to the residuals of an ARIMA model to detect the presence of heteroscedasticty, but I saw also some people applying it on squared returns. Are these two the same?If not, what's the difference between the two?


title: ARIMA unusual forecasts

text: I am trying to implement the ARIMA model using statsmodels. I get really unusual results for my predictions, and was hoping for advice on fixing this.               arima = tsa.ARIMA(train[endogenous], exog=train.drop(endogenous,axis=1), order=(2,2,0),freq='B')     results = arima.fit()     prediction = results.predict(start=1,end=len(x)-1,exog=x.drop(endogenous,axis=1))      My actual dataset is this               2012-01-05    659.010     2012-01-06    650.020     2012-01-09    622.940     ...     2013-11-08    1016.03     2013-11-11    1010.59     2013-11-12    1011.78     2013-11-13    1032.47      Prediction gives me this               2012-01-05   -10.551134     2012-01-06    -8.937889     2012-01-09   -27.941221     ...     2013-11-08    14.739148     2013-11-11    22.567270     2013-11-12     1.844993     2013-11-13   -42.794671      It's unusual how even on examples that I trained on, the predictions aren't even in the same ballpark. When running on order (0,0,1), I got something closer               2011-06-01    525.632183     2011-06-02    528.051265     2011-06-03    523.168162     ...     2013-11-08     609.294394     2013-11-11     613.877047     2013-11-12     607.953117     2013-11-13    1033.130417      In this case, the prediction did really well on the samples that it was trained on, but horribly on the test set. EDIT: I tried looping through a bunch of different orders but it had negligible impact on the end result. I re-ran my model on order (1,1,2). Here were the coefficients:               const                              6.083430e+00     A                                  4.890720e-03     B                                  1.484685e-01     C                                  2.370248e-04     D                                  4.267728e-01     E                                 -4.398669e-08     F                                  5.203415e-01     G                                  5.349448e-03     H                                  2.198789e-03     I                                 -1.088271e-03     J                                 -3.295103e-03     K                                 -8.486572e-03     L                                  2.227050e-03     M                                 -1.844672e-02     N                                  1.561562e-04     O                                  8.622248e-03     P                                  1.560905e-05     Q                                 -1.705549e-02     R                                 -8.615488e-04     S                                  4.247976e-03     T                                  2.029798e-02     ma.L1.ENDOGENOUS                   1.088268e-01      Here was a summary:                                             ARMA Model Results                                   ==============================================================================     Dep. Variable:             ENDOGENOUS   No. Observations:                  449     Model:                     ARMA(0, 1)   Log Likelihood                 130.352     Method:                       css-mle   S.D. of innovations              0.181     Date:                Thu, 28 Nov 2013   AIC                           -214.705     Time:                        10:18:46   BIC                           -120.243     Sample:                    06-01-2011   HQIC                          -177.470                              - 02-18-2013                                              ===================================================================================================                                           coef    std err          z      P>|z|      [95.0% Conf. Int.]     ---------------------------------------------------------------------------------------------------     const                               6.0834        nan        nan      0.000           nan       nan     A                                   0.0049      0.005      0.931      0.352        -0.005     0.015     B                                   0.1485      0.030      4.902      0.000         0.089     0.208     C                                   0.0002   8.36e-05      2.834      0.005      7.31e-05     0.000      D                                   0.4268      0.002    216.200      0.000         0.423     0.431     E                              -4.399e-08        nan        nan      0.000           nan       nan     F                                   0.5203      0.001    356.889      0.000         0.517     0.523     G                                   0.0053      0.003      1.773      0.077        -0.001     0.011     H                                   0.0022      0.002      1.318      0.188        -0.001     0.005     I                                  -0.0011      0.000     -2.421      0.016        -0.002    -0.000     J                                  -0.0033      0.002     -1.937      0.053        -0.007  3.95e-05     K                                  -0.0085      0.002     -4.968      0.000        -0.012    -0.005     L                                   0.0022      0.001      3.398      0.001         0.001     0.004     M                                  -0.0184      0.003     -5.684      0.000        -0.025    -0.012     N                                   0.0002        nan        nan      0.000           nan       nan     O                                   0.0086      0.000     28.311      0.000         0.008     0.009     P                                1.561e-05        nan        nan      0.000           nan       nan     Q                                  -0.0171      0.001    -16.801      0.000        -0.019    -0.015     R                                  -0.0009      0.001     -0.900      0.369        -0.003     0.001     S                                   0.0042      0.007      0.570      0.569        -0.010     0.019     T                                   0.0203      0.013      1.621      0.106        -0.004     0.045     ma.L1.ENDOGENOUS                    0.1088      0.048      2.284      0.023         0.015     0.202                                         Roots                                         =============================================================================                      Real           Imaginary           Modulus         Frequency     -----------------------------------------------------------------------------     MA.1           -9.1889           +0.0000j            9.1889            0.5000     -----------------------------------------------------------------------------


title: Please help with this expectation problem

text: Suppose $y_1$, $y_2$, $\cdots$, $y_n$ are sampled from a normal distribution with expected value $\mu$ and variance $\sigma$. Then how to obtain the expectations and variances of the following expressions of $\mathcal{F}$ and $\mathcal{G}$? ![enter image description here](http://i.stack.imgur.com/JAelX.jpg) where $\rho$ is a constant $\in (0, 1)$, and $c_i, i=1, 2, \cdots$ are also constants. $|D_k|<1$ and $|D_{k-1}|<1$. Actually, this expression was involved in my experiments, and from the experimental results, it seems the expectations should be 0. But I couldn't prove it mathematically, and I have been working on these expressions for a couple of weeks. Please advise! Any suggestions are highly appreciated! Thank you!


title: How to know if my data fits Pareto distribution?

text: I have a sample which is a vector with 220 numbers. Here is a link to a histogram of my data.. And I wish to check if my data fits a Pareto distribution, but I don't want to see QQ plots with that distribution, but I need an exact answer with p-value in R, such as Anderson-Darling test for normality (`ad.test`). How could I do that? Please be as specific as you can.


title: A Question from textbook "Learning with Kernels"

text: I am reading the book "Learning with Kernels: Support Vector Machines, Regularization, Optimization, and Beyond (Adaptive Computation and Machine Learning)". I finished the first chapter and didn't notice any approach taken by the authors, but I assume according to the material in the first chapter it is a frequentist approach. So my question is from the people who has read the book or know about it. Am I right about this assumption?


title: Longitudinal Study into attention bias and insomnia - What test to use?

text: I am doing a study into attentional bias (as measured by the Stroop test) and insomnia (as measured by the Insomnia Severity Index - ISI). I will be comparing the mean reaction times of the Stroop test to the mean ISI scores at 4 different time points (Baseline, 3 month follow up, 6 month follow up and 9 month follow up). If I was looking at the relationship between the 2 means at one time point, I would do a Pearson r, but i know I cannot use this longitudinally as it violates Pearson r assumptions. How do I look at the relationship between the 2 scores over time? Any help GREATLY appreciated!


title: Confusion in Expectation Propagation energy function

text: I have two questions in Energy function of expectation propagation. 1. Seeing the video of [1] (slide 24 of the second part), Minka says that, the evidence is the following: $$ Z = \left( \int q(\mathbf{x}) d\mathbf{x} \right) ^{1-n} \prod_{i=1}^{n} \int \frac{f_i (\mathbf{x})}{\tilde{f}_i(\mathbf{x}) } q(\mathbf{x}) d\mathbf{x} \quad \text{ (Equation a)} $$ But in [2] he introduces the primal energy function as it is shown in equation (6) of [2]. $$ \sum_i \int_{\mathbf{x}} \hat{p}_i(\mathbf{x}) \log \frac{ \hat{p}_i(\mathbf{y}) }{ f_i(\mathbf{x}) p(\mathbf{x}) } d\mathbf{x} - (n-1) \int_\mathbf{x} q(\mathbf{x}) \log \frac{q(\mathbf{x})}{p(\mathbf{x})} d\mathbf{x} \quad \text{ (Equation b)} $$ with $$ \hat{p}_i = \frac{1}{Z} q^{\setminus i} f_i(\mathbf{x}) $$ As I understand, we are trying to minimze $-\log Z$ which must be the same as equation (6) of [2]: $$ \log Z = \sum_{i=1}^{n} \log \int \frac{f_i (\mathbf{x})}{\tilde{f}_i(\mathbf{x}) } q(\mathbf{x}) d\mathbf{x} - (n-1) \log \int q(\mathbf{x}) d\mathbf{x} \quad \text{ (Equation c)} $$ However, I don't see how these are the same ( **UPDATE:** I mean (Equation b) and (Equation c) ). Any idea on this?   1. At [3] after equation 64 it says that : "The scale $s_a$ that minimizes local $\alpha$-divergence is ... ". Let's say we are in the case of KL divergence, and want to find $s_a$ for $\alpha=1$ (correct?). Since $s_a$ is defined to normalize the distribution $\tilde{f}_a(\mathbf{x})$, thus: $$ s_a = \frac{1} {\int \tilde{f}_a(\mathbf{x}) } $$ (right?) and $$ \tilde{f}_a(\mathbf{x})^{new} = \frac{ \arg\min_{q \in \mathcal{Q}} KL \left( q(\mathbf{x}) \frac{f_a(\mathbf{x})}{ \tilde{f}_a(\mathbf{x}) } || q^{\prime} \right) }{ q(\mathbf{x}) \frac{1}{ \tilde{f}_a(\mathbf{x}) } } \quad \text{ (Equation d)} $$ The paper [3] says that (for $\alpha=1$): $$ s_a = \frac{ \int q(\mathbf{x}) \frac{f_a(\mathbf{x})}{ \tilde{f}_a(\mathbf{x}) } d\mathbf{x} }{ \int q(\mathbf{x}) d\mathbf{x} } \quad \text{ (Equation e)} $$ However I don't see why. Any hints?  [1] videolectures.net/mlss09uk_minka_ai/ [2] http://research.microsoft.com/en-us/um/people/minka/papers/ep/minka-ep- energy.pdf [3] http://www.seas.harvard.edu/courses/cs281/papers/minka-divergence.pdf


title: Multidimensional scaling in R with Spearman's rank correlation

text: I'm new to MDS, but I found some good starter code here (http://mhermans.net/static/postdata/r-examples/neighbours-mds/neighbours-mds- example.html) and which I borrowed for below. The code works fine, but I would like to change the way I calculate my dissimilarity matrix. The data I want to analyze is ecological time series, so I'd like to use a Spearman's Rank Correlation, as opposed to the co-occurrence matrix outlined from the above site. So given a simple data set of 10 species from 6 sites:               Site_1 <- c(1, 1, 1, 1, 1, 2, 2, 4, 8, 16)     Site_2 <- c(3, 4, 3, 2, 3, 5, 6, 10, 17, 33)     Site_3 <- c(5, 4, 5, 6, 5, 9, 8, 18, 33, 66)     Site_4 <- c(16, 15, 16, 14, 16, 15, 17, 13, 10, 1)     Site_5 <- c(5, 6, 5, 7, 5, 5, 4, 5, 3 ,5)     Site_6 <- c(16, 9, 3, 3, 2, 2, 1, 1, 1 ,1)          DF <- data.frame(Site_1, Site_2, Site_3, Site_4, Site_5, Site_6)     row.names(DF) <- c('S1','S2','S3','S4','S5','S6','S7',"S8",'S9','S10')      The code that works from the above site looks like this:               X <- as.matrix(DF)  # convert to matrix     C <- t(X) %*% X     # turn dummies into co-occurence matrix     D <- max(C) - C     # subtract max value to turn into dissimilarities     diag(D) <- 0                                # clear values on diagonal     Ds      <- as.dist(D, diag = T)             # convert to distance-object     fit     <- cmdscale(Ds, eig = TRUE, k = 2)  # MDS with 2 dimensions          x <- fit$points[, 1]     y <- fit$points[, 2]          plot(x,y)     text(x,y,labels=names(DF))      I was hoping I could change the dissimilarity matrix like so and still get a working MDS:               C   <- cor(DF, method="spearman")       #Calculate Spearman Rank Correlation     Ds  <- as.dist(C, diag = T)             # convert to distance-object     fit <- cmdscale(Ds, eig = TRUE, k = 2)  # MDS with 2 dimensions          x <- fit$points[, 1]     y <- fit$points[, 2]          plot(x,y)     text(x,y,labels=names(DF))      But the resulting plot no longer groups the sites according to a pattern of distances that makes sense to the original data. So I'm assuming I can't just substitute a Spearman's correlation in, and more globally, that I don't understand the difference between a co-occurrence matrix and correlation matrix. I hope you can help me out plotting a MDS that makes sense with a Spearman's correlation.


title: Comparing continuous predictors for a dichotomous variable

text: I have two continuous predictor variables to predict a dichotomous variable. In addition i have constructed two (interaction) models, based on domain knowledge which use both variables to predict the third. So now i want to compare these predictors, using R. I get these result on the wilcox test: W = 36655, p-value = 3.896e-09 (single predictor) W = 29680.5, p-value < 2.2e-16 (model) But i'm wondering if the area under the ROC curve would not be a good (better?) measure then the wilcox or the t-test. To get a general idea of how well the continuous predictor separates the two groups i plot the two histograms together. But the differences are small and besides, i'd like a 'grade' to tell me which is better and by how much.


title: Is it possible for quartiles to have uneven groups?

text: If I have a data set with many values repeated (example below) how would it be divided into quartiles? I assumed at first that the quartiles would be uneven, but every single definition I've read makes it sound like they must be equal. Example:               { 1, 1, 4, 4, 4, 4, 8, 8}


title: Can RSD be pooled or is it valid to use the uncertainty propagation rule?

text: This might be a strange case, but I'm sure that somebody in here can help me and that I'm perhaps not the only one who wants to pool RSDs correctly. Consider this: The mean weight of a powder bed is `W` mg with a relative standard deviation (RSD) of `Wrsd` %. For each weighted powder bed a precise counting of the particles were done and the mean particle mass is calculated to be `PM` mg/particle and the RSD is found to be `PMrsd` %. Can these two RSD (`Wrsd` and `PMrsd`) be pooled together and thereby be used to estimate the RSD of the number of particles within a weighted powder bed? If so, what does the equation look like? I have tried to treat the RSD as relative uncertainties and thereby using the propagation rule, $e_3 = \sqrt{(e_1)^2 + (e_2)^2 }$ Is this approach violating anything? ## UPDATE (@jonsca thank you for updating my question so it is more readable) (@whuber thank you for your interest in my question) A more detailed description: A powder die was filled multiple times. Each time it was weighted and the particles counted. Putting it in a table looks like this: ![Table with data](http://i.stack.imgur.com/xXTO1.png) Now this could have been done in a different way: I could have taken multiple collections of particles and determined the mean particle mass. So ex. take 7 times 100 random particles (no die involved) and determine the mean particle mass from these measurements. I could then just weighted some new powder beds (same type of particles) within the die and divide that with the mean particle mass to estimate the mean number of particles that fits within the die. I don't just want to know the mean number of particles that fits within the die, but are also interested in the span of this number. Using the latter method I would expect the span to be equal to the RSD of the powder bed mass (`Wrsd`), which clearly the counting shows are too small. What I notices was that combining the two RSDs (powder bed and particle mass) gives a quite "good" estimate (7.3 %) of the RSD found by the counting (8.7%). Does this make sense? [If you read this far, then thank you for you interest :D]


title: SPSS dummy variables in OLS

text: I have a timeseries dataset holding stock data for a large set of companies. Assume the following subset, where `obsDay` is the observation day (148 days in reality) and `weekDay` the day of the week (1 = monday, 2 = tuesday, ...):               company  sector  obsDay weekDay  stockPrice     -------------------------------------------     1        15      1      3        10.40     1        15      2      4         9.42     1        15      3      5         9.66     1        15      4      1        11.00     1        15      5      2        10.21     2        10      1      3        43.55     2        10      2      4        43.50     2        10      3      5        40.31     2        10      4      1        48.43     2        10      5      2        43.00     3        20      1      3        10.00     3        20      2      4        11.00     3        20      3      5        12.00     3        20      4      1        13.00     3        20      5      2        14.00      In an OLS regression, I would like to see the effect, for instance, of day of the week and industry sector, on stockprice. It is suggested that I would then have to make dummies for each variable, as follows (for sector):               company  sector  sector15  sector10  sector20 obsDay weekDay  stockPrice     ------------------------------------------------------------------------     1        15      1         0         0        1      3        10.40     1        15      1         0         0        2      4         9.42      (and the same for weekDay). But with the many variables I have, this is making the dataset hard to read. Can't SPSS create dummy variables from the `sector` and `weekDay` categories on the fly when performing the OLS?


title: Finding the fitted and predicted values for a statistical model

text: Let's say I have the following data and am running a regression model:               df=data.frame(income=c(5,3,47,8,6,5),                   won=c(0,0,1,1,1,0),                   age=c(18,18,23,50,19,39),                   home=c(0,0,1,0,0,1))      On one hand, I run a linear model to predict on income:               md1 = lm(income ~ age + home + home, data=df)      Second, I run a logit model to predict on the won variable:               md2 = glm(factor(won) ~ age + home, data=df, family=binomial(link="logit"))      For both models, I wonder how I can generate a table or data frame with the predictor response category, fitted value, and the model predicted value. So for the linear model, something like:               age  fitted_income  predicted_income     18    3              5      23    3              3     50    4              2     19    5              5     39    6              4          home   fitted_income    predicted_income     0       5               6            1       3               9      Or perhaps it should be for each data point. So for x_i data point, the fitted and predicted values are:               id   age  fitted_income  predicted_income     1     18    3              5      2     23    3              3     3     50    4              2     4     19    5              5     5     39    6              4        1. From a statistical standpoint, is such an undertaking useful? Why or why not?   2. How can this be done in R? (looked at names(md1) and found what I can pull from the model, but haven't proceeded past that) Thanks!


title: Getting p-values from a mixed-effect ANOVA

text: I'm struggling to extract F-values and/or p-values from an ANOVA with all random effects or an ANOVA with both mixed and random effects. I would expect to be able to get F- and p-values regarding the variance of the random effects based on the expected mean squares. Maybe the `aov()` function just doesn't return what I'm looking for? Or I'm specifying my model incorrectly? Here's a simple example where I want to have `cyl` as a fixed effect and `carb` as a random effect:               data(mtcars)     summary(aov(formula = mpg ~ cyl + Error(carb), data = mtcars))      But no p-values for the random effect are reported.               Error: carb         Df Sum Sq Mean Sq     cyl  1  341.8   341.8          Error: Within               Df Sum Sq Mean Sq F value  Pr(>F)         cyl        1  492.1   492.1   48.85 1.1e-07 ***     Residuals 29  292.2    10.1                         ---     Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1


title: Problem with calculating Asymptotic Standard Error for Somer's D

text: I'm trying to calculate Asymptotic Standard Error (ASE) for Somer's D measure of association for some population data. The formula for variation (sd^2) I've found (in the article specified below) is: ![enter image description here](http://i.stack.imgur.com/oYpdw.jpg) ![enter image description here](http://i.stack.imgur.com/XYTDh.jpg) (this is identical with equations presented in _"Analysis of Ordinal Categorical Data" by Alan Agresti_ ) I've tried to implement it in R:               pairsC <- function(x,prob=FALSE) {            l <- nrow(x) # rows       k <- ncol(x) # columns            pc <- 0              for (i in 1:l)         for (j in 1:k) {                row.idx <- (1:l)[which((1:l)>i)]           col.idx <- (1:k)[which((1:k)>j)]                if (length(row.idx)>0 && length(col.idx)>0)                  pc <- pc+x[i,j]*sum(x[row.idx,col.idx])              }            ifelse(prob,2*pc,pc)          }          pairsD <- function(x,prob=FALSE) {            l <- nrow(x) # rows       k <- ncol(x) # columns            pd <- 0            for (i in 1:l)         for (j in 1:k) {                row.idx <- (1:l)[which((1:l)>i)]           col.idx <- (1:k)[which((1:k)<j)]                if (length(row.idx)>0 && length(col.idx)>0)                  pd <- pd+x[i,j]*sum(x[row.idx,col.idx])              }            ifelse(prob,2*pd,pd)          }          prob <- TRUE          pc <- function(x,i,j) {            pcx <- 0            l <- nrow(x) # rows       k <- ncol(x) # columns            for a in 1:l         for b in 1:k                 if (a<i and b<j) pcx <- pcx+x[a,b]            for a in 1:l         for b in 1:k                if (a>i and b>j) pcx <- pcx+x[a,b]            pcx          }          pd <- function(x,i,j) {            pdx <- 0            l <- nrow(x) # rows       k <- ncol(x) # columns            for a in 1:l         for b in 1:k                 if (a<i and b>j) pdx <- pdx+x[a,b]            for a in 1:l       for b in 1:k            if (a>i and b<j) pdx <- pdx+x[a,b]            pdx          }          fi <- function(x,i,j) {            l <- nrow(x) # rows       k <- ncol(x) # columns            tmp <-              -2 * sum(x[i,]) * (pairsC(x,prob)-pairsD(x,prob)) +              -2 * (1-sum(x[setdiff(1:l,i),])^2) * (pc(x,i,j)-pd(x,i,j))              ifelse(length(tmp>0),tmp,0)          }          # nom1 & nom2 - nominators          l <- nrow(x) # rows     k <- ncol(x) # columns          nom1 <- 0     nom2 <- 0          for (i in 1:l)       for (j in 1:k) {              nom1 <- nom1 + x[i,j]*fi(x,i,j)^2              nom2 <- nom2 + x[i,j]*fi(x,i,j)            }          nom2 <- nom2^2          # den - denominator          den <- (1-sum(rowSums(x)^2))^2          # ASE          sqrt(nom1-nom2)/den      The table I'm trying to calculate ASE for is:               > portfolio                      [,1]   [,2]  [,3]   [,4]  [,5]   [,6]   [,7]     defaulted    0.000025 0.0001 0.001 0.0025 0.004 0.0075 0.0125     nondefaulted 0.049975 0.0999 0.199 0.2475 0.196 0.1425 0.0375       According to an article I've taken the table from ( _Oliver Blumke "A proposal for a validation methodology for the discriminatory power of a rating system over time", The Journal of Risk Model Validation, Vol. 5/Num. 1, Spring 2011_ ), the Somer's D should be 0.684 and ASE 2.378. While Somer's D I've calculated is identical to the one in the article, less the sign, i.e. -0.684 instead of +0.684, "my" ASE is 7.5468. Can anyone help me with deriving the correct ASE value? Thanks in advance.


title: Integrating a Gaussian Distribution multiplied by 1/x^(2/3), any ideas?

text: I am trying to integrate the following function over $x$ and $g$: $f(z)= \int_{a/z}^{d} \int_{a}^{zg} \frac{e^{-(g-\mu)^2/2\sigma^2}}{x^{2/3+1}} dx~dg$ When I integrate over $x$ I get: $f(z)= 3 \int_{a/z}^{d} {e^{-(g-\mu)^2/2\sigma^2}}((gz)^{-2/3}-a^{-2/3})dg$ After re-arrangement we get two terms, one of the is easy to solve, the other one is this: $f(z)= 3z^{\frac{-2}{3}} \int_{a/z}^{d} {e^{\frac{-(g-\mu)^2}{2\sigma^2}}}g^{\frac{-2}{3}}dg$ At this point I am not able to solve the integration. Any Idea how to do that, I did it in Mathematica and Matlab but I came up empty. Any ideas recommendations will be highly appreciated.


title: How to measure/clarify a possible relationship between two spatial points using R

text: I'm quite new to the area of spatial statistics, but I'm very interest in some general principles. The last two weeks I've created an example dataset, which contains three datsets.   1. A dataset of ill persons.    2. A dataset of cities with the the overall population.   3. A dataset of points, which visualizes waterfeatures. The whole situation looks like this:   Blue: Water features, yellow: cities, red: persons. Please note that the persons are located on the cities coordinates. ![PLot](http://i.stack.imgur.com/iFi21.png) I've already performed some basic analyses: Person dataset:   * I've calculated the distance between each person and the nearest waterfeature. City dataset:   * I've calculated the number of ill and not_ill people per city.   * I've calculated the rate of ill and not_ill people   * Because persons and cities share the same location, I've also infer the distance between each city and the nearest water feature. Now I want to check a possible correlation between the number of ill persons/rate of ill persons and the proximity to water features. I know, that the datasets are possibly not representative or suitable for my hyptothesis: I hold that there are more ill persons, where the distance to water features are lower than somewhere else. I've already looked for some suitable methods, but there are so many possible ways, so I don't know, which of these could be useful for my notional use case. I've read about semivariogram, variogram, Ripley's K function, G-Function, correlation coefficient. As you can see, I have a broad range of _possible_ methods, but unfortunately not this necessary expertise.   My questions:   1. Do you have any tip, which methods could be the most suitable?    2. And another question: A risk/cluster analysis would be cool. Something, which shows the "areas" with a high risk of becoming ill based on the number of ill persons in a city. But I think there are two problems: I have to interpolate my dataset!? And a cluster analysis is only possible by using polygons, right? Interesting R packages (from my point of view) could be `SpatialEpi` or `DCluster` as well as `spatsat`.   3. And one last question: is there any "general" problem of having many persons at one location? I know that this is not the best dataset, because ideally each point has a own position... To give you a better overview, I've prepared some code in R, which loads my dataset and which plots the picture, which I've included here.               library(RgoogleMaps)     library(ggplot2)     library(ggmap)     library(sp)     library(fossil)          persons = read.csv("http://pastebin.com/raw.php?i=AuAQNqVt", header = TRUE, stringsAsFactors=FALSE)     city= read.csv("http://pastebin.com/raw.php?i=ZfPDFYCK", header = TRUE, stringsAsFactors=FALSE)     water= read.csv("http://pastebin.com/raw.php?i=hQRvMZwE", header = TRUE, stringsAsFactors=FALSE)          # plot data     gc <- geocode('new york, usa')     center <- as.numeric(gc)       G <- ggmap(get_googlemap(center = center, color = 'bw', scale = 1, zoom = 11, maptype = "terrain", frame=T), extent="device")     G1 <- G + geom_point(aes(x=POINT_X, y=POINT_Y ),data=city_all_parameters, shape = 22, color="black", fill = "yellow", size = 4) + geom_point(aes(x=POINT_X, y=POINT_Y ),data=persons, shape = 8, color="red", size=2.5) + geom_point(aes(x=POINT_X, y=POINT_Y ),data=water, color="blue", size=1)     plot(G1)      I'm excited to hear from you :)


title: Comparing two sets of clusters

text: I have some data which I wish to cluster. I would like to see how these clusters compare to the categories that have already been assigned. Is there some kind of metric or visualisation that will tell me how well these line up?   * Any suggestions with each item assigned to just one category is fine, but each data element can actually be assigned multiple categories. Perhaps there's a neat way to deal with this?


title: Why does the EM algorithm have to be iterative?

text: Suppose that you have a population with $N$ units, each with a random variable $X_i \sim \text{Poisson}(\lambda)$. You observe $n = N-n_0$ values for any unit for which $X_i > 0$. We want an estimate of $\lambda$. There are method of moments and conditional maximum likelihood ways of getting the answer, but I wanted to try the EM algorithm. I get the EM algorithm to be $$ Q\left(\lambda_{-1}, \lambda\right) = \lambda \left(n + \frac{n}{\text{exp}(\lambda_{-1}) - 1}\right) + \log(\lambda)\sum_{i=1}^n{x_i} + K, $$ where the $-1$ subscript indicates the value from the previous iteration of the algorithm and $K$ is constant with respect to the parameters. (I actually think that the $n$ in the fraction in parentheses should be $n+1$, but that doesn't seem accurate; a question for another time). To make this concrete, suppose that $n=10$, $\sum{x_i} = 20$. Of course, $N$ and $n_0$ are unobserved and $\lambda$ is to be estimated. When I iterate the following function, plugging in the previous iteration's maximum value, I reach the correct answer (verified by CML, MOM, and a simple simulation):               EmFunc <- function(lambda, lambda0){       -lambda * (10 + 10 / (exp(lambda0) - 1)) + 20 * log(lambda)     }          lambda0 <- 2     lambda  <- 1          while(abs(lambda - lambda0) > 0.0001){       lambda0 <- lambda       iter    <- optimize(EmFunc, lambda0 = lambda0, c(0,4), maximum = TRUE)       lambda  <- iter$maximum     }          > iter     $maximum     [1] 1.593573          $objective     [1] -10.68045      But this is a simple problem; let's just maximize without iterating:               MaxFunc <- function(lambda){       -lambda * (10 + 10 / (exp(lambda) - 1)) + 20 * log(lambda)     }          optimize(MaxFunc, c(0,4), maximum = TRUE)     $maximum     [1] 2.393027          $objective     [1] -8.884968      The value of the function is higher than in the un-iterative procedure and the result is inconsistent with the other methodologies. Why is the second procedure giving a different and (I presume) incorrect answer?


title: Correlation between variables over time

text: I have two variables over two different years. Week by week pageview numbers in 2013 and 2014 and week by week # of posts in 2013 and 2014 (number of articles posted). I'm trying to find out if posting less frequently has caused a dip in pageviews. Right now I am taking the percentage difference for pageviews and posts and looking for correlation between the week over week percentage differences. Can anyone help with me a more accurate analysis?


title: K-means metrics

text: I have read through the scikit learn documentation and Googled to no avail. I have 2000 data sets, clustered as the picture shows. Some of the clusters, as shown, are wrong, here the red cluster. I need a metrics to method to validate all the 2000 cluster-sets. Almost every metric in scikit learn requires the ground truth class labels, which I do not think I have or CAN have for that matter. I have the hourly traffic flow for 30 days and I am clustering them using k-means. The lines are the cluster centers. What should I do? Am I even on the right track? ![enter image description here](http://i.stack.imgur.com/JupUk.png)


title: Exposure classification bias, post-modeling

text: I am building a high-dimensional Bayesian spike-and-slab model to study the association between several organic compounds and a continuous outcome. The goal is to select the most influential compounds that are most explanatory for the continuous outcome. A number of compounds in my data are unclassified, i.e. they do not have a known nomenclature. When I completed the models, some of these unclassified compounds were selected as influential variables. I cannot let the unclassified compounds just stay in my results because they will be tough to interpret by the downstream audience of this analysis. So I am now thinking of running the unclassified compounds that made it to the influential shortlist through another database and retrieve nomenclatures that the original database did not assign. The two ways of classification are fundamentally different; one assigns nomenclature by relying on clustering and structural similarity and another compares it to standards. I can think of a tricky issue that might crop up here. There might be some compounds which are labeled as say, A by the new database. If the old database had assigned some compounds the tag A, then I have a situation where some of the compounds tagged A are influential and some are not. How do I deal with this bias? Any ideas on how I can structure my thoughts about this problem are most welcome.