title: Are Product KDEs proper density functions?

text: I have a fairly simple question about (multivariate) Kernel Density Estimators, but I somehow don't seem to find the answer anywhere: are these estimators supposed to specify a proper probability density function (integrating to 1.0)? I am currently implementing a small software package including kernel density estimators, and in particular a product kernel density estimator for the multivariate case, using the formula found here, and Scott's rule of thumb for the bandwidth selection. The estimator seems to get right the overall shape of the distribution, but the density values are clearly not classical probability densities, as their integration would yield a value far superior to 1. So I am wondering if I made an error in my implementation, or if these values are correct. I have a follow-up question in case the estimator is not supposed to yield regular probability densities: is there an easy way to convert the estimator function into a standard PDF (i.e. to normalise the function)?


title: Dirichlet process mixture model with Bayesian hierarchical clustering

text: I am doing Bayesian hierarchical clustering. From my understanding, there are three basic points for this algorithm.   1. Use marginal likelihoods to decide which clusters to merge   2. Asks what the probability is that all the data in a potential merge were generated from the same mixture component. Compare to exponentially many hypotheses at lower levels of the tree   3. Generative model used is a Dirichlet Process Mixture Model (DPM) The first two points are easy to understand, but why the generative model is DPM? Are there specific characteristics of DPM make it quite fit for this purpose? Are there any other options?


title: Splitting a numeric column for a dataframe

text: I have a dataframe df as shown below                 name  position                 1 HLA   1:1-15                 2 HLA   1:2-16      3 HLA   1:3-17               I would like to split the position column into two more columns based on the ":" character such that i get               name    seq    position                 1 HLA   1       1-15                 2 HLA   1       2-16      3 HLA   1       3-17       So i thought this would do the trick,               df <- transform(df,pos = as.character(position))           df_split<- strsplit(df$pos, split=":")          #found this hack from an old mailing list post     df <- transform(df, seq_name= sapply(df_split, "[[", 1),pos2= sapply(df_split, "[[", 2))      however I get an error               Error in strsplit(df$pos, split = ":") : non-character argument      What could be wrong? How do you achieve this in R. I have simplified my case here, actuality the dataframe runs to over a hundred thousand rows.


title: Mixed-model specification in R - exploring individual differences

text: My current experiment aims to explore the effect of viewing condition, difficulty and depth perception on object grasping. I'm mostly interested in the effect of individual differences in depth perception. I've previously been specifying my model as a mixed-factor ANOVA, with viewing condition (binocular, monocular) and difficulty (5 levels) as repeated- measures factors and ability to perceive depth as a between subjects factor (raw data for this was logarithmic, I separated the participants into two groups ['Good' vs 'Poor'] depending on which side of a threshold their scores fell on). Since I'm interested in individual differences, the idea of using mixed-models is appealing, but I can't quite grasp how to explore the effect of depth perception using this technique. I understand that viewing condition and difficulty would be fixed factors. I think that both ID and depth perception would be random effects, so I can specify the model like so: `fit.1 <- lme(Accuracy ~ Difficulty*View, random=~1|ID/DepthPerception, data=mydata)` Am I correct in thinking that I can leave the depth perception data in its logarithmic form? Having fit this model, I cannot figure out how to interpret the effect of depth perception - only the fixed effects show up in the results.


title: Making sense out of statistics theory and applications

text: I have recently graduated with my masters degree on medical and biological modeling, accompanied with engineering mathematics as a background. Even though my education program included a significant amount of courses on mathematical statistics (see below for a list), which I managed with pretty high grades, I frequently end up completely lost staring down on both theory and applications of statistics. I have to say, compared to "pure" mathematics, statistics really makes little sense to me. Especially the notations and language used by most statisticians (including my past lecturers) is annoyingly convoluted and almost none of the resources I have seen so far (including wikipedia) had simple examples that one could easily relate to, and associate to the theory given... This being the background; I also realize the bitter reality that I cannot have career as a researcher/engineer without a firm grip on statistics, especially within the field of bioinformatics. I was hoping that I could get some tips from more experienced statisticians/mathematicians. How can I overcome this problem I have mentioned above? DO you know of any good resources; such as books, e-books, open courses (via iTunes or OpenCourseware for ex) etc.. Thanks, **EDIT:** As I have mentioned I am quite biased (negatively) towards a majority of the literature under general title of statistics, and since I can't buy a number of large (and expensive) coursebooks per branch of statistics, what I would need in terms of a book is something similar to what Tipler & Mosca is for Physics, but instead for statistics. For those who don't know about Tipler; it's a large textbook that covers a wide majority of the subjects that one might encounter during higher studies, and presents them each from basic introduction to slightly deeper in detail. Basically a perfect reference book, bought it during my first year in uni, still use it every once in a while. * * * The courses I have taken on statistics:   * a large introduction course,    * stationary stochastic processes,   * Markov processes,   * Monte Carlo methods   * Survival analysis


title: Effect size in contrast analysis

text: I'm running 2-way repeated measures ANOVA (3 and 5 levels) with planned contrasts afterwards. I'm interested in comparing 1)levels of the first factor among themselves (e.g. A1B with A2B, A2B with A3B) 2)pairwise comparisons (e.g. A1B1 and A2B4) Question: Some sources claim that I need to report effect size for each contrast. I wonder if this makes sense since the design is repeated measures. If so, how the contrasts should be calculated? I don't think the standard formula for Cohen's d works in this case since it doesn't take into account the correlation.


title: R plots to a connection?

text: Is there a way to have R plot to an in-memory object or connection, rather than a named file? I would like to have a plotting server create many graphs without ever going to a file. The Cairo package documents use of a connection, but it doesn't seem to work. What I would like to do is something like:               library(Cairo)     plot.to.var <- function(data) {        tc = textConnection("output", "w")        CairoPDF(tc)        plot(data)        dev.off()        tc.close()        output     }      When I do this, CairoPDF mentions a connection patch I can't find a reference to, that will allow me to do this, even though the documentation shows this working. I have no particular desire to use Cairo, merely saw that the documentation mentioned this. Any ideas?


title: Concept of a random linear model

text: I want to compare few methods of variable selection. I want to do it using simulations. I'm aware of the fact that it won't be an ultimate answer to question 'Which method is the best one?', but I'm looking just form some hint. To do such simulations I need a method to draw a 'random linear model'. Is there any well accepted algorithm of drawing 'random linear model'? By well accepted I mean method that was used for example in some scientific paper. I was thinking about following simple approach:   1) Choose $n$ and $k$, which denotes number of observations and number of variables.   2) Generate random matrix $X$ by drawing each element using uniform distribution $(0,1)$.   3) Generate parameters using uniform distribution $(0,1)$.   4) Generate residuals using Normal Distribution $(0,\sigma^2)$ for some fixed and arbitrary chosen $\sigma^2$   5) Calculate $Y=X\beta + \epsilon$


title: Algorithm for choosing the number of clusters when using pam in R?

text: I am clustering a dataset using the pam command (from {cluster} package), and I wish to decide on the number of clusters to use. I was able to implement The_Elbow_Method in R (see wiki) for doing that. But that doesn't provide me with any solid criteria (like AIC, for example) for decision. I came by the {clValid} package which looks promising, but I wanted to know if there are **any other R solutions (you know of) for choosing the number of clusters for pam?** Here's some dummy code if someone wants to show examples:               data(iris)     head(iris)     require(cluster)     pam(iris[,1:4], 3)


title: Hierarchical choice modelling

text: I am trying to model students choices on studying from year to year and was thinking that a HLM would be the way to go, however I have become stumped in how to progress. The variable I am trying to predict has 3 categories - Active, Completed or Discontinued. I want to find the effect of their class, how far into their study they are (in years) and whether they are full or part time students on this variable and asses what % of students fall under each category for each faculty, in each year (for Part time and full time students). This problem reminded me of the classic HLM example of test scores varying by the individual and across classrooms, so I have assumed the HLM would also apply in this case. If I have taken the first step in the wrong direction, please let me know To be painstakingly clear, the ideal result would be to have the probability of either staying (active), completing or discontinuing for each class and each year within that class. My attempts so far have been in R using the lme4 package and lmer. I have created the following example of my problem. I have tried to replicate the expected relationships, although only crudely - I am really just trying to show the format of data I am working with.               library(lme4)     class = rep(LETTERS[1:5], each = 20)     P_F_time = sample(c("P", "F"), 100, replace = TRUE)     year = sample(1:7, 100, replace = TRUE)     for(i in 1:100){         x = runif(1)         if(year[i] < 3){         choice[i] = "Active"         next     }     if(x < 0.25){         choice[i] = "Discont"     }else{         choice[i] = "Complete"     }     }     data = data.frame(class, year, choice, P_F_time)     fit = lmer(choice ~ P_F_time + year + (1 + P_F_time + year | class), data = data)      Does this model show the relationship I have described? If so - how would I view the results in a way that shows the relative choices for each faculty in each year for both full time and part time students? as mentioned above. Since this is an area i am not familiar with, please ask for more information or clarification if it is needed. Thanks, Sam


title: Parameter tuning in lars (lasso) matlab

text: I am trying to use lars (matlab implementation:http://www.ece.ubc.ca/~xiaohuic/code/LARS/LARS.htm). I want to do a leave one out cross validation on my data using this code. I have the following doubts: This code has a parameter t, is this the parameter which has to be tuned during cross validation? If it is not given, the default is taken to be infinity. And the coefficients (weight matrix) is calculated to be of arbitrary dimension. Can somebody help me understand whats going on in this code? I would like to perform a leave one out cross validation on my data using lars. **EDIT** I am also trying to use lasso function in-built with MATLAB. It has an option to perform cross validation. It outputs B and Fit Statistics, where B is Fitted coefficients, a p-by-L matrix, where p is the number of predictors (columns) in X, and L is the number of Lambda values. Now given a new test sample, how can I calculate the output using this model? Thanks


title: Weibull distribution with the negative shape parameter

text: Just wondering why in the literature Weibull distribution is always defined for positive shapes, whereas the extension in the negative direction is possible and has many useful properties. Suppose $X \propto \mathrm{Weibull}(\theta, \lambda)$, i.e. is Weibull- distributed with shape $\theta > 0$ and scale $\lambda > 0$ with PDF $$ p_X(x) = \frac{\theta}{\lambda}\left(\frac{x}{\lambda}\right)^{\theta-1}e^{-\left(x/\lambda\right)^\theta}, \quad x > 0\. $$ $Y=1/X$ then follows the inverse Weibull distribution with PDF $$ p_{Y}(y) = \left|\left(1/y\right)'_y\right| p_X\left(1/y\right) = \frac{\theta}{\lambda y^2} \left(\frac{1}{\lambda y} \right)^{\theta-1} e^{-\left(1/\lambda y\right)^\theta} = \frac{\theta}{\lambda^{-1}} \left(\frac{y}{\lambda^{-1}}\right)^{-\theta-1} e^{-\left(y/\lambda^{-1} \right)^{-\theta}}. $$ I.e. if negative shapes would be allowed, we could say $Y \propto \mathrm{Weibull}(-\theta,\lambda^{-1})$ (just $\theta$ in the front of $p_X(x)$ has to be replaced with $|\theta|$). CDF/mean/mode etc also require very minor adaptation for the negative shapes. I guess the same trick is possible for the whole generalized Gamma family, and e.g. inverse Gamma would then become its member.


title: Interpretation of Zero-One inflated Beta Regression with R (GAMLSS)

text: I am not that familiar with the interpretation of a beta regression with the r-package GAMLSS. Papers and package manuals didn`t help me. I modeled a Zero-One inflated Beta Regression. The y(usage)-variable consists of proportion data [0,1]. ![y=usage](http://i.stack.imgur.com/6LBE2.jpg) I had to specify four parameters for the beta distribution (mu, sigma, nu, tau). ![enter image description here](http://i.stack.imgur.com/rxBlN.jpg) The following data shows the summary ![enter image description here](http://i.stack.imgur.com/ZDqVU.jpg) * * * **Question** My intention is to show the relation between the y and x variables. Therefore i have got one question: 1.) How can i guess by reference to the summary if more or less of x1, x2, x3 and x4 influences y positively? (I guess mu just defines the possibility that y is somewhere between one and zero, nu: y is 0 and tau= y is 1)


title: Should I use coefficient correlation or multiple linear regression to build this scoring model?

text: Hi I am a junior Data analyst at a non-profit and I need some help with building a 'donor scoring model' to help us identify potential 'good' donors. I am attempting to build a donor scoring model by identifying potential variables for donors in our donor database and allocating weighted scores for each variable so that, the sum of the weighted scores can give us some indication of how 'good' a donor is (the way I am measuring how 'good' is time adjusted donation total $:=$ total amount donated/number of months as a donor). For example, I have identified 4 factors that may have a correlation to time adjusted donation total: birthdate, email, phone number and has Children. If we had this information on a donor, I gave them a 1 and 0 if we didn't have this information. What I did initially was calculate the correlation coefficient of each variable to the time adjusted donation total, and then adjusted the weighting for each variable based on the correlation coefficient. For example, the 'Has Children' factor had a correlation coefficient of 0.3, so I gave a score of 30, whereas birthdate had a correlation coefficient of 0.02, so I gave that a score of 2. Then adding these scores up for each variable, I calculated the correlation coefficient for the total score. I then manually adjusted the scores to try and yield the highest correlation coefficient for the total score. So my question is this: **Is this a reasonable approach to building a simple scoring system?** Is there a better way to 'optimise' the way I adjust the scores to get the best overall correlation coefficient? **Or should I use multiple regression to use the variables to predict the time adjusted donation total?** (I haven't used this, as my variables don't have a variety of values, only a 1 or 0 for whether or not we have that information for a donor). Any insights or direction on where I can learn more would be greatly appreciated.


title: Dealing with non-normal distribution in "big" datasets, when do we throw out the CLT?

text: Apologies from the go as this question comes from an absolute newbie and will definitely not satisfy a lot of the detail required. Hence, your guidance in providing you the right information to allow for the adequate answering of my question may be inevitable. A brief summary of what's nagging at my mind currently. I have received a dataset from a colleague, which contains ~5000 patients and we're trying to determine how an initial treatment affects various outcome variables. Some of these are categorical (with 3 levels), some are continuous. The categorical variable is currently the least of my concerns, because as far as I can tell, I would have to perform multivariate regression, perhaps using MANOVA, to see how this dependent variable is affected by the other independent variables. Now, what really bothers me is the distribution of the continuous data. This data consists of classic "count" data (number of times the doctor was visited following treatment) but then contains "time to heal". The "count" data would assume a Poisson distribution, but in this case var>mean ... so perhaps a quasi-Poisson would be more adequate. The distribution of that "count" data is depicted in figure 1. ![Figure 1](http://i.stack.imgur.com/2njYo.png) Now the time to heal is what's really bugging me. Some people would argue, that this would require linear regression modelling, but what I can't get my head around, is that time to heal, cannot be zero or take on negative numbers. And if one looks at the distribution of the data, it looks a lot more Poisson than anything else (Figure 2). So would this in fact require Poisson regression? Also biologically, time to heal is probably not completely linear, but more along what we see in cell culture growth, that it has a more exponential pattern, with a linear phase in it. ![Figure 2](http://i.stack.imgur.com/ALZ2d.png) OR - and this brings me to the question in my title - because of the size of the dataset, which is fairly large, could one say that the central limit theorem holds true and I don't have to worry about the data-distribution? Clearly my plots tell me I do, no?!? If I log-transform my data, things start to look quite ok (Figure 3)...but is a linear model adequate as the values cannot be negative or zero? ![Figure 3](http://i.stack.imgur.com/jkUgY.png) I've plotted the residuals of a linear regression model for the time to heal data in Figure 4, which shows that the tails don't really hold true to linearity. Then I plotted the residuals of the regression using a Poisson regression (Figure 5). Then the residuals of a linear regression model using log-transformation (Figure 6). ![Figure 3 - Linear Model](http://i.stack.imgur.com/eGQTO.png) ![Figure 4 - Poisson model](http://i.stack.imgur.com/HSpd1.png) ![Figure 6 - Log transformed linear model](http://i.stack.imgur.com/A0OQP.png) Just by looking at what you have here, what more information would you require? Am I going down a completely wrong path here? And if we assume an exponential curve for time to heal, what do I have to do to correct for this in my modelling process? And just because of the rule of numbers, could I spare myself all this thinking and simply assume that the CLT holds true? I am very new to this, so any guidance would be more than welcome. Hope not to have wasted anyone's time and apologize if I'm not posting crucial information here. Happy to learn, what more is required.


title: Identifying the number of conditions by which variance components will be divided in Generalizability Theory using QME R package

text: I'm using `QME` `R` `Package` (can be downloaded from here) to calculate Generalizability Coefficient. * * * ## MWE * * *               ########################################     #An example in which all facets are considered random.          #Load Brennan's (2001) data set number 4.     data(Brennan.4)     head(Brennan.4)     sapply(Brennan.4, class)          #Estimate generalizability study variance components.     formula.Brennan.4 <- as.formula("Score ~ (1 | Person) + (1 | Item) + (1 | Rater:Item) + (1 | Person:Item)")     gVar.out <- gVar(data.g = Brennan.4, formula.g = formula.Brennan.4)     gVar.out          #Complete G study.     gVar.out$N <- c(3, 12, 1, 3, 12)     gCoef(gVar.out)      I'm confused about `N`. How to determine `N` here? Any help in this regard will be highly appreciated. Thanks


title: Cross Validated Squared Correlation Coefficient in SVM regression

text: I am applying SVM regression on a data set on "Karthikeyan Melting Point Dataset" present here. I am using LIBSVM for it. After applying regression LIBSVM returns the cross validated squared correlation coefficient and RMS value (I am only interested in the 1st value). Is this value the same as the coefficient of determination (COD) or _R_ squared value? If yes, in a few cases I am getting values greater than 1, but as I have read, _R_ ² or COD is between 0 to 1. Please note this value is returned by LIBSVM.


title: Impractically long running time PCA command in R RStudio

text: I am using R in RStudio on OS X ver. 10.9.2 on 1.7 GHz Intel Core i7 with 8 GB RAM. I am trying to run a PCA command (`prcomp`) and plots on a dataset with approximately 200,000 rows and 8 columns. The prompt doesn't return, and most plots don't display (they do work on small ~300 rows datasets). Am I using the wrong tools for the job (R, hardware)? Am I doing something else wrong? Any good suggestions?


title: Why doesn't regression results change after bootstrap?

text: I learned bootstrap is used to treat non-normality of residual and it basically does resampling. I did bootstrapping on Stata and compared the result with normal regression.                reg salary salbegin educ, vce(bootstrap, reps(1000))     ------------------------------------------------------------------------------                  |   Observed   Bootstrap                         Normal-based           salary |      Coef.   Std. Err.      z    P>|z|     [95% Conf. Interval]     -------------+----------------------------------------------------------------         salbegin |   1.672631   .0962783    17.37   0.000     1.483928    1.861333             educ |    1020.39   164.8032     6.19   0.000     697.3818    1343.398            _cons |  -7808.714   1582.694    -4.93   0.000    -10910.74   -4706.692     ------------------------------------------------------------------------------                reg salary salbegin educ     ------------------------------------------------------------------------------           salary |      Coef.   Std. Err.      t    P>|t|     [95% Conf. Interval]     -------------+----------------------------------------------------------------         salbegin |   1.672631    .058847    28.42   0.000     1.556995    1.788266             educ |    1020.39   160.5504     6.36   0.000     704.9064    1335.874            _cons |  -7808.714    1753.86    -4.45   0.000    -11255.07   -4362.355     ------------------------------------------------------------------------------      Surprisingly, coefficient and p-value don't seem to change. If bootstrap does resampling, why doesn't regression result change? What is bootstrap actually doing? Does it just correct standard error?


title: Definition of the deterministic annealing method

text: iI have ran into a shape matching problem and one term which I read about is `deterministic annealing`. I learnt that it would help to convert discrete problems, e.g. **travelling salesman problem** to continuous problems to assist of sticking in the local minimas. I didn't quite understand what is exactly the definition of this method. Can somebody explain it?


title: Software for assisted decision tree construction

text: I have a large dataset based on several thousand surveys consisting of hundreds of questions each. I would like to form a classification tree semi-automatically as follows. Each node of the tree can divide cases based on a **single question only** (which may have multiple possible responses). I would like to build the tree manually, but I would like the software to suggest at each stage which questions might best divide the data[1], and allow me to select the one to use. (For example I may decide that dividing cases according to sex of respondent isn't interesting, at least at the top level, so will pick another question for the root node of the tree). Two questions * Is there any (preferably free) software out there that does this? * Is there any software that can do this without requiring a target variable? (Which most decision tree algorithms seem to require) [1] I am aware 'best' is ambiguous here. Feel free to suggest different ways I might be interested in deciding on 'best', though I am likely to be guided by the ready availability of software that implements whichever method. **UPDATE** I modified title and question based on initial responses.


title: Least Squares Matrix Notation Correlation

text: Consider the model $Y = b_0 + b_1 x_1 + b_2 x_2 + \epsilon$, where the columns $x_1$ and $x_2$ of the design matrix have mean 0 and length 1. That is $x_i' x_i = 1$, and $x_i' J = 0$, for $i = 1, 2$, where $J$ is a vector consisting entirely of ones. Let $p$ be the correlation between $x_1$ and $x_2$. Show that, $$\begin{align}X'X & = \begin{pmatrix} n & 0 & 0 \\\ 0 & 1 & p \\\ 0 & p & 1 \end{pmatrix} \end{align} . $$ **I can easily show that numerically with a $3 \times 3$ matrix of: $$\begin{align}X & = \begin{pmatrix} 1 & -\sqrt{.5} & -\sqrt{.5} \\\ 1 & 0 & 0 \\\ 1 & \sqrt{.5} & \sqrt{.5} \end{pmatrix} \end{align} . $$ However, I'm not sure how to generalize this beyond a $3 \times 3$ (for example, a $10 \times 3$). Any guidance would be appreciated. **


title: How to apply Kalman filter to one dimensional data?

text: I asked a question on StackOverflow for which I was suggested to use Kalman Filter. The question is as follows: http://stackoverflow.com/questions/5726358/what-class-of-algorithms-reduce- margin-of-error-in-continuous-stream-of-input/5728373#5728373 > A machine is taking measurements and giving me discrete numbers continuously > like so: > > 1 2 5 7 8 10 11 12 13 14 18 > > Let us say these measurements can be off by 2 points and a measurement is > generated every 5 seconds. I want to ignore the measurements that may > potentially be same > > Like continuous 2 and 3 could be same because margin of error is 2 so how do > I partition the data such that I get only distinct measurements but I would > also want to handle the situation in which the measurements are continuously > increasing like so: > > 1 2 3 4 5 6 7 8 9 10 > > In this case if we keep ignoring the consecutive numbers with difference of > less than 2 then we might lose actual measurements. Now how do I apply Kalman Filter to solve this? All examples I see take multiple error estimations while I know a single thing that each measurement can be off by a value Q thats it and all the examples also work on multi dimensional vectors, that too multiple vectors.


title: Space-efficient clustering

text: Most clustering algorithms I've seen start with creating a each-to-each distances among all points, which becomes problematic on larger datasets. Is there one that doesn't do it? Or does it in some sort of partial/approximate/staggered approach? Which clustering algorithm/implemention takes less than O(n^2) space? Is there a list of algorithms and their Time and Space requirements somewhere?


title: Manipulating and searching data-frames

text: Given a data-frame:               d1 <-c("A","B","C","A")     d2 <-c("A","V","C","F")     d3 <-c("B","V","E","F")     d4 <-c("A","B","C","A")     data.frame(d1,d2,d3,d4)       d1 d2 d3 d4     1  A  A  D  A     2  B  V  B  B     3  C  C  C  C     4  A  F  A  A      Also given that each row may have a unique pattern such that the occurrence of the values A,D,A (first row) represents a unique pattern assigned to a class 1 and F,A,A last row also represents a unique pattern assigned a class 4. I would like to manipulate the data-frame to search for rows that contain such 'unique patterns' and return a new column that classifies them such that, 0 represents rows that do not have any of the patterns. The pattern has to occur exactly as indicated.                  d1 d2 d3 d4  class     1  A  A  D  A     1     2  B  V  B  B     0     3  C  C  C  C     0     4  A  F  A  A     4      I tried to use a select statement with a concat qualifier using package sqldf, but it does not provide a useful approach. I would appreciate ideas on how to perform the search or if there are relevant packages to perform this type of search. Thank you


title: Given two responses for two groups, how to decide what to test on response or profit?

text: For many years I have been conducting t-tests on response to mailing activity. Recently I was challenged that we should infact be conducting tests on profit rather than response. So, let me put this in context. If you have two groups of customers of sample size 10,000 each that you were mailing two different offers. One receives free delivery (group a) on goods and one receives 25% off their next purchase (group b). And assume the two groups responded as follows.   * Group a - 10% of the 10,000 customers mailed responded    * Group b - 15% of the 10,000 customers mailed responded I would conduct a t-test on the response rates (10% and 15%), which is correct. I have been asked to do a test on the profitability, which i am reluctant to do. Again to put this into context:   * Group a - generated a profit of £1 per customer mailed (i.e. £10,000)   * Group b - generated a profit of £1.50 per customer mailed (i.e. £15,000) I have been asked to do a test on the £1 and £1.50 profit, but my concern is that there are many variable factors that contribute to the final profit figures, eg, margins, cost of the mailing the customer received etc, that are not really considered when setting the groups up. Am I correct in thinking that for this particular scenario that testing on profit is not the right thing to do, instead test on response as we are currently doing and then factor in all the other variable costs when making a final commercial decision?


title: How to apply Kalman filter to geographical coordinates?

text: I have a constant stream of geographical coordinates like so: (lat1, long1),(lat2, long2),(lat3, long3) I know that device taking measurement can be off by 100 meters. How do I filter out noise and get the actual measurements? I mean if two positions are with in 100 meters of range then they can be considered as a single position. How can I use Kalman filter to solve this problem?


title: Solving intercept / equivalence point / crossing point of where two Weibull regressions meet

text: I'm looking to solve the point where 2 opposite **Weibull** functions meet. I'm using the drc package, with a type 2 Weibull having 2 parameters in `R`. I've fit both lines to my data, and now I have the _e_ and _b_ intercepts for the functions defining each line. I would like to now solve for where the x values are equivalent. The **Weibull (type2) 2 parameter** function is: $$ f(x)= 1-\exp\left[-\exp\left(b\times (\log(x) - \log(e)\right))\right] $$ exponential decay line:                          Estimate Std. Error   t-value p-value      b:(Intercept) -2.226194   0.225339 -9.879323       0      e:(Intercept)  1.209326   0.072042 16.786444       0      exponential increase line:                              Estimate Std. Error   t-value p-value     b:(Intercept)  1.616248   0.145047 11.142956       0     e:(Intercept)  1.837511   0.072107 25.482970       0      Ideally I would also like to get a confidence interval for this point. I have this run in R so if I could input a suggested code that would be great. This is for a paper I'm writing, and I would add the solver of this question to this publication. In the paper, I'm using the Weibull function to model a PCR reaction. Here is the code for the exponential increasing instance:                 # first import the raw data:       # I pull it from Excel          > CCrelative <- read.xlsx('CC relative fold increase.xlsx', 1)     > colnames(CCrelative) <- c("MolOffTarget", "x1", "x2")     > print(CCrelative)        #here it is after the import                MolOffTarget    x1    x2         1      50001.0 8.474 8.372         2       5001.0 7.795 7.617         3        501.0 7.090 7.291         4         56.0 6.258 4.803         5          6.0 2.093 1.890         6          1.5 0.911 0.679         7          1.1 0.480 0.508                 #now set the minimum of each column to 0 by subtracting the lowest value from each     > CClog0 <- transform(CCrelative, x1 = (x1-(CCrelative[7,2])), x2=(x2-(CCrelative[7,3])))           #now normalize to the maximum value      > CClogT <- transform(CClog0, x1 = (x1/(CClog0[1,2])), x2=(x2/(CClog0[1,3])))           # now make the first column log10 scale     > CClogT <- transform(CClogT, MolOffTarget = log10(MolOffTarget))           #now merge columns x1 and x2 into a new data frame     > CClogT2 <- data.frame(rep(CClogT$MolOffTarget, 2), c(CClogT$x1,CClogT$x2))     > colnames(CClogT2) <- c("MolOffTarget", "FAM")     > CClogT2.W2.2 <- drm(FAM ~ MolOffTarget, data = CClogT2, fct = W2.2())      Here is the code for the exponential decay:               > CCvicRelative <- read.xlsx('CCvicRelative fold increase.xlsx', 1)     > colnames(CCvicRelative) <- c("MolOffTarget", "x1", "x2")     > print(CCvicRelative)            MolOffTarget     x1     x2     1      50001.0 -0.717 -0.706     2       5001.0 -0.565 -0.567     3        501.0 -0.360 -0.349     4         56.0 -0.001  0.584     5          6.0  1.568  1.582     6          1.5  1.767  1.822     7          1.1  1.802  1.844           #now set the minimum of each column to 0 by subtracting the lowest value from each     > CCvic0 <- transform(CCvicRelative, x1 = (x1-(CCvicRelative[1,2])), x2=(x2-(CCvicRelative[1,3])))           #now normalize to the maximum value      > CCvicT <- transform(CCvic0, x1 = (x1/(CCvic0[7,2])), x2=(x2/(CCvic0[7,3])))           # now make the first column log10 scale     > CCvicTlog <- transform(CCvicT, MolOffTarget = log10(MolOffTarget))           #now merge columns x1 and x2 into a new data frame     > CCvicT2 <- data.frame(rep(CCvicTlog$MolOffTarget, 2), c(CCvicTlog$x1,CCvicTlog$x2))     > colnames(CCvicT2) <- c("MolOffTarget", "VIC")     > CCvicT2.W2.2 <- drm(VIC ~ MolOffTarget, data = CCvicT2, fct = W2.2())      Apologies it's coming up a bit funny on this HTML. It took me a bit to get used to this package. I was using Bioconductor packages before. I made a mistake with the Weibull function earlier. The correct one is now present.


title: Leave-one-out cross-validation for Gaussian graphical models inverse covariance estimation

text: I want to calculate the test log likelihood of an estimated Gaussian graphical model inverse covariance matrix in leave-one-out cross validation test in R. So, what I want to calculate is $log(det(\theta))-trace(S\theta)$ where $S$ is the sample covariance matrix of the CV test data and $\theta$ is the inverse covariance matrix estimated using CV training data. The problem is that, since test data contains just one sample (dimension of $1\times p$ where $p$ is the number of variables), the empirical covariance matrix in this formulation ($S$) is not $p \times p$ but it is $1 \times 1$ (calculated using $var()$ or $cov()$ functions in R). But $\theta$ is estimated using $n-1$ samples of training data, so it is $p\times p$, so I am confused about how to calculate $trace(S\theta)$ since the dimensions of $S$ and $\theta$ are not suitable to multiply. If the sample size of test data is bigger than 1 (other cross validations than leave-one-out), there is no such a problem. So, how to calculate GGM test log likelihood for leave-one-out CV? Thanks.


title: Have Narrow Distribution - Need Standard Distribution

text: I am building a model in R using GLM based on this predictor variable. As you can see, the data is concentrated in the center of the distribution, and then falls off sharply. Is there a transformation to make this distribution wider? Another option I'm considering is to weight the data points on the ends, using the weights argument in GLM. Is there any reason why that would be a bad idea? ![enter image description here](http://i.stack.imgur.com/pQhxZ.png)


title: Longitudinal data format in WinBUGS or OpenBUGS, different number of items per time point

text: I am working with longitudinal data, where there are a different number of "items" per time point. I am working with the `MCMCdynamicIRT1d` function in `MCMCpack` in R to estimate a dynamic IRT model with binary items. However, now I am interested in modeling ordinal items within the same model. I found an example of a dynamic IRT model with ordinal data in BUGS and the code is based on a CSV data file looks like this:               CTRY            YEAR    Item1 Item2 Item3 Item4     Afghanistan     1981        0     1     0     0     Afghanistan     1982        2     0     1     2     ....     Albania         1981        0     2     0     0     Albania         1982        1     0     0     1     ....      The issue is that my data have a different number of "items" per timepoint. I am working with US Supreme Court cases, and the Court sees a different number of cases per year. So my data might look something like this:               Justice1 Justice2   year  time            0        0   1960     1            1        0   1960     1            0        1   1960     1            0        0   1961     2            0        0   1961     2            0        0   1961     2            0        0   1961     2            0        0   1962     3            0        0   1962     3      Any thoughts on how I can accommodate my data structure with a different number of "items" per time point? I appreciate any advice you can offer!!! ===== BUGS code I am trying to work with ====               model{     for(i in 1:n){# n is the number of obs     for(item in 1:4){     logit(Z[i, item, 1]) <- alpha[item, 1] - beta[item]*x[i]     logit(Z[i, item, 2]) <- alpha[item, 2] - beta[item]*x[i]     Pi[i, item, 1] <- Z[i, item, 1]     Pi[i, item, 2] <- Z[i, item, 2] - Z[i, item, 1]     Pi[i, item, 3] <- 1 - Z[i, item, 2]     y[i, item] ~ dcat(Pi[i, item, 1:3])     }          x[i] <- mu[country[i], year[i]]          }          sigma ~ dunif(0,1)     kappa <- pow(sigma, -1)          for(c in 1:n.country){     mu[c, 1] ~ dnorm(0, 1)          for(t in 2:n.year){ #n.year is number of years     mu[c, t] ~ dnorm(mu[c, t-1], kappa)          }     }          for(j in 1:4){     beta[j] ~ dgamma(4, 3)     alpha0[j, 1] ~ dnorm(0, .25)     alpha0[j, 2] ~ dnorm(0, .25)     alpha[j, 1:2] <- sort(alpha0[j, 1:2])     }     }


title: SEM variances of residuals fixed to 1?

text: I'm trying to perform structural equations with two second-order latent variables and five first order latent variables (since "a picture is worth a thousand words" I pasted my model below). To perform SEM correctly, I follow steps decribed in a book by Barbara Byrne(2010). My questions are:   1. Is it reasonable to fix variances of first order factor residuals (res_F1--res_F5) to 1?   2. Does the choice of reference class (for setting regression weights to 1) is so crucial? Because, sometimes I get some results or sometimes AMOS "says" that variance of res_F4 is negative, which I don't know how to interpret it. ![enter image description here](http://i.stack.imgur.com/qMShv.jpg)


title: How to pick the "sigma" when sigma clipping?

text: In astronomy (and I'm sure in other places as well) it is very common to use sigma-clipping as the outlier rejection scheme. The idea is that you regress the data, subtract the fit from the points, and then reject all points with residuals exceeding some sigma times either the standard deviation or median absolute deviation. My question: is there a good way to automatically select sigma based on the data? I have some points that have very high quality observations -- the measurement might be 15 magnitudes and the the error on measurement is something like 0.000001 -- and some with much lower quality observations -- error around 0.05. I'd think that one would want different values of sigma for these different cases. Any ideas? Heuristics are fine for my purposes.


title: Simple regression proof using general formula

text: I want to derive the least square estimation of the coefficients for $y=\beta_0+\beta_1*x_1+\varepsilon$ using the general formula. Can someone walk me through how to get from B to D in the image below? ![http://i.imgur.com/9zwvO.jpg](http://i.stack.imgur.com/hw2xa.jpg)


title: Convert continues number to integer number in optimization algorithms in MATLAB

text: I'm using a continuous optimization algorithm for optimizing neural network's number of neurons in first and second layers besides feature selection so I used this structure for converting continues output (between 0 and 1) to positive integer values in MATLAB:               Neuron_Layer_I=(min(floor(1+21*Neuron1),21))+3;       For example in above equation `Neuron1` is output of optimization algorithm so if we have `0` for `Neuron1`, number of neurons in first layer is `3` and if we have `1` for `Neuron1`, number of neurons in first layer is `24`. I have `[0 24]` range for second layer too. My questions :   * Is this equation is best equation for converting my values?   * My system tendency is toward using `24 neurons` in first layer and `24 neurons` for second layer. when I change (decrease or increase) ranges for number of neurons in `layer 1` and `layer 2` my system tendency is toward maximum values of that ranges. Why? I have a binary classification problem and using `patternnet` function in MATLAB. I'm using `Imperialistic Competitive Algorithm (ICA)` for optimization. Some times ago I have optimization of `early stopping (max iterations)` in my system besides this structure and the best early stopping max iteration was towards maximum optimization rage as we discoursed above, so I removed it. Now I have same problem with number of neurons. Besides above questions, What is your recommendation for improving my system. My system doesn't have a specific sample-size, In my testing cases it is between `200~700` . This is a financial problem. Thanks.


title: Linear model- Understanding performances on training and test sets

text: I have a small normalized data set, 30 observations and 18 Predictors. All are continuous and some variable are related. I ran linear regression on it using Weka. The model automatically dropped some collinear predictors reducing the number of predictors from 18 to 7. I then plotted model's prediction for target variable versus actual target variable data (training+test set). The correlation coefficient was 0.92. ![prediction vs actual target variable data](http://i.stack.imgur.com/bWiAo.png) The mean squared error on test was 12% if I split the data so that 10% (3 instances) of it was used as test set. Error increases as I increase the test points which is expected as the data set is small. QUESTIONS-   * A quick look at the figure below, seems to convey that the model should perform well on most of the test set points too as the points don't seem that scattered and the model is linear. Why does this graph looks so 'deceptionally' good, as just changing the 'random seed' in Weka, that is using a different set of 3 training set points using the same the model, drops the correlation coefficient to 0.33 and increases the Mean Squared Error to 82%.    * In some cases the corr was still high (+0.80) but mean error was +100%! I believe that a 'linear model' cannot possibly 'overfit' data since it is just a straight line trying to fit the data points (correct me if I am wrong). Overfitting is usually an issue when dealing with high order polynomials. It cannot underfit too as it is performing quite well on the training set for this case, as indicated by high corr. So what is it doing? (The figures below are not from the data above. For illustration purpose only) ![enter image description here](http://i.stack.imgur.com/EYE9J.png) Many thanks for any help!


title: Why would one want to carry out an orthogonalized (gram schmit) regression?

text: A regression in which all the IVs are replaced by their residuals after regressing them with the rest of the IVs. I have seen (for example in Elements of Statistical Learning) how doing this in a multivariate regression yields the multivariate coefficients from univariate regressions. But other than this, why would one benefit from using orthogonalized IVs rather than the regular ones? What are some reasons to be doing this in a multivariate regression?


title: Goodness of fit likelihood ratio test with zero values

text: I have a vector of observed frequencies that have zero values in some cells and a vector of expected frequencies generated by a model. I would like to do a likelihood ratio test rather than a chi squared test because I have different nested models and would like to use $2 \ln\cal{L}$ difference to evaluate the models. My question is how to deal with the zero cells? Is it sensible to remove the contribution from the 0 cell like in the R-code below? Or should I combine cells or do a continuity correction? Example:               observed <- c(142,34,18,15,12,136,6,7,2,3,2,1,1,0,2)     expected <- c(141.99,53.33,32.59,23.66,18.78,15.77,13.75,12.33,11.3,10.54,                    9.98,9.57,9.29,9.1,9.02)          2 * (sum(observed * log(observed/expected), na.rm=TRUE))          [1] 455.8816


title: Optimization problem

text: A friend of mine sells $k$ models of blenders. Some of the blenders are very simple and cheap, others are very sophisticated and more expensive. His data consists, for each month, of the prices of each blender (which are fixed by him), and the number of sold units for each model. To establish some notation, he knows for months $j=1,\dots,n$ the vectors $$ (p_{1j},\dots,p_{kj}) \qquad \textrm{and} \qquad (n_{1j},\dots,n_{kj}) \, , $$ where $p_{ij}$ is the price of blender model $i$ during month $j$, and $n_{ij}$ is the number of sold units of blender model $i$ during month $j$. Given the data, he wants to determine prices $(p^*_1,\dots,p^*_k)$ which maximize the value of his expected future sales. I have some ideas on how to start modeling this problem with some sort of Poisson regression, but I really don't want to reinvent the wheel. It would also be nice to prove that the desired maximum exists under certain conditions. Would someone please give me pointers to the literature of this kind of problem?


title: Is it compulsory for a linear regression analysis that a dependent as well as independent variable have equal variance?

text: The literature suggests that we need to have dataset that meets the condition of homoscedasticity. However, it seems that such a condition is not proper.


title: Cannot determine the P-value

text: I have a chi squared of 0.1185 with 1 degree of freedom. How can i determine the P-value?


title: Testing the variance part of a Generalized Linear Model out of sample

text: Suppose I have a response vector and an ANOVA design (for simplicity, assume it’s a one-way ANOVA with two treatments). A few Generalized Linear Models (Poisson, Negative Binomial, etc) are fitted to the data. This is done separately for each of K experimental units. Each unit has a distinct response vector, but the design matrix is the same across units. There are a few ways to decide what particular GLM (i.e. the distribution of response) to use for each unit. E.g. one can say that the choice is determined by AIC separately for each unit. Apparently, using such unit-specific GLM amounts to a higher number of effective parameters, so it may be the case that using the same GLM for all units works better overall. The problem is that there is no treatment effect in any of the units, so I can’t construct a ROC curve because it requires having both true positives and true negatives, and the former are absent here. What I can and will do is to check whether the test size is preserved (i.e. a “good” strategy should result in calling about 5% of the units when testing for the treatment effect with the cutoff p-value = 0.05). I am wondering whether it’s possible to perform a more direct, out-of-sample test. Note that the value of predicted response is equal to the corresponding cell mean regardless of what GLM is used. Therefore, I will have to compare not the accuracy of predicting the future response, but the accuracy of predicting the variance of response or some other statistic. E.g. suppose for each unit I have 20 observations per treatment. I use 10 observations per treatment to fit a few GLMs. Each GLM produces an estimate of response variance in each cell. Using the other half of the sample, I compute the observed variance in each cell and the discrepancy between the observed and predicted values. The discrepancy measures are summed up across all of the units. If you have seen something like that before, please provide suggestions and references. Thanks, James


title: Monte Carlo Optimisation

text: I have been doing a LOT of research about this but I can't figure out exactly how to do it. I need to find a problem that can be solved using Monte Carlo Optimisation (it is important that it is an optimisation not a simulation). I have came across - calculating the area of complicated shapes and the TSP problem. Does anybody know if there is any code that uses random number generation (Monte Carlo method) for any of these optimisations (or other optimisations that can be solved by Monte Carlo?) I need to analyze the situation given some code. Thanks


title: What are differences between the terms "time series analysis" and "longitudinal data analysis"

text: When talking about longitudinal data, we may refer to data collected over time from the same subject / study unit repeatedly, thus there are correlations for the observations within the same subject, i.e., within-subject similarity. When talking about time-series data, we also refer to the data collected over a series of time and it seems very similar to the longitudinal setting mentioned above. I'm wondering if someone can provide a clear clarification between these two terms, what is the relationship and what are the differences?


title: Being uniform in log of some parameter?

text: Alright, I am reading this book called "Bayesian Data Analysis" and it states this idea of being uniform in log of something which I don't have any clue what that might mean?? Example 1: Introduction of Dirichlet distribution as a conjugate for Multinomial distribution: $$p(\theta|\alpha) \propto \prod\limits_{j=1}^{k}{\theta_{j}^{\alpha_{j}-1}}$$ where the distribution is restricted to nonnegative $\theta_{j}$'s with $\sum\limits_{i=1}^{k}{\theta_j}=1$. ... Setting $\alpha_j=0$ for all $j$ results in an improper prior distribution that is uniform in the $log(\theta_j)$. Example 2: Or somewhere else in the book about normal distributions it mentions that considering $p(\sigma^2)\propto \frac{1}{\sigma^2}$ is makes a uniform distribution with regard to $\sigma$ on log scale. Could somebody please explain the situation?


title: Comparing file sizes day on day for outliers, what test to use?

text: I am quite new to statistical analysis, and have to crunch through several thousand folders/files each day to, at a first glance, see if the file size indicates an error (too large/small etc). The file sizes are not normally distributed, and some files have upward/downward trends over time. At the moment, since I am programming this using matlab and I am new to both statistical and matlab scripting, I am taking 2* the standard deviation of the data, but this does not always result in an accurate test. Does anyone have any recommendations? Attached are two examples of outliers (green) and my current testing range (red). ![trending file size](http://i.stack.imgur.com/X8ci3.png) ![outlier](http://i.stack.imgur.com/kaJWc.png) I thought a moving average may be a good idea, but since this will take time to program I want to consult the experts before proceeding and I am unsure how to implement a test to check for outliers in this case.


title: Concerns about the size of odds-ratio estimates in binary logistic regression model

text: I am creating some binary regressions and am wondering if I should be surprised by the size of the odds ratios. In this particular example, the dependent variable indicates whether a student can meet a specified target in a second language. The explanatory variables are age as a continuous covariate and a measure of literacy in the first language. I am using SPSS to calculate the models, and they seem to be working fine (I have a large dataset of around 40,000 cases and very few empty cells in the data, no warning signs etc.) I have chosen to have 0 as the reference category for both outcome and categorical explanatory variable. As I understand it, the estimates I am seeing below tell me that if the student demonstrates L1 reading at level 4 they are 372 times as likely to be able to demonstrate the required skill in the second language. The pattern of estimates is exactly as I would expect, my issue is with whether I should be concerned at the magnitude of the estimates, or whether this is usual in cases where there is clearly a strong connection between the two variables. Many thanks! KB               > Parameter          B         Exp(B)     >      > (Intercept)      -3.077         > Age                .075        1.078      > [L1 Reading=4]    5.921      372.784      > [L1 Reading=3]    5.330      206.438      > [L1 Reading=2]    4.288       72.821      > [L1 Reading=1]    3.025       20.594     > [L1 Reading=0]        0            0      Crosstabulation of L1 and DV: > >     +-------+-------+-------+------+------+-------+-------+ >     |       |     0 |     1 |    2 |    3 |     4 |  Total| >     | 0.00  | 17537 |  6487 |  784 |  100 | 44    | 24952 | >     | 1.00  |  1284 |  9891 | 4287 | 1576 | 1328  | 18366 | >     | Total | 18821 | 16378 | 5071 | 1676 | 1372  | 43318 | >     +-------+-------+-------+------+------+-------+-------+ >


title: Ordinal Logistic Model

text: I am new in research and have just started. My anticipated study is on ''determinants of community participation in planning for HIV and AIDS non medical interventions. Now the question is if the model is relevant to my study. The index scale for my dependent variable is 1 and 2 which means low and high participation in planning respectively. There are several factors which are thought to influence the participation positively of negatively.....i.e J=0,1....j. Pls kindly advise if it's relevant to use.


title: Good data example needed with covariate affected by treatments

text: I have looked at a lot of R datasets, postings in DASL, and elsewhere, and am not finding very many good examples of interesting datasets illustrating analysis of covariance for experimental data. There are numerous "toy" datasets with contrived data in stat textbooks. **I'd like to have an example where:**   * The data are real, with an interesting story   * There is at least one treatment factor and two covariates   * At least one covariate is affected by one or more of the treatment factors, and one is not affected by treatments.   * Experimental rather than observational, preferably ## Background My real goal is to find a good example to put in the vignette for my R package. But a larger goal is that people need to see good examples to illustrate some important concerns in covariance analysis. Consider the following made-up scenario (and please understand that my knowledge of agriculture is superficial at best).   * We do an experiment where fertilizers are randomized to plots, and a crop is planted. After a suitable growing period, we harvest the crop and measure some quality characteristic - that's the response variable. But we also record total rainfall during the growing period, and soil acidity at time of harvest -- and, of course, which fertilizer was used. Thus we have two covariates and a treatment. The usual way to analyze the resulting data would be to fit a linear model with the treatment as a factor, and additive effects for the covariates. Then to summarize the results, one computes "adjusted means" (AKA least-squares means), which are predictions from the model for each fertilizer, at the average rainfall and the3 average soil acidity. This puts everything on an equal footing, because then when we compare these results, we are holding rainfall and acidity constant. But this is probably the wrong thing to do -- because the fertilizer probably affects the soil acidity as well as the response. This makes the adjusted means misleading, because the treatment effect includes its effect on acidity. One way to handle this would be to take acidity out of the model, then the rainfall-adjusted means would provide a fair comparison. But if acidity is important, this fairness comes at great cost, in the increase in residual variation. There are ways to work around this by using an adjusted version of acidity in the model instead of its original values. The upcoming update to my R package **lsmeans** will make this downright easy. But I want to have a good example to illustrate it. I will be very thankful to, and will duly acknowledge, anyone who can point me to some good illustrative datasets.


title: Sample size question - ill posed question?

text: I'm working on a project in which we are performing clustering on high dimensional data (~1000 variables) and looking for subpopulations of observations that result from clustering. Think gene expression analysis. More specifically, I have paired samples in which I've calculated the ratio of each gene between two conditions to give a set of relative values. I then reduce the number of variables by identifying those that have a given magnitude of change for a minimum number of observations - for example, keep only the variables that have 10X change (up or down) for at least 25% of the observations. Clustering is then done on that reduced set to identify subpopulations that have different levels of expression. The idea is to find different subpopulations that may respond differently to various experimental conditions. The question of sample size has come up in this form: Assuming our experimental conditions have an effect on 10% of the overall population, how many samples do we need in order to identify real subpopulations. I've looked at power calculations for t-test in which I assume that the two populations are the different subpopulations that may emerge from clustering, and I assume alpha=.05 is the significance desired on group separation of those two. I don't know how to determine effect size, though. And, this just seems like a contrived way to get to an answer. Thoughts?


title: What do you call a measurement that's halfway between ordinal and interval?

text: Suppose you have a measurement $Y_i$. Now suppose that the differences between two measurements $\delta_{ij}=Y_i-Y_j$ are an ordinal measurement. From my understanding, $Y_i$ is more informative than an ordinal measure because the differences are meaningful. But it also is less informative than an interval scale because the intervals are not a ratio measurement. What do you call $Y_i$? An example that I was thinking of is a points system for a ladder in a sporting competition. Suppose the top three teams in some sport have $72,60,58$ points respectively at the end of a season. Then I would expect the first team to be "better" than the second team _by a larger amount_ than the second team is "better" than the third team as $(72-60)>(60-58)$. But I wouldn't say that they are $6$ times "better" ($6=\frac{72-60}{60-58}$) using this measurement - I would just say (using poor grammar) "more better".


title: Conditional Expectaction 3 variables

text: Suppose $X,Y$ and $Z$ are multivariate normal with means and full covariance matrix. The conditional expectation $E(X | Y)$ is well know. What is the conditional expectation of $E(X | Y,Z)$ if $Y$ and $Z$ (and $X$) are correlated? Standard textbooks only seem to cover the case when $Y$ and $Z$ are uncorrelated.


title: Choosing an appropriate minibatch size for stochastic gradient descent (SGD)

text: Is there any literature that examines the choice of minibatch size when performing stochastic gradient descent? In my experience, it seems to be an empirical choice, usually found via cross-validation or using varying rules of thumb. Is it a good idea to slowly increase the minibatch size as validation error decreases? What effects would this have on generalization error? Am I better- off using an extremely small minibatch and updating my model hundreds of thousands of times? Would I be better off with a balanced number somewhere between extremely small, and batch?   Should I scale the size of my minibatch with the size of the dataset, or the expected number of features within the dataset? I obviously have a lot of questions about implementing minibatch learning schemes. Unfortunately, most papers I read don't really specify how they chose this hyperparameter. I've had some success from authors such as Yann LeCun, especially from the Tricks of the Trade collection of papers. However, I still haven't seen these questions fully addressed. Does anyone have any recommendations for papers, or advice as to what criteria I can use to determine good minibatch sizes when trying to learn features?


title: Features for time series classification

text: I consider the problem of (multiclass) **classification** based on **time series** of variable length $T$, that is, to find a function $$f(X_T) = y \in [1..K]\\\ \text{for } X_T = (x_1, \dots, x_T)\\\ \text{with } x_t \in \mathbb{R}^d ~,$$ via a global representation of the time serie by a set of selected **features** $v_i$ of fixed size $D$ independent of $T$, $$\phi(X_T) = v_1, \dots, v_D \in \mathbb{R}~,$$ and then use standard classification methods on this feature set. _I'm **not** interested in forecasting, i.e. predicting $x_{T+1}$_. For example, we may analyse the way a person walks to predict the gender of the person. **What are the standard features that I may take into account ?** In example, we can obviously use the _mean_ and _variance_ of the serie (or higher order moments) and also look into the frequency domain, like the energy contained in some interval of the _Discrete Fourier Transform_ of the serie (or _Discrete Wavelet Transform_ ).


title: Definition of Conditional Probability with multiple conditions

text: Specifically, say I have two events, A and B, and some distribution parameters $ \theta $, and I'd like to look at $P(A | B,\theta)$. So, the simplest definition of conditional probability is, given some events A and B, then $P(A|B) = \frac{P(A \cap B)}{P(B)}$. So if there are multiple events to condition on, like I have above, could I say that $P(A | B,\theta) \stackrel{?}{=} \frac{P((A | \theta)\cap(B | \theta))}{P(B|\theta)}$ or am I looking at the in totally the wrong way? I tend to psych myself out when I deal with probability sometimes, I'm not really sure why.


title: Transform log posterior probabilities for Bayesian model averaging

text: I have four models, each one is run with MCMC on a large number phylogenetic trees. Each resulting model has an associated log posterior probability (sum of the log likelihood and log prior). Since I need to average the models weighting by their posteriors (bayesian model averaging) I need the original posterior probabilities. How can I standardize them from 0 to 1? Should I normalize among the four model for each tree? Also consider that my log posteriors are negative and very low (<-700), then when I calculate their exponential I get 0. Thanks in advance


title: PLS-DA with binary predictors in R (package mixOmics)

text: I am trying to analyse a dataset with at minimum 50 explanatory variables coded as 0 and 1 for presence/absence and a binary response variable (case/control). The goal is to see how the variables can predict the separation between case and control. As there are more variables than observations I applied a partial least square discriminant analysis (PLS-DA) using the package mixOmics in R. However, when I want to test the significance of the analysis with PLSDA.test (package RVAideMemoire) I get a lot of warnings : `1: In pls(X, ind.mat, ncomp = ncomp, mode = "regression", ... : Zero- or near-zero variance predictors. Reset predictors matrix to not near-zero variance predictors. See $nzv for problematic predictors.` I guess the problem with near-zero variance results from the 0/1 coding of the predictor variables. I tried to convert the variables to factors, but this doesn't help. Is there a different analysis more suitable? How can I deal with presence/ absence variables as predictors? Thanks for your answer in advance!!!


title: How to determine if my sample size is large enough

text: I currently have 146 participants in 4 groups as follows:   * group 1 = 46   * group 2 = 40   * group 3 = 29   * group 4 = 37 Each group completes the same 20 trials, but the groups have different characteristics (which is why they are separated). I'd like to know if I have enough participants per category to get statistically significant results or if I need to keep running the study and get more participants.


title: How do you constrain variogram estimates to be a percentage?

text: I am interpolating soil texture data with GeoR and ARCGIS. The ordinary kriging process is straight forward and GeoR facilitates cross validation. Has anyone else produced variogram parameters, nugget and/or sill, greater than the range of clay 0-100%? This seems to ruin Kriging variance results. The output variance is greater than the range of observations. Thank you for your time.


title: Find out how unlikely Y is than X -- when X is 50% more likely to cure cancer than Y

text: So 50% likely to me sounds like 150% percent that X will cure cancer, and 100% that Y will cure cancer. Other than this, I have NO IDEA without googling how to figure this out. Detailed explanation is good. **edit: 66% unlikely (or likely???) that Y will cure cancer when compared to X --** **now What are the basic rules that is needed to figure this out, along with other similar problems?**


title: When making inferences about group means, are credible Intervals sensitive to within-subject variance while confidence intervals are not?

text: This is a spin off of this question: How to compare two groups with multiple measurements for each individual with R? In the answers there (if I understood correctly) I learned that within-subject variance does not effect inferences made about group means and it is ok to simply take the averages of averages to calculate group mean, then calculate within-group variance and use that to perform significance tests. I would like to use a method where the larger the within subject variance the less sure I am about the group means or understand why it does not make sense to desire that. Here is a plot of the original data along with some simulated data that used the same subject means, but sampled the individual measurements for each subject from a normal distribution using those means and a small within- subject variance (sd=.1). As can be seen the group level confidence intervals (bottom row) are unaffected by this (at least the way I calculated them). ![enter image description here](http://i.stack.imgur.com/eiRJ9.png) I also used rjags to estimate the group means in three ways. 1) Use the raw original data 2) Use only the Subject means 3) Use the simulated data with small within-subject sd The results are below. Using this method we see that the 95% credible intervals are narrower in cases #2 and #3. This meets my intuition of what I would like to occur when making inferences about group means, but I am not sure if this is just some artifact of my model or a property of credible intervals. Note. To use rjags you need to first install JAGS from here: http://sourceforge.net/projects/mcmc-jags/files/ ![enter image description here](http://i.stack.imgur.com/H2rnX.png) The various code is below. The original data:               structure(c(1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,      1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,      1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2,      2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2,      2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3,      3, 3, 3, 3, 3, 3, 3, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 3, 3,      3, 3, 3, 3, 4, 4, 4, 4, 4, 4, 5, 5, 5, 5, 5, 5, 6, 6, 6, 6, 6,      6, 7, 7, 7, 7, 7, 7, 8, 8, 8, 8, 8, 8, 9, 9, 9, 9, 9, 9, 10,      10, 10, 10, 10, 10, 11, 11, 11, 11, 11, 11, 12, 12, 12, 12, 12,      12, 13, 13, 13, 13, 13, 13, 14, 14, 14, 14, 14, 14, 15, 15, 15,      15, 15, 15, 16, 16, 16, 16, 16, 16, 17, 17, 17, 17, 17, 17, 18,      18, 18, 18, 18, 18, 2, 0, 16, 2, 16, 2, 8, 10, 8, 6, 4, 4, 8,      22, 12, 24, 16, 8, 24, 22, 6, 10, 10, 14, 8, 18, 8, 14, 8, 20,      6, 16, 6, 6, 16, 4, 2, 14, 12, 10, 4, 10, 10, 8, 4, 10, 16, 16,      2, 8, 4, 0, 0, 2, 16, 10, 16, 12, 14, 12, 8, 10, 12, 8, 14, 8,      12, 20, 8, 14, 2, 4, 8, 16, 10, 14, 8, 14, 12, 8, 14, 4, 8, 8,      10, 4, 8, 20, 8, 12, 12, 22, 14, 12, 26, 32, 22, 10, 16, 26,      20, 12, 16, 20, 18, 8, 10, 26), .Dim = c(108L, 3L), .Dimnames = list(         NULL, c("Group", "Subject", "Value")))      Get subject Means and simulate the data with small within-subject variance:               #Get Subject Means     means<-aggregate(Value~Group+Subject, data=dat, FUN=mean)          #Initialize "dat2" dataframe     dat2<-dat          #Sample individual measurements for each subject     temp=NULL     for(i in 1:nrow(means)){       temp<-c(temp,rnorm(6,means[i,3], .1))     }          #Set Simulated values     dat2[,3]<-temp      The function to fit the JAGS model:                require(rjags)           #Jags fit function     jags.fit<-function(dat2){            #Create JAGS model       modelstring = "            model{       for(n in 1:Ndata){       y[n]~dnorm(mu[subj[n]],tau[subj[n]]) T(0, )       }            for(s in 1:Nsubj){       mu[s]~dnorm(muG,tauG) T(0, )       tau[s] ~ dgamma(5,5)       }                 muG~dnorm(10,.01) T(0, )       tauG~dgamma(1,1)            }       "       writeLines(modelstring,con="model.txt")          #############              #Format Data       Ndata = nrow(dat2)       subj = as.integer( factor( dat2$Subject ,                                      levels=unique(dat2$Subject ) ) )       Nsubj = length(unique(subj))       y = as.numeric(dat2$Value)            dataList = list(         Ndata = Ndata ,         Nsubj = Nsubj ,         subj = subj ,         y = y       )            #Nodes to monitor       parameters=c("muG","tauG","mu","tau")                 #MCMC Settings       adaptSteps = 1000                    burnInSteps = 1000                   nChains = 1                          numSavedSteps= nChains*10000                 thinSteps=20                             nPerChain = ceiling( ( numSavedSteps * thinSteps ) / nChains )                             #Create Model       jagsModel = jags.model( "model.txt" , data=dataList,                                n.chains=nChains , n.adapt=adaptSteps , quiet=FALSE )       # Burn-in:       cat( "Burning in the MCMC chain...\n" )       update( jagsModel , n.iter=burnInSteps )            # Getting DIC data:       load.module("dic")                 # The saved MCMC chain:       cat( "Sampling final MCMC chain...\n" )       codaSamples = coda.samples( jagsModel , variable.names=parameters ,                                    n.iter=nPerChain , thin=thinSteps )              mcmcChain = as.matrix( codaSamples )            result = list(codaSamples=codaSamples, mcmcChain=mcmcChain)          }      Fit the model to each group of each dataset:               #Fit to raw data     groupA<-jags.fit(dat[which(dat[,1]==1),])     groupB<-jags.fit(dat[which(dat[,1]==2),])     groupC<-jags.fit(dat[which(dat[,1]==3),])          #Fit to subject mean data     groupA2<-jags.fit(means[which(means[,1]==1),])     groupB2<-jags.fit(means[which(means[,1]==2),])     groupC2<-jags.fit(means[which(means[,1]==3),])          #Fit to simulated raw data (within-subject sd=.1)     groupA3<-jags.fit(dat2[which(dat2[,1]==1),])     groupB3<-jags.fit(dat2[which(dat2[,1]==2),])     groupC3<-jags.fit(dat2[which(dat2[,1]==3),])      Credible interval/highest density interval function:               #HDI Function     get.HDI<-function(sampleVec,credMass){        sortedPts = sort( sampleVec )       ciIdxInc = floor( credMass * length( sortedPts ) )       nCIs = length( sortedPts ) - ciIdxInc       ciWidth = rep( 0 , nCIs )       for ( i in 1:nCIs ) {         ciWidth[ i ] = sortedPts[ i + ciIdxInc ] - sortedPts[ i ]       }       HDImin = sortedPts[ which.min( ciWidth ) ]       HDImax = sortedPts[ which.min( ciWidth ) + ciIdxInc ]       HDIlim = c( HDImin , HDImax, credMass )       return( HDIlim )     }      First Plot:               layout(matrix(c(1,1,2,2,3,4),nrow=3,ncol=2, byrow=T))          boxplot(dat[,3]~dat[,2],      xlab="Subject", ylab="Value", ylim=c(0, 1.2*max(dat[,3])),     col=c(rep("Red",length(which(dat[,1]==unique(dat[,1])[1]))/6),     rep("Green",length(which(dat[,1]==unique(dat[,1])[2]))/6),     rep("Blue",length(which(dat[,1]==unique(dat[,1])[3]))/6)     ),     main="Original Data"     )     stripchart(dat[,3]~dat[,2], vert=T, add=T, pch=16)     legend("topleft", legend=c("Group A", "Group B", "Group C", "Individual Means +/- 95% CI"),     col=c("Red","Green","Blue", "Grey"), lwd=3, bty="n", pch=c(15),     pt.cex=c(rep(0.1,3),1),     ncol=3)          for(i in 1:length(unique(dat[,2]))){       m<-mean(examp[which(dat[,2]==unique(dat[,2])[i]),3])       ci<-t.test(dat[which(dat[,2]==unique(dat[,2])[i]),3])$conf.int[1:2]            points(i-.3,m, pch=15,cex=1.5, col="Grey")       segments(i-.3,                ci[1],i-.3,                ci[2], lwd=4, col="Grey"       )     }                    boxplot(dat2[,3]~dat2[,2],      xlab="Subject", ylab="Value", ylim=c(0, 1.2*max(dat2[,3])),     col=c(rep("Red",length(which(dat2[,1]==unique(dat2[,1])[1]))/6),     rep("Green",length(which(dat2[,1]==unique(dat2[,1])[2]))/6),     rep("Blue",length(which(dat2[,1]==unique(dat2[,1])[3]))/6)     ),     main=c("Simulated Data", "Same Subject Means but Within-Subject SD=.1")     )     stripchart(dat2[,3]~dat2[,2], vert=T, add=T, pch=16)     legend("topleft", legend=c("Group A", "Group B", "Group C", "Individual Means +/- 95% CI"),     col=c("Red","Green","Blue", "Grey"), lwd=3, bty="n", pch=c(15),     pt.cex=c(rep(0.1,3),1),     ncol=3)          for(i in 1:length(unique(dat2[,2]))){       m<-mean(examp[which(dat2[,2]==unique(dat2[,2])[i]),3])       ci<-t.test(dat2[which(dat2[,2]==unique(dat2[,2])[i]),3])$conf.int[1:2]            points(i-.3,m, pch=15,cex=1.5, col="Grey")       segments(i-.3,                ci[1],i-.3,                ci[2], lwd=4, col="Grey"       )     }               means<-aggregate(Value~Group+Subject, data=dat, FUN=mean)          boxplot(means[,3]~means[,1], col=c("Red","Green","Blue"),     ylim=c(0,1.2*max(means[,3])), ylab="Value", xlab="Group",     main="Original Data"     )     stripchart(means[,3]~means[,1], pch=16, vert=T, add=T)          for(i in 1:length(unique(means[,1]))){       m<-mean(means[which(means[,1]==unique(means[,1])[i]),3])       ci<-t.test(means[which(means[,1]==unique(means[,1])[i]),3])$conf.int[1:2]            points(i-.3,m, pch=15,cex=1.5, col="Grey")       segments(i-.3,                ci[1],i-.3,                ci[2], lwd=4, col="Grey"       )     }     legend("topleft", legend=c("Group Means +/- 95% CI"), bty="n", pch=15, lwd=3, col="Grey")               means2<-aggregate(Value~Group+Subject, data=dat2, FUN=mean)          boxplot(means2[,3]~means2[,1], col=c("Red","Green","Blue"),     ylim=c(0,1.2*max(means2[,3])), ylab="Value", xlab="Group",     main="Simulated Data Group Averages"     )     stripchart(means2[,3]~means2[,1], pch=16, vert=T, add=T)          for(i in 1:length(unique(means2[,1]))){       m<-mean(means[which(means2[,1]==unique(means2[,1])[i]),3])       ci<-t.test(means[which(means2[,1]==unique(means2[,1])[i]),3])$conf.int[1:2]            points(i-.3,m, pch=15,cex=1.5, col="Grey")       segments(i-.3,                ci[1],i-.3,                ci[2], lwd=4, col="Grey"       )     }     legend("topleft", legend=c("Group Means +/- 95% CI"), bty="n", pch=15, lwd=3,   col="Grey")      Second Plot:               layout(matrix(c(1,2,3,4,4,4,5,5,5,6,6,6),nrow=4,ncol=3, byrow=T))          #Plot priors     plot(seq(0,10,by=.01),dgamma(seq(0,10,by=.01),5,5), type="l", lwd=4,          xlab="Value", ylab="Density",          main="Prior on Within-Subject Precision"     )     plot(seq(0,10,by=.01),dgamma(seq(0,10,by=.01),1,1), type="l", lwd=4,          xlab="Value", ylab="Density",          main="Prior on Within-Group Precision"     )     plot(seq(0,300,by=.01),dnorm(seq(0,300,by=.01),10,100), type="l", lwd=4,          xlab="Value", ylab="Density",          main="Prior on Group Means"     )               #Set overall xmax value     x.max<-1.1*max(groupA$mcmcChain[,"muG"],groupB$mcmcChain[,"muG"],groupC$mcmcChain[,"muG"],                        groupA2$mcmcChain[,"muG"],groupB2$mcmcChain[,"muG"],groupC2$mcmcChain[,"muG"],                    groupA3$mcmcChain[,"muG"],groupB3$mcmcChain[,"muG"],groupC3$mcmcChain[,"muG"]     )               #Plot result for raw data     #Set ymax     y.max<-1.1*max(density(groupA$mcmcChain[,"muG"])$y,density(groupB$mcmcChain[,"muG"])$y,density(groupC$mcmcChain[,"muG"])$y)          plot(density(groupA$mcmcChain[,"muG"]),xlim=c(0,x.max),               ylim=c(-.1*y.max,y.max), lwd=3, col="Red",              main="Group Mean Estimates: Fit to Raw Data", xlab="Value"         )         lines(density(groupB$mcmcChain[,"muG"]), lwd=3, col="Green")     lines(density(groupC$mcmcChain[,"muG"]), lwd=3, col="Blue")          hdi<-get.HDI(groupA$mcmcChain[,"muG"], .95)     segments(hdi[1],-.033*y.max,hdi[2],-.033*y.max, lwd=3, col="Red")          hdi<-get.HDI(groupB$mcmcChain[,"muG"], .95)     segments(hdi[1],-.066*y.max,hdi[2],-.066*y.max, lwd=3, col="Green")          hdi<-get.HDI(groupC$mcmcChain[,"muG"], .95)     segments(hdi[1],-.099*y.max,hdi[2],-.099*y.max, lwd=3, col="Blue")          ####          #Plot result for mean data          #x.max<-1.1*max(groupA2$mcmcChain[,"muG"],groupB2$mcmcChain[,"muG"],groupC2$mcmcChain[,"muG"])         y.max<-1.1*max(density(groupA2$mcmcChain[,"muG"])$y,density(groupB2$mcmcChain[,"muG"])$y,density(groupC2$mcmcChain[,"muG"])$y)          plot(density(groupA2$mcmcChain[,"muG"]),xlim=c(0,x.max),               ylim=c(-.1*y.max,y.max), lwd=3, col="Red",              main="Group Mean Estimates: Fit to Subject Means", xlab="Value"         )         lines(density(groupB2$mcmcChain[,"muG"]), lwd=3, col="Green")     lines(density(groupC2$mcmcChain[,"muG"]), lwd=3, col="Blue")          hdi<-get.HDI(groupA2$mcmcChain[,"muG"], .95)     segments(hdi[1],-.033*y.max,hdi[2],-.033*y.max, lwd=3, col="Red")          hdi<-get.HDI(groupB2$mcmcChain[,"muG"], .95)     segments(hdi[1],-.066*y.max,hdi[2],-.066*y.max, lwd=3, col="Green")          hdi<-get.HDI(groupC2$mcmcChain[,"muG"], .95)     segments(hdi[1],-.099*y.max,hdi[2],-.099*y.max, lwd=3, col="Blue")                         ####     #Plot result for simulated data     #Set ymax     #x.max<-1.1*max(groupA3$mcmcChain[,"muG"],groupB3$mcmcChain[,"muG"],groupC3$mcmcChain[,"muG"])         y.max<-1.1*max(density(groupA3$mcmcChain[,"muG"])$y,density(groupB3$mcmcChain[,"muG"])$y,density(groupC3$mcmcChain[,"muG"])$y)          plot(density(groupA3$mcmcChain[,"muG"]),xlim=c(0,x.max),               ylim=c(-.1*y.max,y.max), lwd=3, col="Red",              main=c("Group Mean Estimates: Fit to Simulated data", "(Within-Subject SD=0.1)"), xlab="Value"         )         lines(density(groupB3$mcmcChain[,"muG"]), lwd=3, col="Green")     lines(density(groupC3$mcmcChain[,"muG"]), lwd=3, col="Blue")          hdi<-get.HDI(groupA3$mcmcChain[,"muG"], .95)     segments(hdi[1],-.033*y.max,hdi[2],-.033*y.max, lwd=3, col="Red")          hdi<-get.HDI(groupB3$mcmcChain[,"muG"], .95)     segments(hdi[1],-.066*y.max,hdi[2],-.066*y.max, lwd=3, col="Green")          hdi<-get.HDI(groupC3$mcmcChain[,"muG"], .95)     segments(hdi[1],-.099*y.max,hdi[2],-.099*y.max, lwd=3, col="Blue")      **EDIT with my personal version of the answer from @StéphaneLaurent** I used the model he described to sample from a normal distribution with mean=0, between subject variance =1 and within subject error/variance= 0.1,1,10,100. A subset of the confidence intervals are shown in the left panels while the distribution of their widths is shown by the corresponding right panels. This has convinced me that he is 100% correct. However, I am still confused by my example above but will follow this up with a new more focused question. ![enter image description here](http://i.stack.imgur.com/qVyVm.png) The code for the above simulation and charts:               dev.new()     par(mfrow=c(4,2))               num.sims<-10000     sigmaWvals<-c(.1,1,10,100)     muG<-0  #Grand Mean     sigma.between<-1  #Between Experiment sd          for(sigma.w in sigmaWvals){            sigma.within<-sigma.w #Within Experiment sd            out=matrix(nrow=num.sims,ncol=2)       for(i in 1:num.sims){              #Sample the three experiment means (mui, i=1:3)         mui<-rnorm(3,muG,sigma.between)              #Sample the three obersvations for each experiment (muij, i=1:3, j=1:3)         y1j<-rnorm(3,mui[1],sigma.within)         y2j<-rnorm(3,mui[2],sigma.within)         y3j<-rnorm(3,mui[3],sigma.within)                   #Put results in data frame         d<-as.data.frame(cbind(           c(rep(1,3),rep(2,3),rep(3,3)),           c(y1j, y2j, y3j )         ))         d[,1]<-as.factor(d[,1])              #Calculate means for each experiment         dmean<-aggregate(d[,2]~d[,1], data=d, FUN=mean)              #Add new confidence interval data to output         out[i,]<-t.test(dmean[,2])$conf.int[1:2]            }            #Calculate % of intervals that contained muG       cover<-matrix(nrow=nrow(out),ncol=1)       for(i in 1:nrow(out)){         cover[i]<-out[i,1]<muG & out[i,2]>muG       }                      sub<-floor(seq(1,nrow(out),length=100))       plot(out[sub,1], ylim=c(min(out[sub,1]),max(out[sub,2])),            xlab="Simulation #", ylab="Value", xaxt="n",            main=c(paste("# of Sims=",num.sims),                   paste("% CIs Including muG=",100*round(length(which(cover==T))/nrow(cover),3)))       )       axis(side=1, at=1:100, labels=sub)       points(out[sub,2])            cnt<-1       for(i in sub){         segments(cnt, out[i,1],cnt,out[i,2])         cnt<-cnt+1       }       abline(h=0, col="Red", lwd=3)            hist(out[,2]-out[,1], freq=F, xlab="Width of 95% CI",            main=c(paste("muG=", muG),                    paste("Sigma Between=",sigma.between),                    paste("Sigma Within=",sigma.within))       )          }


title: Rule of thumb - number of predictors - Poisson regression  rates

text: I am interested in estimating a Poisson regression for mortality rates, with number of deaths as the dependent variable and log(population size) as the offset. I have 50 observations (states). I am aware of rules of thumbs in varions regression models to figure out the maximum number of predictors. For example: Cox model: at least 5 events per predictor, logsitic regression: at least 5 observations where Y=1 per predictor _etc_... Is there a similar rule of thumb for Poisson models? MP


title: Iterative method to find Ridge Regression Parameter

text: I have seen a method whereby instead of trying to estimate the ridge parameter (k) directly from the data (using one of the many many ridge parameter estimators in the literature) you solve for it iteratively. The method is simple enough: You simply increase k (in suitably small steps) until the condition number is reduced blow 10. At first blush this seems like quite a nice solution to me but I've never seen a Ridge Regression paper/book that uses it. Update OK this is basically the method suggested by Marquardt "Generalized inverses, Ridge Regression, Biased Linear Estimation and Non-linear Estimation" the only difference being he used VIF's to measure the MC while this method uses the condition number. McDonald and Galrneau "A Monte-Carlo Evaluation of some Ridge-Type Estimators" note that this method is may not be appropriate for all data sets as it does not include the y values (observations). I still have not found a paper where the Marquardt method is tested against other estimators for the ridge parameter does anybody know of such a paper? Is this method theoretically sound though? Even if (as I suspect) it isn't does it really matter for the average practitioner who just want to produce more stable estimates of their Beta's (the weights in the regression) rather than having them "blow up" to grossly unrealistic values when they experience severe MC? Truly I would like to find a better method than this ideally with a solid theoretical underpinning but its hard to see from a practical view point it can be improved upon?


title: Advantages of doing "double lasso" or performing lasso twice?

text: I once heard a method of using the lasso twice (like a double-lasso) where you perform lasso on the original set of variables, say S1, obtain a sparse set called S2, and then perform lasso again on set S2 to obtain set S3. Is there a methodological term for this? Also, what are the advantages of doing lasso twice?


title: Why is Poisson regression used for count data?

text: I understand that for certain datasets such as voting it performs better. Why is Poisson regression used over ordinary linear regression or logistic regression? What is the mathematical motivation for it?


title: bias of peeking at AB test data and adjusting minimum detectable effect

text: Let's say we're running an A/B test on my website comparing blue button clicks (baseline) to green button clicks.   * I use http://www.evanmiller.org/ab-testing/sample-size.html to calculate my required number of subjects per branch with the following parameters:     * significance level of **5%**     * statistical power of **80%**     * an observed historical baseline conversion rate of **5%**     * a desired minimum detectable effect of **1%** (ie. conversions between 4% and 6% will be indistinguishable from the baseline) Using the calculator, I determine that we need **7,663 pageviews** to declare a result. Now let's say everyone gets impatient and decides to check in on the experiment after **only 900 pageviews**. The Game Plan: 1) If it turns out that the green button is **at least 3% better** than baseline, we will decide to conclude the experiment and declare the green button as the winner (a 3% MDE given the same other initial parameters requires only 894 pageviews according to the calculator). 2) If it turns out that the green button is **less than 3% better** than baseline after 900 pageviews, we will decide to keep the experiment running to it's full course of 7,663 pageviews and then make a conclusion at that time. Are we introducing bias with this Game Plan?


title: Marginal Likelihood Prior

text: I have a model with probability matrix for a distribution of $x$, $y=\\{0,1\\}$, $p(x,y|w)$ where $w=[w_1,w_2,w_3,w_4]$ $p(x=0,y=0)=w_1$, $p(x=0,y=1)=w_2$ $p(x=1,y=0)=w_3$, $p(x=1,y=1)=(1−w_1−w_2−w_3)$ I want to use Bayes Factors to compare it with another model, so I want to find the marginal likelihood of my model. Hence, I need a prior. **What distribution would you use for $p(w)$ in this case and what would be its form?**


title: A contrast suggested by the data

text: One approach to testing hypotheses suggested by the data is Scheffé's method. According to the Wikipedia article, the contrast $C = \sum_{i=1}^r c_i \mu_i$ (where $\sum_{i=1}^r c_i=0$) is estimated by the corresponding sample contrast $\hat C$, and the variance of the latter is estimated by $s^2_{\hat C} = \hat\sigma^2_e \sum_{i=1}^r c_i^2/n_i$ (where $n_i$ is the size of the sample from the $i$th population and $n=\sum_{i=1}^r n_i$), and with probability $1-\alpha$ all confidence intervals of the form $\hat C\pm s_C\sqrt{(r-1)F_{\alpha,r-1,n-r}}$ are valid (they all contain the parameters that they estimate). "All" means for all contrasts. So I'm looking at a variation on this situation. Instead of multiple comparisons, I've got this simple linear regression: $y_i = \alpha+\beta x_i+\text{error}_i$, with 22 data points, no two of which have a common $x$ value. 20 of them look as if fitting a line through them is reasonable. The other two are set apart from that line, and the line through those two points is _very_ close to the same in slope as is the line you get from fitting the model with only the other 20 points. So I've hazarded a quick guess as to how to emend Scheffé's method for this situation. (Maybe later I'll work through the details carefully and see if I can prove that my guess makes sense.)   * Introduce an indicator variable $w$ equal to $1$ in the 20 cases and equal to $2$ in the other two. Fit the model $y_i=\alpha + \gamma w_i+ \beta x_i+\text{error}_i$.   * Instead of $s^2_C$ use the usual mean squared error from this fit.    * Instead of contrasts among $Y_i$, use contrasts among the RESIDUALS from the model WITHOUT THE INDICATOR.   * Instead of $(r-1)F_{\alpha,r-1,N-r}$, use $(22-1)F_{\alpha,22-1,22-3}$ since there are $22-3$ degrees of freedom for error. Fools rush in . . . . Will this still look plausible after I've checked everything? PS: Note my comment below. And if this is the wrong thing to do, is there a right thing?


title: Comparing two classifiers on separate pairs of train and test datasets

text: In my problem I have 2 classifiers, C1 and C2. Both C1 and C2 are Naive Bayes classifiers but the difference between them is that they use different feature selection methods. Both classifiers are trained over a dataset of 10,000 instances (noisy-labeled), and tested over a different dataset of 1,000 instances (manually labeled), both datasets being balanced. Now, I have plotted accuracy of both classifiers on increasing number of instances, and I found by visual inspection that C2 has generally better accuracy and recall than C1 . I would like to know whether such difference is statistically significant or not to assess that C2 is better than C1. Previously, I used the same dataset for k-cross validation, got the mean and variation of the accuracies of both classifiers and computed student t-test on a specific amount of features. However, now I have 2 different datasets for training and testing. How could I perform a test in such situation? Should I get the mean of accuracies for all different feature amounts? Thanks in advance... **EDIT** Regarding the domain, I am dealing with sentiment analysis (SA), classifying text data in 3 classes : positive, negative and neutral. Regarding error cost, at this stage I suppose that all error costs are the same (although I understand that the cost of classifying negative as positive would higher than negative as neutral). Regarding the "practical significant difference" when dealing with SA I am assuming that an improvement of 1% is significant, since previous papers usually present such kind of improvements. I want to test the accuracy of C1 and C2 when trained over automatic-labeled data, and tested over manually-labeled data.


title: Parallelization of an iterative model

text: What I have: An iterative process based on the application of a very simple algebra to represent a rate, with total dependence among any iteration and its predecessor (one input parameter for (n)th iteration is the output parameter from the (n-1)th iteration). What I want: Some kind of shift on algebra allowing some level of independence among sequential iterations, this way allowing parallelization. What I know (or believe to) so far: There is extensive published literature on parallel iterative methods, but most of them are about classical methods like: http://machinelearning.wustl.edu/mlpapers/paper_files/NIPS2006_725.pdf But I'm failing to recognize which of them I could "use" to help me with my simple algebra. Question: Can someone help me pointing literature on the basics of algebra transformation aiming parallelization? (any thoughts on that matter will be very helpful)


title: Calculate confidence interval from point estimate and standard error

text: Im having a problem with calculating a confidence interval for an incidence rate ratio between two populations. Using PROC STDRATE from SAS, I've obtained the Incidence Rate Ratio (the difference in incidence rate between two time periods), the log rate ratio, the standard error and the z-value. I wonder (1) why on earth did SAS not provide confidence intervals for the rate ratio? (2) How can I myself calculate these? I've read: http://www.ncbi.nlm.nih.gov/pmc/articles/PMC1255808/ But this is not a mean, which is why the simple formula: SE = SD/√(sample size), probably isnt correct. Many thanks in advance for help!


title: Suitable number of classes for SVM in text categorization

text: I'm doing text categorization with R and SVM in the package e1071. I have around 30000 text files for training, and 10000 for test. The goal is to **_hierarchically_** categorize these files. For example, I have 13 categories in level 1, such as sports, literature, politics, etc, and in the second level, there are more than 300 categories. For instance, below sports category, there are sub-categories, like football, basketball, rugby, etc. There are **two strategies** to reach the categorization in level 2. **First** is to classify the files in first level (13 categories), and then recursively, classify the files among its own subcategories. **Second** strategy is more direct, i.e. we assign different labels to all categories (more than 300) in level 2, then we train the model with SVM. For the second strategy, although I have used SVD to doc-term matrix, reducing its dimension to 30,000 * 10. The **_svm_** function in package e1071 still breaks down, giving the error **_cannot allocate vector of size 12.4 Gb_**. So I'd like to ask you gurus, whether the large number of categories is a real problem for SVM? Specifically, in my case, which strategy will produce better results and is more feasible in practical ? Thanks very much for any helpful comments and solutions !


title: Understanding hypothesis testing, confidence intervals and frequencies

text: Sorry if this is a silly question but I seem to be mixing a few things together and would love if someone could straighten it all out for me. Please correct me! This is somewhat related to my post here: Proportion Composite score for DV in rANOVA- Need Help Interpreting Results I conducted a study where all participants answered 3 binary items in each of 4 conditions (2x2 design where each fixed, within subjects factor was crossed with the other). I am not including any between subject factors . The 3 items were averaged to form a proportion (yes, yes, I have seen the Dixon and Jaeger articles but I don't want to get into multi-level modeling so bear with me, please. If you do think I should do GEE- please tell me how since I have never used it before). My first hypothesis is that in the 1,1 and 1,2 conditions, people will be more likely to choose the less expensive product (A in 1,1 and B in 1,2). My second hypothesis is that in both 2,1 and 2,2 people will be more likely to choose B. I ran a repeated measures ANOVA and saw that the 2 main effects and the interaction are significant. I also got a table with the means and confidence intervals.   1. Instead of running an analysis, should I have compared each level separately to the criterion 0.5? Which test should I use since the data across conditions is correlated?   2. How do I explain main effects and levels in the interaction? Do I compare the levels to each other (i.e. are they different) and then look at the CI to see which direction each one is going (i.e. > or < than 0.5).    3. How do I interpret a CI that contains 0.5? Can I say that that level does not the preference for a product? I feel like I have a lot of things swimming around here. I appreciate all the help so thanks in advance.


title: Significance of variance components in Stata output

text: This might be trivial, but I'm used to HLM7 software output and now I'm switching to Stata (`xtmixed`). To give an example imagine I have students (level-1) nested within schools (level-2). Running an empty model, in HLM, I can easily see the variance component associated to each level, to see how much variation is at level-1 and how much is at level-2. Starting from that, I also calculate the intra- class correlation coefficient. Now I see the variance components, in the `xtmixed` output, these are reported as the standard deviation estimates of the intercept `sd(_cons)` and of the residuals `sd(Residual)`. How do I calculate the associated p-value to see if there is significance?