title: Given the below dependent variable description, should I chose either Ordered or Multinomial logistic regressions?

text: My dependent variable looks like a range of ranks. It actually might be considered that way. But the ranking is based on subjective non-quantifiable cutoffs. We assessed the behaviors of a group of people (confidentiality doesn't allow me to explain further) and we have categorized their behavior from non-cooperative (the score zero) to fully collaborative (the top score). There are some cutoffs defined within this range of zero and the top score, making it look like an order of ranks. The cutoffs really exist and are important behavioral factors, based on which we have determined if a subject should get a score of 1, 2, etc. But the problem is that the distance between 1 and 2 is not necessarily equal to the distance between 2 and 3. And by distance, I mean the extent of the improvement in the behavior pattern. That distance is not actually quantifiable. So I don't know should I use an ordered logit or a multinomial one? My preference would be the ordinal one since its interpretation would be easier! But I don't know if its assumptions allow me or not. The Wikipedia page has bad news apparently "The model only applies to data that meet the proportional odds assumption, that the relationship[clarification needed] between any two pairs of outcome groups is statistically the same. This means that the coefficients that describe the relationship between, say, the lowest versus all higher categories of the response variable are the same as those that describe the relationship between the next lowest category and all higher categories, etc. Because the relationship between all pairs of groups is the same, there is only one set of coefficients." But still I think I should ask, maybe there is some way around. I also much appreciate any other hints you might have for me.


title: Extended Cox-model with ancillary variable

text: I have survival data (coming from a labor market) looking like this                  gender  nation          profession year age event   survtime     1       F Austria  blue-/white-collar 1997  30  TRUE   707 days     2       M Austria          freelancer 2001  26  TRUE  1765 days     3       F Austria marginal employment 2006  57  TRUE    57 days     4       F Austria  blue-/white-collar 2011  33  TRUE   465 days     5       F Austria  blue-/white-collar 1997  22  TRUE    35 days     6       M Germany  blue-/white-collar 1996  38  TRUE   301 days     7       M Austria  blue-/white-collar 2007  24 FALSE  2251 days     8       F Germany  blue-/white-collar 2011  23  TRUE   172 days     9       M Austria  blue-/white-collar 2004  35  TRUE   249 days     10      M Austria          freelancer 2000  32  TRUE   366 days     11      M Germany  blue-/white-collar 2008  22  TRUE    64 days     12      F Austria  blue-/white-collar 1998  21  TRUE    56 days     13      M Austria  blue-/white-collar 2012  36 FALSE   200 days     14      F Austria  blue-/white-collar 2000  47  TRUE   841 days     15      F Austria  blue-/white-collar 2008  27  TRUE   117 days     16      F   other  blue-/white-collar 1971  24  TRUE 11324 days     17      M Austria          freelancer 2001  34  TRUE   456 days     18      M Austria  blue-/white-collar 2007  25  TRUE    31 days     19      F Austria  blue-/white-collar 1992  25  TRUE  2010 days     20      M   other  blue-/white-collar 2003  24  TRUE   347 days      (get it with `df <\- source("http://pastebin.com/raw.php?i=WjQUatKy")$value`) I would expect the `year`, which describes the year at which a labor contract started, to be an ancillary variable, meaning that I think that `survtime` is changing with the economic conditions (like when there is a recession, survtime is maybe shorter). How can I account for that in an extended Cox proportional hazard model? I use R and so far I have               df <- source("http://pastebin.com/raw.php?i=WjQUatKy")$value     s <- Surv(as.numeric(df$survtime), df$event)     cox <- coxph(s ~ gender + nation + strata(profession) + age + year, data=df)      which however doesn't take this into account.


title: Time series forecasting using autoregressive and linear terms in R

text: I have real daily market data which I'm looking at to create a model for forecasting. The model that I created (below) used autoregressive terms within a linear regression. I was sharing this with a colleague and he said "autoregressive variables are correlated with the other variables in multiple linear setting which creates multicollinarity problem, creating unreliable result." So I'm turning to the group for help. Here is the data and the analysis that I performed in R.               #Read in Data     MarketData = read.table('http://sharp-waterfall-3397.herokuapp.com/MarketCategories6.txt', header=TRUE,na.strings = "NA", sep=",")     MarketData$Month <- as.factor(MarketData$Month)     MarketData$Weekday <- as.factor(MarketData$Weekday)          str(MarketData)      I'm using the PERregress library to help with the autoregression using the `back()` function and to help with the residual diagnostics:               library(PERregress)     descStat(MarketData)      Subsetting the data for model building and prediction purposes:               Total = MarketData     MarketData = MarketData[1:268,]     attach(MarketData)      Here is a regression with everything that I can think of. Note you can have higher autoregressive terms but this will start to mask events since R will ignore the first several rows. Also just an FYI for some reason the residual analysis is breaking which I liked to look for points with undue leverage.               #Market1Category1 Regression for the markets with everything that I can think of it     Market1Category1Output=lm(Market1Category1 ~ Trend+Month2+Month3+Month4+                               Month5+Month6+Month7+Month8+Month9+Monday+                               Tuesday+Wednesday+Thursday+Friday+Saturday+                               Holiday1+Holiday2+Holiday3+Holiday4+                               Event1+Event2+Event3+Event4+Event5+Event6+Event7+                               Event8+Event9+Event10+Event11+Event12+Event13+                               Event14+Event15+Event16+Event17+Event18+Event19+                               Event20+Event21+Event22+Event23+Event24+Event25+                               Event26+Event27+Event28+                               back(Market1Category1)+back(Market1Category1, 2))     summary(Market1Category1Output)     acf(Market1Category1Output$residuals)     residualPlots(Market1Category1Output)      Here is the final equation. I'd like to say that I reduced the variables using partial f-test but I couldn't find an easy way to do this so if you know a function please let me know. Basically I looked at the change in adjusted $R^2$.               #Final regression equation      Market1Category1Output=lm(Market1Category1 ~ Month5+Month6+Month7+                               Tuesday+Wednesday+Thursday+Friday+Saturday+                               Holiday1+Holiday3+Event2+Event7+Event10+                               Event13+Event16+Event25+Event28+                               back(Market1Category1)+back(Market1Category1, 6))     summary(Market1Category1Output)     acf(Market1Category1Output$residuals)     residualPlots(Market1Category1Output)      Here is a plot of the actuals in green vs the predictions in blue but there's a problem:               plot(Time, Market1Category1, col='green')     points(Time, predict(Market1Category1Output, MarketData), col='blue', pch=20)      The issue is that predict will use the data values instead of it's predicted values for the autoregressive terms. In order to make it use predicted terms I created this loop. If you know a better way let me know.               dataSet2 <- Total     dataSet2[8:length(dataSet2$Time),"Market1Category1"] <- NA         for (i in (1:(length(dataSet2$Time)-7))) {       dataSet2[6+i+1,"Market1Category1"] <- 1       dataSet2[6+i+1,"Market1Category1"] <- predict(Market1Category1Output,                                                      dataSet2[0:6+i+1,])[6+1]      }      Here is the plot again with the results in blue using the predicted results for the autoregressive terms (with the exception of the first 7 since the model needs those to `predict`):               plot(Total$Time, Total$Market1Category1, col='green')     points(dataSet2$Time, dataSet2$Market1Category1, col='blue', pch=20)      So here are my questions in order of importance:   1. Does using autoregressive and linear terms violate any fundamental assumptions?   2. What issues can this cause and what analysis/steps should I do take to avoid these problems?   3. Is there a better approach to modeling this timeseries?   4. Is there a more efficient approach?   5. Given the residuals what steps would you take? Finally two questions which is just causing me more work than possibly necessary:   1. As you can see instead of using the factors for weekday and month I'm using separate conditional variables. I'm doing this because if I use the factor and a level turns out to be insignificant (e.g., Monday for days of the week). I can't remove it. Perhaps there's a way?   2. Is there a quick way to run a partial F-statistic to understand whether removing a variable makes sense?


title: Conditional Maximum Likelihood Estimation for ARMA(p,q)

text: I'm making a study on STARIMA (Space-Time Autoregressive Integrated Moving Average) Time Series. And one problem I have is the parameter estimation since there is no software yet that has a STARIMA function, so far from what I have searched. Now, I want to compute it manually using estimation procedure. And I prefer Conditional Maximum Likelihood Estimation as it is used in STARIMA. I want to start with the univariate time series, the ARMA(p,q). I wish to estimate first the parameters of ARMA(p,q), then I'll just apply the procedure in STARMA. According to William W.S. Wei author of _Time Series Analaysis: Univariate and Multivariate Methods_ that (Original text from his book): * * * For the general stationary ARMA(p,q) model$$Z_t^*=\phi_1Z_{t-1}^*+\cdots+\phi_pZ_{t-p}^*+a_t-\theta_1a_{t-1}-\cdots-\theta_qa_{t-q}\quad(7.2.1)$$where $Z_t^*=Z-\mu$ and $\{a_t\}$ are i.i.d. $N(0,\sigma_a^2)$ white noise, the joint probability density of $\mathbf{a}=(a_1,a_2,\cdots,a_n)'$ is given by$$P(\mathbf{a}|\phi,\mu,\theta,\sigma_a^2)=(2\pi\sigma_a^2)^{-n/2}exp\left[-\frac{1}{2\sigma_a^2}\sum_{t=1}^na_t^2\right].\quad(7.2.2)$$Rewriting $(7.2.1)$ as$$a_t=\theta_1a_{t-1}+\cdots+\theta_qa_{t-q}+Z_t^*-\phi Z_{t-1}^*-\cdots-\phi_pZ_{t-p}^*\quad(7.2.3)$$we can write down the likelihood function of the parameters $(\phi,\mu,\theta,\sigma_a^2)$. Let $\mathbf{Z}=(Z_1,Z_2,\cdots,Z_n)'$ and assume the initial conditions $\mathbf{Z_*}=(Z_{1-p},\cdots,Z_{-1},Z_0)'$ and $\mathbf{a_*}=(a_{1-q},\cdots,a_{-1},a_0)$. The conditional log-likelihood function$$\ln L_*(\phi,\mu,\theta,\sigma_a^2)=-\frac{n}{2}\ln 2\pi\sigma_a^2-\frac{S_*(\phi,\mu,\theta)}{2\sigma_a^2}\quad(7.2.4)$$where$$S_*(\phi,\mu,\theta)=\sum_{t=1}^na_t^2(\phi,\mu,\theta|\mathbf{Z_*},\mathbf{a_*},\mathbf{Z})\quad(7.2.5)$$is the conditional sum of squares function. The quantities of $\widehat{\phi},\widehat{\mu}$ and $\widehat{\theta}$, which maximize Equation $(7.2.4)$, are called the conditional maximum likelihood estimators. * * * Now, I don't have idea how to find the quantities $\widehat{\phi},\widehat{\mu}$ and $\widehat{\theta}$ that maximizes $(7.2.4)$. Is this the same with the usual Maximum Likelihood Estimation like for parameters of say Geometric Distribution? If so, how will I obtain the formula of each parameters ($\widehat{\phi},\widehat{\mu}$ and $\widehat{\theta}$)? I know its tedious to perform MLE, but just give me an idea how to set up the likelihood function, and I'll just perform the usual MLE for parameters of ARMA(p,q). Well, if you could give me other source for parameter estimation of STARIMA then that would be great. Thanks in advance!


title: Examples for Type I and Type II errors

text: I was checking on Type I (reject a true H$_{0}$) and Type II (fail to reject a false H$_{0}$) errors during hypothesis testing and got to to know the definitions. But I was looking for where and how do these errors occur in real time scenarios. It would be great if someone came up with an example and explained the process where these errors occur.


title: Is subtracting the mean from PCA necessary when using an SVD result that is feature scaled?

text: I've applied SVD to the original data matrix and eliminated insignificant columns and rows from U and V^T respectively using the Sigma values. I multiplied together my optimized U, Sigma, and V^T matrices to obtain a close approximation of the original matrix and normalized the matrix by applying a feature-scaling formula. I took the max and min of the entire matrix and for each value, applied this formula: (x - min) / (max - min) Now I used this feature-scaled matrix to get the PCA (getting covariance matrix, selecting significant eigenvectors, multiplying the original matrix by the transpose of the eigenmatrix). At which step do I subtract the mean and from what? For that matter, do I still have to subtract the mean since the data is already feature-scaled and SVD'd?


title: Distributions on the simplex with correlated components

text: I'm looking for some kind of distribution over the simplex in which components are correlated in an ordinal way. That is, if $p = (p_1, ..., p_J)$ is drawn from our distribution on the simplex, I would like $p_i$ to be positively correlated with its neighbours $p_{i + 1}$ and $p_{i - 1}$, say. A vanilla Dirichlet clearly cannot satisfy this requirement. One option I suppose is a mixture of Dirichlet distributions; for example, when $J = 4$ one could take $\mathcal D(1, 1, 0, 0) + \mathcal D(0, 1, 1, 0) + \mathcal D(0, 0, 1, 1)$ or something similar to induce correlation, but I'm wondering if there is something a little more natural. Another option I suppose is to take any distribution on $\{1, 2, ..., J\}$, say $f(j | \eta)$, put a distribution on $\eta$ take $p_j = f(j | \eta)$. So I could take, for example, $\eta \sim \mbox{Beta}(\alpha, \beta)$ and let $f(j | \eta) = {J \choose j} \eta^j (1 - \eta)^{J - j}$. At any rate, I'd like whatever I end up with to be as tractable as possible. The mixture of Dirichlet's is appealing because I could get some nice conditional conjugacy going for me, but it's not clear how to set things up. This question talks about the logistic normal distribution, but I don't know much about it; is it tractable for Bayesian inference? Of course, the components of a Dirichlet are already negatively correlated, and asking for "positive correlation" probably isn't totally coherent since if $p_i$ is large then it is, by nature, taking up most of the mass and hence forcing the probability of its neighbours to be small. Perhaps what I mean is that $p_i$ is positively correlated with $p_{i + 1} / \sum_{j \ne i} p_j$. Hopefully the question as stated is enough for people to know what I want and be able to help me.


title: How to set up a non-linear mixed effects model with random effects in R using nlme?

text: I have some data with predictor variables, A and B and response variable C. I have a grouping factor SITE.               df <- data.frame(A = c(0.4, 0.4, 0.2, 0.2, 0.2, 0.2, 0.2), B = c(0.3, 0.3, 0.1, 0.1, 0.1, 0.1, 0.1),     C = c(4.4, 4.3, 5.6, 4.7, 5.1, 4.5, 4.9), SITE = c("south","south","east","east", "east", "north", "north"))      The relationships between C and predictor variables are non-linear. I have not shown all the dataframe. I would like to fit a non-linear mixed effects model with random effects in `R`. I can use `nlme()` to fit this model. However, I am struggling to find any 'accessible' (I don't have a strong statistics background) documentation that goes through how I can set up the model. Could someone please assist me in sitting up my model? My 'attempt' to set up the model is:               m1 <- nlme(C ~ A + B, data = df, random = ~1 | SITE)      However, when I run the above, I get an error: argument "start" is missing, with no default.


title: Under what circumstances is an MA process or AR process appropriate?

text: I have a very basic question. Please let me know if this has been asked before, but in my defence I haven't seen it on Cross Validated. I understand that if a process depends on previous values of itself then it is an AR process. If it depends on previous errors than it is an MA process. My question is, when would one of either of these two situations occur? Does anyone have a solid example that illuminates the underlying issue regarding what it means for a process to be best modelled as MA vs AR? In my search for some sort of intuitive or dare I say philosophical insights into the difference between these two models, I keep running into material where an author makes a statement along the line of the above, and then immediately moves to discussing diagnostics: acf, pacf, etc. Although I have the gut feeling that there is no simple answer to my question, if anyone can either enlighten me or point me to a good resource, I'd be grateful.


title: What is "the" correct way of calculating Cp/Cpk values for Statistical Process Control?

text: There seems to be many different methods of calculating Cp/Cpk values in SPC, and I need help in determining which one to use at any given situation. Currently, we need to show Cp/Cpk values when drawing two different chart.   * X-MR Chart   * X-Bar R Chart Let's assume we have an LSL/USL of 3.0/3.2 respecitvely and we have the following sample: 3.55 2.41 3.61 2.71 2.91 3.42 3.52 2.92 2.98 3.02 3.14 3.25 2.98 3.33 3.31 So according to this website, here's how I would go about determining CPK value.   1. Calculate sample mean mean = 3.137333333   2. Calculate Standard Deviation σ = 0.335249729   3. Calculate Cpu, Cpk and take the smaller value as Cpk Cpu = (USL-sample mean)/3x standard deviation. Cpl = (LSL-sample mean)/3x standard deviation. Cpu = 0.062308444 Cpl = 0.460409628 **Cpk = 0.4604** This works fine and all for X-MR Chart, since we don't have any subgroups, but let's assume for X-Bar R Chart, we're looking at a subgroup of 3. Now according to this website, we're supposed to use a different formula to calculate standard deviation, namely: ![enter image description here](http://i.stack.imgur.com/BLsvV.png) This will lead to a different value of Cpk. My question is, which formula am I supposed to use to calculate Cp/Cpk value?


title: How to compare the results of two classifiers are statistically significant different?

text: I am applying kNN and SVM classifiers to my classification problem. Both of the classifiers get over 95% cross-validation accuracy (leave one out cross validation used). Not sure how to tell if the results of two classifiers are statistically significant different? I did a little bit search that Discussion1 Discussion2 these posts discussed this topic; however, I am still confused how they use hypothetical test for this purpose. For example, I have two cross-validation accuracies and how to form the 2 × 2 contingency tables of McNemar's test? If my performance measure is AUC of ROC curve, then how to perform the McNemar's test? Thanks a lot. A.


title: Addressing Non-response in a Convenience Sample

text: I am studying customer satisfaction in a large hierarchical organization. I plan to administer a voluntary survey to customers across the organization, and need to address non-response in my analysis. I know:   * How many customers there are (population size)   * Some demographic/market segment data for all customers   * The products/services they used   * The account manager and department that handles their account   * 50% response rate expected I want to draw conclusions about departments, account managers and products. It'd also be nice to be able to break things out by customer demographics, but that's secondary. The naive approach would assume that non-response is random and treat the data as a simple random sample of the population. My concern is that non-response is probably not completely random. So here are my questions:   * Is the 50% response rate sufficiently high that the naive approach isn't terrible?   * If not, is resampling or post-stratification a sufficient remedy?   * If not, does it make sense to try to model the sampling probability mechanism using the known customer demographic information?   * If none of these are sufficient, would I be better off designing a smaller stratified sample and aggressively pursuing non-responders? Related:   * Does SurveyMonkey ignore the fact that you get a non-random sample?   * Recommend references on survey sample weighting   * Ways to overcome small number of survey responses   * Are benchmarking and post-stratification equivalent? EDIT 4/22/2014: Reviewing @Steve Samuels response and additional independent research, I think what I'm dealing with is a census with non-response. The population (all current customers) is sufficiently small and well-known to serve as the sample frame. Thus by definition we have a census. It's actually less feasible to address a random sub-population than the whole population due to the platform we're using to deliver the questionnaire. My plan moving forward is to execute the census, then study the difference between the responding & non-responding populations in a variety of ways. Moving forward, follow-up efforts will be adjusted to address specific problems found in the non-response analysis. Some good resources for this problem include:   * Analysis of Customer Surveys: with applications using R (Kenett and Salini, 2011) has a good chapter on sample design and continues with a solid overview of survey analysis techniques, all with a focus on customer satisfaction surveys.   * The causes and consequences of response rates in surveys by the news media and government contractor survey research firms (Holbrook, Krosnick, and Pfent, 2007) does a meta-analysis of 114 phone surveys conducted by professional survey firms and "challenges the assumptions that response rates are a key indicator of survey data quality and that efforts to increase response rates will necessarily be worth the expense"   * Survey response rate levels and trends in organizational research (Bruch and Holtom, 2008) provides benchmarks for response rates in organizational surveys (50% is typical) and guidelines for responsibly reporting response statistics.


title: Is there a way to specify a lme model with more than one within-subjects factor?

text:
## The data Suppose we have a dataset `d` with two between-subject factors (i.e., groups), `group` and `condition`, and two within-subject factors (i.e., repeated- measures factors), `topic` and `problem` (I uploaded the data to pastebin, so everybody should be able to obtain it):               > d <- read.table(url("http://pastebin.com/raw.php?i=4hRFyaRj"), colClasses = c(rep("factor", 6), "numeric"))     > str(d)     'data.frame':   2928 obs. of  6 variables:       $ code     : Factor w/ 183 levels "A03U","A08C",..: 1 1 1 1 1 1 1 1 1 1 ...
       $ group    : Factor w/ 2 levels "control","experimental": 2 2 2 2 2 2 2 2 2 2 ...       $ condition: Factor w/ 3 levels "alternatives",..: 3 3 3 3 3 3 3 3 3 3 ...
       $ topic    : Factor w/ 4 levels "1","2","3","4": 1 1 1 1 2 2 2 2 3 3 ...       $ problem  : Factor w/ 4 levels "AC","DA","MP",..: 3 4 1 2 3 4 1 2 3 4 ...
       $ mean     : num  94.5 94.5 86.5 84.5 80 46.5 73.5 43.5 51 39 ...      The data is from a behavioral experiment in which participants in six groups (2 levels of `group` times 3 levels of `condition`) worked on 16 tasks (for each of 4 `topics` 4 different `problems`). Allocation of participants to group/condition was fully random. Presentation of tasks was random insofar that problem was blocked within topic (i.e., for each topic all problems where presented sequentially), but order of problem and topic was random.   **Update:** The factor identifying the participant (in which topic and problem are nested) is `code`. ## The Problem How can I fit this dataset using `lme`?   (Sidenote: I would also consider using `lme4`, but I am kind of afraid of not having p-values, if there is something easily digestible as p-values, I would also consider `lme4` an option). So far I managed to fit an `lme` model with only one within-subject factor, but not two (see below). ## What I tried I can fit an `lme` model if I have just one within-subject factor:               require(nlme)      m1 <- lme(mean ~ condition*group*problem, random = ~1|code/problem,                 data = d, subset = topic == "1")          anova(m1)                             numDF denDF F-value p-value     (Intercept)                 1   531   12101  <.0001     condition                   2   177      31  <.0001     group                       1   177       2  0.2178     problem                     3   531      35  <.0001     condition:group             2   177       1  0.3672     condition:problem           6   531      24  <.0001     group:problem               3   531       1  0.2180     condition:group:problem     6   531       2  0.0281      This (especially the df) nicely correspond with the results from an standard ANOVA (using `ez`):               require(ez)     ezANOVA(subset(d, topic == "1"), dv = .(mean), wid = .(code), between = .(condition, group), within = .(problem))$ANOVA          Warning: Data is unbalanced (unequal N per group). Make sure you specified a well-considered value for the type argument to ezANOVA().                        Effect DFn DFd     F                             p p<.05     ges     2               condition   2 177 30.69 0.000000000003611248905859672     * 0.13079     3                   group   1 177  1.53 0.217821969825403999321267179       0.00374     5                 problem   3 531 34.85 0.000000000000000000014254103     * 0.10028     4         condition:group   2 177  1.01 0.367225806638525886782531416       0.00492     6       condition:problem   6 531 24.40 0.000000000000000000000000142     * 0.13503     7           group:problem   3 531  1.48 0.217959293081550348203379031       0.00472     8 condition:group:problem   6 531  2.38 0.028119961573665430004664856     * 0.01499      Trying to fit this data with two within-subject factors in `lme` fails (either per code, or per dfs):               m2 <- lme(mean ~ condition*group*problem*topic, random = ~1|code/(problem*topic), data = d)     # fails: Error in getGroups.data.frame(dataMix, groups) :      #  Invalid formula for groups          m3 <- lme(mean ~ condition*group*problem*topic, random = ~1|code/problem/topic, data = d)     # the next model takes some time (probably already an indicator, that it is the wrong model)     # and produces wrong denominator df!          # with both factors as ANOVA     m4 <- ezANOVA(d, dv = .(mean), wid = .(code), between = .(condition, group), within = .(problem, topic))          #effects are the same     all(row.names(anova(m3))[-1] == m4$ANOVA$Effect)          #denominator dfs are not:     anova(m3)$denDF[-1] == m4$ANOVA$DFd          # only for effects with topic:     row.names(anova(m3))[-1][!(anova(m3)$denDF[-1] == m4$ANOVA$DFd)]      **UPDATE** : As the precise error or nesting is somewhat unclear I here provide the equivalent `aov` call (this is the "standard" model via `aov`), which matches the results from `ezANOVA`. The critical error term is `Error(code/(problem*topic))`:               m5 <- aov(mean ~ (condition*group*problem*topic) + Error(code/(problem*topic)), d)     summary(m5)


title: Choosing covariance function

text: Is there a systematic approach or a specific methodology by which we can choose covariance functions that best suit a particular data set? Was carl in his GPML evaluating the suitability of the covariance function when he computed the value of nlml?


title: A MATLAB package to calculate Entropy of a discrete non-stationary (hydrologic) time series

text: I am looking to survey the community to find out any prefered packages researchers, practitioners, or interested modelers are using to calculate entropy (cross-entropy, conditional entropy, etc...) for a given discrete non- stationary time series, let's say for example maximum air temperature measured at an airport at the daily scale for 10 years. The most important aspect of this entropy calculation method I am searching for should consider as few assumptions as possbile in calculating the entropy of a particular time series. I am interested in using entropy metrics to score the utility of a variable within an input variable selection procedure. Unfortunately, many of the packages I have came across only address entropy calculations for discrete categorical inputs, for which my application is not suited. I hope my question is clear, but please do let me know if further information about my specific research problem or any important details is required..


title: What are some good (ideally free) tools to give laymen access to basic statistical techniques?

text: I have been a contract worker for two different companies now, which means I show up at a place for 3-6 months and blow their minds with my awesomeness, but then, to everyone's disappointment, I leave. This is great for me, since I get a new workplace with new problems, but it can be troublesome for the place I left, since they no longer have a person that does something more than look at pivot tables. In an effort to add more value to the places I work, I have started looking into ways to help business people step beyond just pivot tables. Obviously my goal is to not send my nice coworkers off into the woods without any knowledge, but rather to teach them the basics of regression, some graphics options (really pulls a lot of weight in the business world) and maybe some kind of classification (thoughts on what you think is most relevant would be great) in a weekly/monthly meeting. I am picturing an intro level statistics course for an undergrad, which while long ago for many of us still offers quite a bit of power over something like little green arrows if there was growth at all and little red arrows if there was decline next to an excel spreadsheet. I am an R person myself, so I have been looking at the different options there, but I am not going to limit myself to just R, so don't let my list limit you. What I have come across so far: http://en.wikipedia.org/wiki/RExcel I am waiting on IT to give me permission to install this, but it looked promising to me. Excel is the life blood of the places I have worked, so if there is an easy way to get them access to R as an add-in it would help get over the hurdle of learning something new. http://rattle.togaware.com/ I have played with Rattle some and generally had a good time with it. I think it is a bit beyond what I had in mind in terms of statistical difficulty, but it is put together in a very appealing way.


title: Determining which pair of empirical distributions is the most similar

text: I'm looking at three discrete (seven-option) survey questions that were asked on two different surveys (100 < N < 150 for each question). I've paired up the six answer distributions so that each distribution is paired with its counterpart on the other survey (i.e., the one that was generated with the same question). I want to measure which pair is the most similar to each other. Right now I am using a bootstrap Kolmogorov-Smirnov Test (the bootstrap is to deal with the discreteness of the variable) and can generate test statistics and p-values. However, I'm not sure what to do next to figure out which question generated the most similar distributions across the two surveys. This is my first time fiddling with a KS test, so I'm quite in the dark. I understand that confidence intervals might make it hard or impossible to tell which of the two are the most similar, but I wanted to check in case there was a "clear winner." Is this the right procedure to be using? If so, how should I proceed from here? If not, what would you suggest in its stead? (P.S.: I'm using R, so if there's a command or package you would recommend I'd love to hear it.)


title: Model with truncated count response

text: I am trying to select the correct model for analyzing data from a behavioral sciences experiment. The experiment consists of eight trials for each participant, and each experiment generates a score on the interval [0--8] for each participant. Each trial within the experiment can be conceived of as a binomial choice. The ideal model will give a prediction equation for the estimated score, using values for a small number of binary categorical explanatory variables. I am inclined to assume a Poisson distribution for the response variable, with a log link, or a negative binomial distribution if overdispersion is an issue. However, the true distribution is truncated, since any score above 8 is not possible. Is such an approach on the right track? Would it be wiser to treat each trial as a separate outcome, then do a binomial regression with a random effect for participant? (I have seen a similar question on this site, but would appreciate any advice beyond what is given there)


How do I compute a confidence level for a forced ranking survey? For example, say I have 10 categories, each with four associated statements that must be ranked either 1, 2, 3 or 4 with no duplicate ranking within the category (forced ranking). Say I get fifty responses. How do I compute the confidence level/margin of error?


title: Bayesian estimation using Gibbs sampling for financial models

text: I am trying to do Gibbs sampling, from this paper,   www.jds-online.com/file_download/353/JDS-746.pdf This is a CIR financial model, I want to do Gibbs on its parameters:   $$y(t+{\Delta}^{+})=y(t)+(\alpha-\beta y(t)){\Delta}^{+}+\sigma \sqrt{y(t)}{\epsilon}_{t}$$ The conditionals are: Step 1: ${y}_{t,j+1}^{*}|{y}_{t,j}^{*},\theta\ \tilde{\ }\ N({y}_{t,j}^{*}+(\alpha-\beta {y}_{t,j}^{*})\Delta, {\sigma}^{2}\Delta {y}_{t,j}^{*})$ Step 2: $\psi|Y,{Y}^{*},{\sigma}^{2}\ \tilde{}\ N(\mu, {\Lambda}^{-1})$ Step 3: $\sigma|Y,{Y}^{*}, \psi\ \tilde{}\ inverse-gamma(E, F)$   $E=(T-1)(M+1)/2$, and $F=\sum\limits_{t=1}^{T-1}{\sum\limits_{j=0}^{M}{\frac{{{\left( y_{t,j+1}^{*}-\left[ y_{t,j}^{*}+(\alpha -\beta y_{t,j}^{*})\Delta \right] \right)}^{2}}}{2y_{t,j}^{*}}}}$   Repeat Gibbs between step 2 and 3. Notations: $\theta=(\psi,{\sigma}^{2})$ and $\psi=(\alpha, \beta)$ I get negative numbers when sampling from inverse-gamma, I use rinvgamma(). What am I doing wrong?


title: Extrapolating a lowess model

text: Those lines are based on a `lowess` fit of points I have. I would like to be able to extrapolate those lines up to `x=500` for instance. I would like to extend this plot and see if those lines reach a plateau or not. Is there a way to do this ? The CSV data for the blue line is here: http://db.tt/GHhAGLtM ![enter image description here](http://i.stack.imgur.com/iDSx1.png) The idea is fill this plot with extrapolated dashed lines. ![enter image description here](http://i.stack.imgur.com/uHIS6.png)


title: What qualifies as a time series?

text: I have a simple question: does any data gathered over time could be treated as a time series? For instance, if I measure CPU load at certain intervals, would that be a time series? Actually, I'm not interested on actual time of the measurement. I could simply mark the points along x-axis as 1, 2, 3, ... -- what is relevant to me is certain values at certain points of time. Further, could I do time series analysis on the data set if I have irregular observation intervals?


title: How to use scikit-learn's cross validation functions on multi-label classifiers

text: I'm testing different classifiers on a data set where there are 5 classes and each instance can belong to one or more of these classes, so I'm using scikit- learn's multi-label classifiers, specifically `sklearn.multiclass.OneVsRestClassifier`. Now I want to perform cross- validation using the `sklearn.cross_validation.StratifiedKFold`. This produces the following error:               Traceback (most recent call last):       File "mlfromcsv.py", line 93, in <module>         main()       File "mlfromcsv.py", line 77, in main         test_classifier_multilabel(svm.LinearSVC(), X, Y, 'Linear Support Vector Machine')       File "mlfromcsv.py", line 44, in test_classifier_multilabel         scores = cross_validation.cross_val_score(clf_ml, X, Y_list, cv=cv, score_func=metrics.precision_recall_fscore_support, n_jobs=jobs)       File "/usr/lib/pymodules/python2.7/sklearn/cross_validation.py", line 1046, in cross_val_score         X, y = check_arrays(X, y, sparse_format='csr')       File "/usr/lib/pymodules/python2.7/sklearn/utils/validation.py", line 144, in check_arrays         size, n_samples))     ValueError: Found array with dim 5. Expected 98816      Note that training the multi-label classifier does not crash, but the cross- validation does. How must I perform cross-validation for this multi-label classifier? I have also written a second version that breaks down the problem into training and cross-validating 5 separate classifiers. This works just fine. Here is my code. The function `test_classifier_multilabel` is the one giving problems. `test_classifier` is my other attempt (breaking up the problem into 5 classifiers and 5 cross-validations).               import numpy as np     from sklearn import *     from sklearn.multiclass import OneVsRestClassifier     from sklearn.neighbors import KNeighborsClassifier     import time          def test_classifier(clf, X, Y, description, jobs=1):         print '=== Testing classifier {0} ==='.format(description)         for class_idx in xrange(Y.shape[1]):             print ' > Cross-validating for class {:d}'.format(class_idx)             n_samples = X.shape[0]             cv = cross_validation.StratifiedKFold(Y[:,class_idx], 3)             t_start = time.clock()             scores = cross_validation.cross_val_score(clf, X, Y[:,class_idx], cv=cv, score_func=metrics.precision_recall_fscore_support, n_jobs=jobs)             t_end = time.clock();             print 'Cross validation time: {:0.3f}s.'.format(t_end-t_start)             str_tbl_fmt = '{:>15s}{:>15s}{:>15s}{:>15s}{:>15s}'             str_tbl_entry_fmt = '{:0.2f} +/- {:0.2f}'             print str_tbl_fmt.format('', 'Precision', 'Recall', 'F1 score', 'Support')             for (score_class, lbl) in [(0, 'Negative'), (1, 'Positive')]:                 mean_precision = scores[:,0,score_class].mean()                 std_precision = scores[:,0,score_class].std()                 mean_recall = scores[:,1,score_class].mean()                 std_recall = scores[:,1,score_class].std()                 mean_f1_score = scores[:,2,score_class].mean()                 std_f1_score = scores[:,2,score_class].std()                 support = scores[:,3,score_class].mean()                 print str_tbl_fmt.format(                     lbl,                     str_tbl_entry_fmt.format(mean_precision, std_precision),                     str_tbl_entry_fmt.format(mean_recall, std_recall),                     str_tbl_entry_fmt.format(mean_f1_score, std_f1_score),                     '{:0.2f}'.format(support))          def test_classifier_multilabel(clf, X, Y, description, jobs=1):         print '=== Testing multi-label classifier {0} ==='.format(description)         n_samples = X.shape[0]         Y_list = [value for value in Y.T]         print 'Y_list[0].shape:', Y_list[0].shape, 'len(Y_list):', len(Y_list)         cv = cross_validation.StratifiedKFold(Y_list, 3)         clf_ml = OneVsRestClassifier(clf)         accuracy = (clf_ml.fit(X, Y).predict(X) != Y).sum()         print 'Accuracy: {:0.2f}'.format(accuracy)         scores = cross_validation.cross_val_score(clf_ml, X, Y_list, cv=cv, score_func=metrics.precision_recall_fscore_support, n_jobs=jobs)         str_tbl_fmt = '{:>15s}{:>15s}{:>15s}{:>15s}{:>15s}'         str_tbl_entry_fmt = '{:0.2f} +/- {:0.2f}'         print str_tbl_fmt.format('', 'Precision', 'Recall', 'F1 score', 'Support')         for (score_class, lbl) in [(0, 'Negative'), (1, 'Positive')]:             mean_precision = scores[:,0,score_class].mean()             std_precision = scores[:,0,score_class].std()             mean_recall = scores[:,1,score_class].mean()             std_recall = scores[:,1,score_class].std()             mean_f1_score = scores[:,2,score_class].mean()             std_f1_score = scores[:,2,score_class].std()             support = scores[:,3,score_class].mean()             print str_tbl_fmt.format(                 lbl,                 str_tbl_entry_fmt.format(mean_precision, std_precision),                 str_tbl_entry_fmt.format(mean_recall, std_recall),                 str_tbl_entry_fmt.format(mean_f1_score, std_f1_score),                 '{:0.2f}'.format(support))          def main():         nfeatures = 13         nclasses = 5         ncolumns = nfeatures + nclasses              data = np.loadtxt('./feature_db.csv', delimiter=',', usecols=range(ncolumns))              print data, data.shape         X = np.hstack((data[:,0:3], data[:,(nfeatures-1):nfeatures]))         print 'X.shape:', X.shape         Y = data[:,nfeatures:ncolumns]         print 'Y.shape:', Y.shape              test_classifier(svm.LinearSVC(), X, Y, 'Linear Support Vector Machine', jobs=-1)         test_classifier_multilabel(svm.LinearSVC(), X, Y, 'Linear Support Vector Machine')          if  __name__ =='__main__':         main()      I am using Ubuntu 13.04 and scikit-learn 0.12. My data is in the form of two arrays (X and Y) that have shapes (98816, 4) and (98816, 5), i.e. 4 features per instance and 5 class labels. The labels are either 1 or 0 to indicated membership within that class. Am I using the correct format as I don't see much documentation about that?


I have just learned principal component analysis, which seems to be an effective way to do dimension reduction. Now, I would like to know: **Can PCA be extended to account for nonlinear dependencies?** It seems this is possible, based on a few searches to Wikipedia, Google, ScienceDirect, etc. However, the mathematics are quite difficult. PCA starts by finding covariance matrix $C$ of some multivariate data $X$. If I somehow replaced the entries of $C$ so that nonlinear dependencies were taken into account, would the technique then effectively become "nonlinear PCA"? Now if I understand this correctly, it would be so, but in practise finding the entries is way too difficult --> one reverts to approximation/optimization techniques.


title: Lack of memory (memorylessness)

text: The purpose of this question is to gather material about "lack of memory" and to add new ideas about that. If the conditional distribution of a given distribution is equal to unconditional distribution then that distribution possesses a "lack of memory" property. The exponential distribution contains "lack of memory". So my questions are 1) Is there any other distribution which has "lack of memory"? 2) What special properties are associated with exponential or other distributions (if exist) because of "lack of memory"? Please provide interesting examples on applied side which will be helpful for learners to understand lack of memory. One example is here.


title: How to measure the smoothness of a line?

text: I have a function that produces a rising line chart (bank account balance). The ideal case is when the chart is a straight rising line, but in the most cases the line is rising and falling in different degree. What statistical parameter to use to measure how smooth the line is? Currently I take several sampling points and measure the distance from the function to the theoretical straight line between the start and the end point. This method works more or less, but it misses eventual picks or drops between the sampling points. May I use Standard Deviation in some way? Current code:               private bool IsBalanceLineSmooth(int checkPoints, double maxPercentDeviation )     {         for (int i = 1; i <= checkPoints; i++)         {             double netBalance        = Backtester.NetMoneyBalance;             double startBalance      = Backtester.MoneyBalance(0);             double checkPointBalance = Backtester.MoneyBalance(i*Bars/(checkPoints + 1));             double targetBalance     = startBalance + i*(netBalance-startBalance)/(checkPoints+1);             double minBalance        = targetBalance*(1 - maxPercentDeviation);             double maxBalance        = targetBalance*(1 + maxPercentDeviation);             if (checkPointBalance < minBalance || checkPointBalance > maxBalance)                 return false;         }         return true;     }      ![Example Chart](http://i.stack.imgur.com/pS7zC.png) Generally I need to evaluate the line with an integer number. Lets say if 100 is a straight line and 0 is hmm.. fully jagged line. I need a number (for example 80) that shows the minimum acceptable smoothness value. I need users of the software to be able to set this number as a an acceptance criteria. I'm not sure if the smoothness criteria have to represent the maximum deviation (like in my code) or an average value for the whole line. Probably the second is the better solution. A C# code example is appreciated.


Can the general linear model be solved by generalized least squares? The model that generalized least squares (GLS) can solve doesn't assume anything on the covariance matrix of the errors, which I think is the same as a general linear model?


title: Which statistical test to use with multiple response variables and continuous predictors?

text: I need help to define which test to use. I want to know which biochemical parameters are influencing the reproductive gonad stage in a scallop. Gonads have a great variability therefore each gonad stage was determined by taking 50 random points in each gonad and scoring the percent in three categories: "spawning", "mature" and "atresia". My explanatory variables are a set of environmental variables: temperature and chlorophyll and biochemical composition inside the gonad and the muscle of the scallop: glycogen, protein, lipids. Can anyone please tell me how can I analyze these data?


title: Time variation in coefficients

text: Given $x_t, y_t$ ($t=1,\ldots,240$), I want to estimate $y_t = \alpha_t + \beta_t x_t$ and test $H_0: \alpha_1=\ldots=\alpha_T=0$. It is crucial to allow for time variation in the regression coefficients. For instance, $x_t$ is the monthly market return in excess of the riskfree rate, and $y_t$ are excess returns of a stock portfolio. I consider to use one of the following two methods.   * **Kalman filter** This requires a parametric model for the processes of $\alpha_t$ and $\beta_t$. Is normality for the coëfficient processes ok?   * **Feasible GLS** Estimate for $s = 60,\ldots,240$ the overlapping rolling estimation windows $$ \left[\begin{array}\\ y_{s-59}\\ \vdots\\ y_{s}\end{array}\right] = \left[\begin{array}\\ 1\\ \vdots\\ 1\\ \end{array} \right] \alpha_s + \left[\begin{array}\\ x_{s-59}'\\ \vdots\\ x_s'\\ \end{array} \right] \beta_s + \left[\begin{array}\\ e_{s,s-59}\\ \vdots\\ e_{s,s}\\ \end{array} \right] $$ Since the error terms are dependent for the same $y_t$ in different windows ending at, say $s1$ and $s2$, define $\sigma_{s_1s_2}:=Cov(\varepsilon_{s1,t},\varepsilon_{s2,t})$ for all common $t$ in the two regression windows. Further, the variance $\sigma_{t,t}:=Cov(\varepsilon_{s,t},\varepsilon_{s,t})$ is allowed to be nonzero and identical for all regression windows $s$. All other covariance terms are supposed to be zero. First, the covariance matrix is estimated with residuals from OLS estimates. Second, the GLS estimator is used with this covariance matrix. Remarkably, I could not find papers in the finance literature which compare the two methods. Is either of the two methods preferred? Or do better methods exist? Thanks in advance.


title: Are GAMMs/GLM the best choice for calculating number of germs on hands?

text: We would like your opinion on whether GAMMs are a good option and how best to go about implementing for the following: During a period of patient care, a nurse will accrue $Y$ colonies of bacteria on their hands. During every hand contact ($n$) (with surface area $A$) with a surface this number may increase a percentage ($\lambda$ %) of what she touches ($V$). However potentially this may also decrease as some percentage ($\beta$ %) of the bacteria are removed. After they finish the care episode, they may wash their hands with efficacy $h$. The model is hence dependent on $Y=Y(n,A,V,\lambda,\beta,h)$ **We think:** > $n$=empirical, different for every nurse, $h=lognormal~(1.5,0.1)$, > $\lambda=\Gamma(15,3)$, $\beta=$empirical non-negative. > $A=lognormal~(7,1.9)$ and $V$=empirical, $h=\Gamma(5.91,0.4)$. We'd like to know if GAM fitting is an appropriate way of estimating $Y$. Please let us know if you require further clarification. Regards,


title: Significant difference from regression confidence intervals

text: I have a question about statistical significance in relation to confidence intervals from linear regression. I'm obviously far from a stats expert, and I've been searching for the answer to this, probably simple, question for a while now without any luck. I've made an example to clarify my question: I'm interested in looking at the treatment effect of doing a change (e.g. spraying with pesticide) on one area, and use an untreated area for control. Before the treatment a "calibration line" is established between the two areas as the correlation between some observed response (e.g. annual crop yield, white open circles below, with regression line and 95% confidence intervals drawn) ![Link to example plot: \(sorry, not rep to post image\)http://imgur.com/M95dqEk](http://i.stack.imgur.com/43d0P.jpg) [Fig: X and Y axis show the same response (e.g. annual crop yield), but for the control area(X axis) and for the treatment area (Y axis). Each data point is then the annual crop yield for both the control and the treatment area, a total of 10 years of data] After the treatment this response is measured again (red open triangles), and the "treatment effect" is defined as the difference between the observed response in the treatment and the predicted response (the regression line). My question is if you can say that the treatment effect is statistically significant if the data point is outside of the 95% confidence interval of the calibration regression? And why/why not ? (so in the plot example 3 of the observed treatment effects are significantly different from the predicted response to a 95% confidence level (p=0.05)?) Thanks Edit1, additional question: Would prediction intervals, instead of confidence intervals, be more suitable to describe whether there has been a change in the relationship between the two areas after treatment? Edit2: Would it be right so say that the confidence intervals can be used to check if the mean of the treatment effects are significantly different (and, as @Glen_b suggests, use the regression line/confidence line for the treatment points instead of single points. But when talking about whether a single sample is significantly different (as in my comment below to Glen_b) it is better to use the prediction interval?


title: Assumptions of generalised linear model

text: I have made a generalised linear model with a single response variable (continuous/normally distributed) and 4 explanatory variables (3 of which are factors and the fourth is an integer). I have used a Gaussian error distribution with an identity link function. I am currently checking that the model satisfies the assumptions of the generalised linear model, which are:   1. independence of Y   2. correct link function   3. correct scale of measurement of explanatory variables   4. no influential observations My question is: **how can I check that the model satisfies these assumptions?** The best starting point would seem to be plotting the response variable against each explanatory variable. However, 3 of the explanatory variables are categorical (with 1-4 levels), so what should I be looking for in the plots? Also, do I need to check for multicollinearity and interactions amongst explanatory variables? If yes, how do I do this with categorical explanatory variables?


title: Why is clutter problem intractable for large sample sizes?

text: Suppose we have a set of points $\mathbf{y} = \{y_1, y_2, \ldots, y_N \}$. Each point $y_i$ is generated using distribution $$ p(y_i| x) = \frac12 \mathcal{N}(x, 1) + \frac12 \mathcal{N}(0, 10). $$ To obtain posterior for $x$ we write $$ p(x| \mathbf{y}) \propto p(\mathbf{y}| x) p(x) = p(x) \prod_{i = 1}^N p(y_i | x). $$ According to Minka's paper on Expectation Propagation we need $2^N$ calculations to obtain posterior $p(x| \mathbf{y})$ and, so, problem becomes intractable for large sample sizes $N$. However, I can't figure out why do we need such amount of calculations in this case, because for single $y_i$ likelihood has the form $$ p(y_i| x) = \frac{1}{2 \sqrt{2 \pi}} \left( \exp \left\{-\frac12 (y_i - x)^2\right\} + \frac{1}{\sqrt{10}} \exp \left\{-\frac1{20} y_i^2\right\} \right). $$ Using this formula we obtain posterior by simple multiplication of $p(y_i| x)$, so we need only $N$ operations, and, so we can solve this problem for large sample sizes exactly. I make numerical experiment to compare do I really obtain the same posterior in case I calculate each term separately and in case I use product of densities for each $y_i$. Posteriors are same. See ![enter image description here](http://i.stack.imgur.com/3ZgRn.png) Where am I wrong? Can anyone make it clear to me why do we need $2^N$ operations to calculate posterior for given $x$ and sample $\mathbf{y}$?


title: Kolmogorov-Smirnov and bootstrap

text: The KS test uses the statistic $$ D_n=\sup_x |\hat{F}_n(x)-F_0(x)| $$ where $F_0(x)$ is the distribution to be tested and $\hat{F}_n(x)$ the empirical distribution. Under the null hypothesis $D_n$ is distributed according to the Kolmogorov distribution and the KS-test is performed testing $D_n$ against the Kolmogorov distribution. My question concerns the bootstrapped version , i.e., $$ D^*_n=\sup_x |\hat{F}_n^*(x)-F_0(x)| $$ where $\hat{F}_n^*(x)$ is obtained through bootstrapping from $\hat{F}_n(x)$. Given the data, which is the distribution of $D^*_n$ generated by bootstrap? Which is the relationship between this distribution and the kolmogorov distribution for $n\rightarrow \infty$? Can I derive a KS-test based on the distribution $D^*_n$ given the data? Is there some known result about that? Thank you for any help.


title: How to combine step(x) and dbern(p) in JAGS/WinBUGS?

text: The problem involves trying to characterize the probability:               P.f = Pd*Pr{t1 < t2}      using jags or WinBUGS. The issue is the last term where both t1 and t2 are random variables. A sample, stripped down, model using this is given below and hopefully provides insight into what I am trying to do. As expected, I get an error indicating that I am redefining the variable 'y'. Searching for a trick to get past this is proving difficult. Any insight would be appreciated! _FWIW, I've also posted a similar question on the WinBUGS list, but cross- posted here with the hope of reaching a wider audience._ * * *               model {                  for(j in 1:N) {             t.1[j] ~ dweib(alpha.1,lambda.1)              t.2[j] ~ dweib(alpha.2, lambda.2);             p.det[j] ~ dbeta(a,b);             y[j] <- step(t.2[j]-t.1[j]);                         y[j] ~ dbern(py);         }            alpha.1 ~ dgamma(0.3,0.0001);          lambda.1 ~ dnorm(0., 10000.);         alpha.2 ~ dgamma(0.3,0.0001);          lambda.2 ~ dnorm(0., 10000.);            py ~ dbeta(0.3,0.3);         a ~ dgamma(1, 0.01)          b ~ dgamma(1, 0.01)          pd ~ dbeta(a,b)          p.f<- pd*py     }


title: Classify unlabeled data

text: I want to classify data into two classes based on three parameters given. My data is web server log and I want to classify it into potential user class and non potential user. The parameters are the time spent by the user in the website, the number of pages visited, and the number of transaction per session. All the parameters are discretized into categorical value. what is the best and easiest method to do that?


title: Convergence in probability, $X_i$ IID with finite second moment

text: Let $\{X_i\}_{i\geq 1}$ be IID with finite second moment, and $$ Y_n = \frac{2}{n(n+1)}\sum_{i=1}^n \,i\cdot X_i \, , \qquad n\geq 1 \, . $$ Could you please tell me how can I show that $Y_n$ converges in probability to $\mathrm{E}[X_1]$? I'm thinking Kolmogorov Convergence criterion. But seems like I cannot prove it using that. Any suggestion?


title: When is the IID Assumption too strict for the bootstrap?

text: The Bootstrap (Efron 1979) assumes that the data are IID. Obviously, if we have time series data then we probably cannot make that assumption unless we have a special case that we a time series of IID noise. Under what other, non-obvious conditions is this assumption invalid? How about the following scenario? I have a table of events, which were triggered by different individuals at different points in time. At certain times, it is more likely that an event will be triggered. Some events are more common that others. However if individual A triggers an event, this has no influence on whether or not individual B will trigger an event. Some types of events are more popular among most individuals. It is attributes of the events that drive their popularity. Some events only occur at certain times. Events are assigned "minute buckets", so that their actual time is not used. These minute buckets are of varying sizes. I am trying to answer the following question: "what attributes of these events attract more users to trigger them?" An individual may trigger multiple events. There could be two cases here:   1. Multiple events could be so compelling that they will both be triggered   2. An individual could have persistence to their actions This is probably due to two kinds of individuals. The data described above would not be called "panel data", right?


title: Which type of residuals to use for the Durbin-Watson test (autocorrelation)

text: I want to check if there is residual autocorrelation in my model and the test for this is the Durbin-Watson test. I am using R and my question is if it makes a difference which type of residuals one uses when applying the test, since there are several types of residuals like deviance, pearson, working...


title: Dealing with least and maximum score in IRT

text: Generally when I do IRT modelling (using R software, 2PL model, ltm package) I remove records of students who have scored either 0 or 100%. The logic being that we do not have enough information about the ability of these students. But this way I cannot allocate any score to students who scored 0 or 100% I have tried running the model with such records also. The model converges. But I am not sure if it is correct to do it. I am getting some values for zero and 100% scorers but I am not sure how it is being calculated in the ltm package.


title: How to use auto.arima to impute missing values

text: I have a zoo series with many missing values. I read that `auto.arima` can impute these missing values? Can anyone can teach me how to do it? thanks a lot! This is what I have tried, but without success:               fit <- auto.arima(tsx)     plot(forecast(fit))


title: Determining Whether Values Within One Standard Deviation Meet Some Condition

text: Suppose I have a set of objects, each of which can take on values within some unknown range, and each of which only has one sample of 50+ values associated with it. The values of each observation is not known. In order to determine how likely it is for the value of an object in a set of objects $O$ to have a negative value I opted to use the difference $\sigma(X_i)-|E(X_i)|$ with the assumption that $\forall$ $X_i$ $\in $ $ O$, $X_i$ is i.i.d. Gaussian. In this way I can find out whether one standard deviation from a known mean value for each $X_i$ is negative. If it is, I can use normality to determine how likely it might be for this object to be negative. Is this a good way to go about determining the likelihood, or is it possibly too naïve? I worry that I may be making some errors.


title: Line of best fit (Linear regression) over vertical line

text: I want to get a line of the best fit which is a line that passes as close as possible to a set of points defined by coordinates point_i = (X_i, Y_i). When I apply linear regression, I have a special case where the line do not correspond to what I want. It is the case where the data-points are aligned over a vertical line (i.e. Y values do not depend on X values) like the following example for instance (or see the points on image bellow), e.g:               Y = [1,2,3,4,5,8,10,11,15,13,21,22,23,24,25,28,30,31,35,33]     X = [1,1,1,1,1,1,1 ,1 ,1 ,1 ,1 ,1 ,1 ,1 ,1 ,1 ,1 ,1 ,1 ,1 ];      In this case I expect to get a line of best fit which is vertical and passes as close to a set of points, but when I get from regression is a horizontal line ! However, what I'm interested in, is not the error over the Y axis, but I'm interested in the error to minimize which indicate how far (according to distance) are the points from the line of best fit that we get (e.g. the sum of shortest distance from each point to the line). Is there any why to get such a line and to get the corresponding score which indicates how far are the points this line ? I think "linear regression" is not the same thing as "line of best fit", and what I want to get is "line of best fit". ![enter image description here](http://i.stack.imgur.com/aTW7i.png) EDIT: I want to estimate line of best fit for another purpose, not to evaluate how much X and Y are correlated (it is more for a pattern recognition problem, where my data-points are some pixels defined by their coordinates x,y on the image). So I want to fit the line even if X values are not correlated to Y values. And if I just always swap X and Y to be able to fit the line, the next time if I have points lying on a horizontal line (Y = [3,3,3,3], X=[3, 6, 23, 30] for example) the if they are swapped I'll get again the same problem. I want to always fit the line whatever is the correlation between W and Y (even if there is no correlation, like in the example I provided). So to be more clear: Whatever is the correlation between X and Y, my objective is not to predict some y values based on new values of x, my objective is rather just to get a line which passes as close as possible to all points that we have (even if Y are independent of X ...).


title: Test for significance of correlation matrix

text: If one wishes to test if the correlations in a correlation matrix are statistically significant as a whole group, one can perform a likelihood ratio test of the hypothesis that the correlation matrix is equal to the identity matrix. The ratio of the restricted and unrestricted likelihood functions is $\alpha = |R|^{N/2}$ , where |R| is the determinant of the correlation matrix (Morrison, 1967). The test statistic is therefore $-2\log(\alpha)$, which is distributed as $\chi^2$ with $\frac{1}{2}p(p-1)$ df. **My question is, how do you perform the calculation of the observed value $-2\log( \alpha)$?**


title: how to estimate CTR (ctr-click-through-rate)?

text: How many times should a banner be shown to estimate its click-through-rate (CTR)? For example, if a banner was shown $x$ times, and was clicked $y$ times. $$\text{CTR} = \frac{y}{x}$$ How could I evaluate inaccuracy of this value?


title: Quantify strength of association of two continuous variables while controlling for random effects

text: I have a data set from a repeated measures experimental design with different sets of stimuli. I want to know how strong the association between the continuous dependent variable and the continuous predictor is while accounting for the interindividual and interstimulus variation. My `lmer` model description in `R` looks like this               dv ~ pred + (1 | subject) + (1 | stimulus)       _Question 1_ : I understand that it is non-trivial to calculate R squared for random intercept-slope models. Is the same true for random intercept models? Is there an R-implementation of any of the available methods? _Question 2_ : If I z-transform my dependent variable and predictor, will the parameter estimate for the fixed effect reflect the strength of the association such as it would in an ordinary regression model? **Update** : I think I phrased this questions too vaguely. I was wondering if scaling variables would yield standardized regression estimates. I found this question has been answered before in another question. _Question 3_ : Is there an entirely different/more appropriate way to quantify the association strength while controlling for the interindividual and interstimulus variation?


title: Longitudinal item response theory models in R

text: I'm trying to fit longitudinal item response theory (IRT) models in R. I have a test that was administered at multiple measurement occasions. I'd like to examine individuals' growth curves of factor scores (i.e., ability levels) from graded response models (GRMs). I have used the ltm package in R to fit cross-sectional GRM models in IRT, but it's not clear to me how (or whether it's even possible in ltm) to extend the models to handle repeated measures of the same items across time. How can I fit growth curves to longitudinal GRM factor scores to see changes in means/variances of ability levels over time? If this isn't possible with the ltm package, what packages/functions permit this? Specific code examples would be especially appreciated. Here are empirical examples similar to what I'm looking to do (except the examples use a Rasch model for binary items where I'm using GRMs for polytomous data): McArdle, J. J., Grimm, K. J., Hamagami, F., Bowles, R. P., & Meredith, W. (2009). Modeling life-span growth curves of cognition using longitudinal data with multiple samples and changing scales of measurement. _Psychological Methods_ , _14_ , 126-149. doi: 10.1037/a0015857 McArdle, J. J., & Grimm, K. J. (2011). An empirical example of change analysis by linking longitudinal item response data from multiple tests. In A. A. Davier (Ed.), _Statistical models for test equating, scaling, and linking_ (pp. 71-88): Springer: New York. For example, I'd like to be able to estimate and plot individuals' growth curves to examine mean-level change in ability levels over time (from McArdle & Grimm, 2011): ![Individuals' Growth Curves In Ability Levels Over Time](http://i.stack.imgur.com/YFXf0.jpg) And, I'd like to be able to estimate an average or prototypical growth curve for the sample (from McArdle & Grimm, 2011): ![Average/Prototypical Growth Curve In Ability Levels Over Time](http://i.stack.imgur.com/wjH3w.jpg) Here's a simulated data set with 20 polytomous items (1-3 response scale) at 3 different time points:               library(mirt)     library(mvtnorm)          set.seed(1)          numberItems <- 20     numberItemLevels <- 2     sampleSize <- 1000          a <- matrix(rlnorm(numberItems, .2, .2))     d <- matrix(rnorm(numberItems*numberItemLevels), numberItems)     d <- t(apply(d, 1, sort, decreasing=TRUE))          Theta <- mvtnorm::rmvnorm(n=sampleSize, 0, matrix(1))          t1 <- simdata(a, d, N=sampleSize, itemtype="graded", Theta=Theta)     t2 <- simdata(a, d, N=sampleSize, itemtype="graded", Theta=Theta+.5)     t3 <- simdata(a, d, N=sampleSize, itemtype="graded", Theta=Theta+1)          dat <- data.frame(t1, t2, t3)


title: ROC vs Accuracy

text: I've designed a 4 classifiers which perform pretty decently (all of them are above 90% in accuracy). However, they don't have similar AUC for their respective ROC curves (obviously, it doesn't have to be). If I were to use these classifiers in real-time data, which one do I choose based on the following result **Classifier A** : Accuracy: 100%, AUC: 84% **Classifier B** : Accuracy: 95%, AUC: 83% **Classifier C** : Accuracy: 100%, AUC: 69% **Classifier D** : Accuracy: 100%, AUC: 77%


title: Derivation of equations in kriging

text: I have some confusion regarding some derivations in the equation of kriging of wiki article http://en.wikipedia.org/wiki/Kriging. It says that kriging error is given by $$ \sigma_k^2(x_0)=Var(\hat{Z}(x_0)-Z(x_0)).......(1) $$ $$ =E((Z(x_0) - \hat{Z}(x_0))^2)....(2) $$ $$ =\sum_{i=1}^n\sum_{j=1}^nw_i(x_0)w_j(x_0)c(x_i,x_j) + Var(Z(x_0)) -2\sum_{i=1}^{n}w_i(x_0)c(x_i,x_0) ......(3) $$ I am confused how the equivalent expression 3 from 2 was derived.


title: Creating censored DV for survival analysis (in long form)

text: I am working on a project that requires me to use survival analysis. Specifically, I need to use discrete-unit survival analysis (so Cox regression or other methods that assume my DV is continuous are out of the question). I have transformed my dataset into long form, and I need to create a censored DV ("dv") by (STEP 1) finding the last time a variable ("q") =1; (STEP 2) recoding my "dv" as "0" at time+1; (STEP 3) recoding my "dv" as "NA" at time+2...n; (STEP 4) deleting all of the observations (rows) for when "dv"=NA. Before I do this, it seems like I should create to be that I should create a DV that is assigned to "1" by default. My data would look like this after step 3:                id  time q  dv         1    1   1   1      1    2   0   1      1    3   1   1      1    4   0   0      1    5   0   NA      1    6   0   NA      2    1   1   1      2    2   1   1      2    3   0   0      2    4   0   NA      2    5   0   NA      2    6   0   NA      Please note that "q" can vary between "0" or "1" over time (See the observation for id=1 at time=2), but due to the nature of survival analysis, "dv" cannot. What I need to do is create a variable that finds the LAST time that "q" changes between "1" and "0", and then is censored accordingly. After Step 4, my data should look like this:                id  time q  dv         1    1   1   1      1    2   0   1      1    3   1   1      1    4   0   0      2    1   1   1      2    2   1   1      2    3   0   0      Thanks for any and all help possible!


title: SVD of correlated matrix should be additive but doesn't appear to be

text: I'm just trying to replicate a claim made in the following paper, Finding Correlated Biclusters from Gene Expression Data, which is: > Proposition 4. If $X_{IJ}=R_{I}C^{T}_{J}$. then we have: > > i. If $R_{I}$ is a perfect bicluster with additive model, then $X_{IJ}$ is a > perfect bicluster with correlation on columns;   >  ii. If $C_J$ is a perfect bicluster with additive model, then $X_{IJ}$ is a > perfect bicluster with correlation on rows;   >  iii. If both $R_I$ and $C_J$ are perfect biclusters with additive model, > then $X_{IJ}$ is a perfect correlated bicluster. > > These propositions can be easily proved... ... but of course, they don't prove it. I'm using some of the simple examples in the paper plus base + custom R code to see if I can demonstrate this Proposition.               corbic <- matrix(c(0,4,-4,2,2,-2,6,0,4,-8,16,-2,-2,10,-14,4), ncol=4)      (from Table 1F) some custom code to convert standard X = $UdV^T$ svd form to $X=RC^{T}$ as described in the paper:               svdToRC <- function(x, ignoreRank = FALSE, r = length(x$d), zerothresh=1e-9) {     #convert standard SVD decomposed matrices UEV' to RC' form     #x -> output of svd(M)     #r -> rank of matrix (defaults to length of singular values vector)                 # but really is the number of non-zero singular values     #ignoreRank -> return the full decomposition (ignore zero singular values)     #zerothresh -> how small is zero?              R <- with(x, t(t(u) * sqrt(d)))         C <- with(x, t(t(v) * sqrt(d)))              if (!ignoreRank) {             ind <- which(x$d >= zerothresh)         } else {             ind <- 1:r         }              return(list(R=as.matrix(R[,ind]), C=as.matrix(C[,ind])))     }      apply this function to the dataset:                > svdToRC(svd(corbic))     $R                [,1]       [,2]     [1,]  0.8727254 -0.9497284     [2,] -2.5789775 -1.1784221     [3,]  4.3244283 -0.7210346     [4,] -0.8531261 -1.0640752          $C               [,1]       [,2]     [1,] -1.092343 -1.0037767     [2,]  1.223860 -0.9812343     [3,]  3.540063 -0.9586919     [4,] -3.408546 -1.0263191      Unless I'm hallucinating, this matrices are not additive, even though corbic exhibits perfect correlation between rows and columns. It seems strange that the example they provide does exhibit the property they said it should... unless I'm missing some kind of pre- or post- svd transformation step?


title: How to compare first to second half of measurements in a counterbalanced behavioral experiment?

text: I want to compare measurements in two halves of a behavioral experiment to examine for learning effects. My first course of action is to split the measurements in the two halves and then repeat the original MANOVA for this comparison. I am asking to find out if there are standard procedures to conduct such a comparison. Or if there are known pitfalls to avoid in such a comparison. I usually work in R.


title: How do you use the EM algorithm to calculate MLEs for a latent variable formulation of a zero inflated Poisson model?

text: The log likelihood of the zero-inflated poisson regression model is $$L(\gamma,\mathbf{\beta}; \mathbf{\gamma}) = \sum_{y_i=0} \log(e^{G_i \gamma}+\exp(-e^{\textbf{B}_i \mathbf{\beta}})) +\sum_{y_i >0} (y_i \textbf{B}_i \mathbf{\beta}-e^{\textbf{B}_i \mathbf{\beta}})-\sum_{i=1}^{n} \log(1+e^{G_{i} \gamma})-\sum_{y_i >0} \log(y_{i}!)$$ Assuming that we could observe $Z_i = 1$ when $Y_i$ is from the perfect, zero state and $Z_i = 0$ when $Y_i$ is from the Poisson state the log-likelihood would be $$L(\gamma,\mathbf{\beta}; \mathbf{\gamma}, \textbf{z}) = \sum_{i=1}^{n} \log(f(z_i|\mathbf{\gamma}))+\sum_{i=1}^{n} \log(f(y_i|z_i, \mathbf{\beta}))$$ $$ = \sum_{i=1}^{n} (z_{i} \textbf{G}_i \gamma-\log(1+e^{G_{i} \gamma}))+ \sum_{i=1}^{n} (1-z_i)(y_{i} \textbf{B}_i \beta-e^{\textbf{B}_i \beta}) -\sum_{i=1}^{n}(1-z_i) \log(y_{i}!)$$ But aren't latent variables unobservable? The purpose is the maximize the first log-likelihood. But we have to introduce latent variables and derive a new log-likelihood. Then using the EM algorithm, we can maximize the second log-likelihood. But this assumes that we know that either $Z_i = 0$ or $Z_i = 1$? **Added.** Suppose $ \textbf{Y} = (Y_1, \dots, Y_n)'$ are independent and $$\eqalign{ Y_i = 0 & \text{with probability} \ p_i+(1-p_i)e^{-\lambda_i}\\ Y_i = k & \text{with probability} \ (1-p_i)e^{-\lambda_i} \lambda_{i}^{k}/k! }$$ Also suppose the parameters $\mathbf{\lambda} = (\lambda_1, \dots, \lambda_n)'$ and $\textbf{p} = (p_1, \dots, p_n)$ satisfy $$\eqalign{ \log(\mathbf{\lambda}) &= \textbf{B} \beta \\ \text{logit}(\textbf{p}) &= \log(\textbf{p}/(1-\textbf{p})) = \textbf{G} \mathbf{\lambda}. }$$


title: Cross entropy and the fitness function

text: I'm training some neural networks using NEAT C++, and would like to use the _cross-entropy_ error function: $$E = -t\ln(y) -(1-t)\ln(1-y)$$ to train a network for a two-class classification problem. NEAT C++ is a C++ library that allows one to use a special kind of neuro evolution (namely that of _augmenting topologies_ ) to train the networks. Yet to be able to use it, one needs to define a fitness function (and not an error function, such as the one above). ### Question: How can I construct a fitness function from the cross-entropy error function? _Everything below this line is wrong I guess because $E\in [0,\infty]$, and thus it's impossible to define a maximum. How should I do it then?_ * * * I was thinking about this: calculate the maximum error $E_{\text{max}}$ somehow and define the fitness, given a certain network output $y$, as: $$f(y)=E_{\text{max}} - E(y)$$ This way the fitness will be $0$ when $E(y)=E_{\text{max}}$ and it should reach a maximum if $E(y)=0$. Is this the right method? If so how do I calculate the maximum error $E_{\text{max}}$? _note: NEAT does not allow the use of negative fitness functions, so using $f(y)=-E(y)$ would not work. Also I would like to refrain from using arbitrary stuff such as $f(y)=10000000 - E(y)$_.


title: How to bootstrap prediction intervals for "customized" regression models in R

text: Are there functions in R that could help me do the following? We have a special type of regression which is called Geometric Mean Regression. We have done some search and found the following: https://stat.ethz.ch/pipermail/r-help/2011-July/285022.html The question is: how to do the statistical inference on GMR results? More specifically, we are looking for the prediction interval: Lets say we regress y1, y2, ..., yn onto x1, x2, ..., xn: we would like to know what's the prediction interval for a new data point: x_new=x1+x2+x3 (i.e. the new data point is the sum of the existing first three data points) In ordinary linear regression, we could derive prediction interval for an in- sample data point as well as a new data point... For our x_new=x1+x2+x3, we can derive formulas for the prediction interval. But for the above customized regression, how do we obtain the prediction intervals? * * * Are there functions in R that can help us do this? We are thinking of using bootstrapping, etc. Are there functions in R help us on this? Thanks a lot! * * * I did some googling and researching... Reading the following article, http://www.ecd.bnl.gov/pubs/BNL-79819-2008-JA.pdf It seems that once we estimate the parameters of the bivariate normal distribution, then we can plug into the formula of conditional distribution of Y|X=x1+x2+x3 ? http://en.wikipedia.org/wiki/Multivariate_normal_distribution My question is: Is it a correct procedure to do the following: Step 1: estimate the parameters of the bivariate normal distribution; Step 2: plug the estimated parameters into the Y|X=x1+x2+x3 formula and get the 95% quantile of it? Do I need to repeat Step 2 many times following the bootstrapping procedure? Or one shot of Step 2 is enough? I got very much confused... Any thoughts? Thanks a lot!


title: Politically correct way of ordering stacked barplots?

text: What is the politically correct way of ordering stacked barplots that contain gender statistics? For example, in a stacked barplot for the number of grantees of a fellowship program, should one put male or female on the bottom? How about left/right if it's on horizontal bars?


title: How to check if a time series is I(1) in R?

text: I'm testing the cointegration of two time series of stock prices using adfTest from fUnitRoots, but first I need to check if the series are I(1). How can I check if a time series is I(1)?


title: PCA and varimax references?

text: After doing PCA it is common that the first component describes the largest part of variability. This is important in i.e. study of body measurements where it is commonly known (Jolliffe, 2002) that PC1 axis captures size variation. My question is whether PCA scores after varimax rotation retain the same properties or are they different as mentioned in this topic? Since I need PCA scores for further statistical analyses I am wondering if varimax is needed and does it in fact disrupt the representation of real sample variability so that individual scores on rotated axes are uninformative or lead to miss-interpretation of reality? Also could someone suggest some other references on this topic? **WORKFLOWS in R**   1. PCA (FactoMineR or prcomp) -> Extract individual scores -> Enter scores in the lm   2. PCA (FactoMiner or prcomp) -> Varimax on loadings matrix -> calculate the individual scores -> enter scores in the lm   3. FA (psych, varimax and pca extraction method) -> extract individual scores -> Enter scores in the lm Now, without rotation (1.) percentages of explained variability are i.e. 29.32, 5.6, 3.2, on the first three axes. 2. and 3. solutions yield similar percentages on the first three factors i.e. 12.2, 12.1, 8.2. Off course 1. solution tends to push all high variable loadings on the first axis, while 2. and 3. tend to distribute loadings between axes (which is the reason for rotation). I wanted to know if these three workflows are essential the same since individual scores are different on rotated vs. unrotated axes?


title: Most efficient way of sampling a product of distributions $k \cdot f(x) \cdot g(x)$, given that PDFs $f(x)$ and $g(x)$ can easily sampled?

text: Is there an efficient approach tailored to sampling from a PDF $k \cdot f(x) \cdot g(x)$ (where $k$ is a normalizing constant) that would perform better than naive Metropolis, slice sampling, etc.? Can we exploit the fact that both $f(x)$ and $g(x)$ can already easily be sampled from? An example of this might be sampling $X \sim k \cdot f(x|a) \cdot f(e^x|b)$ where $f(x|\theta)$ is the normal distribution, and $a$ and $b$ are known. Note that the terms aren't independent, so Gibbs sampling wouldn't necessarily help (this is actually meant to be part of a Gibbs sampler).


title: Duplicate detection/classification using TF/IDF and cosine similarity metrics

text: I am working on text classification using TF/IDF and cosine similarity metrics. In computing the similarity between two records (objects) in a database, TF/IDF assigns weight vector for each record by assigning weight score to each attributes (strings) in a record. Cosine similarity is then used to measure the similarity of the two records from their weight vectors generated by TF/IDF. I realised that TF/IDF does not give accurate weight score when typographical (spelling) errors are present in strings, this is because TF/IDF assigns accurate scores to words present in the dictionary (derived from the training set), hence wrongly spelt words are poorly scored, and results in the computation of inaccurate similarity between records. I also realised that character-based similiarity metrics (affine gap, Jaro- winkler, levenshtien distance etc) are effective in computing the similarity between strings with typographical (spelling) errors, but are not effective when words are out of order in strings. Hence, in order to measure the similarity between documents (records), I would like to enhance the effectiveness of TF/IDF to deal with typographical errors by using the error correcting ability of a character based similarity metric such as Levenshtein distance or affine gap or even Jaro distance. Any suggestions will be highly appreciated.


title: Visualizing results from multiple latent class models

text: I am using latent class analysis to cluster a sample of observations based on a set of binary variables. I am using R and the package poLCA. In LCA, you must specify the number of clusters you want to find. In practice, people usually run several models, each specifying a different number of classes, and then use various criteria to determine which is the "best" explanation of the data. I often find it very useful to look across the various models to try to understand how observations classified in model with class=(i) are distributed by the model with class = (i+1). At the very least you can sometimes find very robust clusters that exist regardless of the number of classes in the model. I would like a way to graph these relationships, to more easily communicate these complex results in papers and to colleagues who aren't statistically oriented. I imagine this is very easy to do in R using some kind of simple network graphics package, but I simply don't know how. Could anyone please point me in the right direction. Below is code to reproduce an example dataset. Each vector xi represents the classification of 100 observations, in a model with i possible classes. I want to graph how observations (rows) move from class to class across the columns.               x1 <- sample(1:1, 100, replace=T)     x2 <- sample(1:2, 100, replace=T)     x3 <- sample(1:3, 100, replace=T)     x4 <- sample(1:4, 100, replace=T)     x5 <- sample(1:5, 100, replace=T)          results <- cbind (x1, x2, x3, x4, x5)      I imagine there is a way to produce a graph where the nodes are classifications and the edges reflect (by weights, or color maybe) the % of observations moving from classifications from one model to the next. E.g. ![enter image description here](http://i.stack.imgur.com/muEii.png) UPDATE: Having some progress with the igraph package. Starting from the code above... poLCA results recycle the same numbers to describe class membership, so you need to do a bit of recoding.               N<-ncol(results)      n<-0     for(i in 2:N) {     results[,i]<- (results[,i])+((i-1)+n)     n<-((i-1)+n)     }      Then you need to get all the cross-tabulations and their frequencies, and rbind them into one matrix defining all the edges. There is probably a much more elegant way to do this.               results <-as.data.frame(results)          g1           <- count(results,c("x1", "x2"))          g2           <- count(results,c("x2", "x3"))     colnames(g2) <- c("x1", "x2", "freq")          g3           <- count(results,c("x3", "x4"))     colnames(g3) <- c("x1", "x2", "freq")          g4           <- count(results,c("x4", "x5"))     colnames(g4) <- c("x1", "x2", "freq")          results <- rbind(g1, g2, g3, g4)          library(igraph)          g1 <- graph.data.frame(results, directed=TRUE)          plot.igraph(g1, layout=layout.reingold.tilford)      ![enter image description here](http://i.stack.imgur.com/iCJ2Z.png) Time to play more with the igraph options I guess.


title: Can you run a multiple regression where variables use different Likert scales?

text: I have 8 predictor variables and 1 criterion variable, all variables are measured by 9 different constructs, some are 0-4 Likert scales, some are 1-5 Likert scales, can I use their mean scores and run a multiple regression?


title: Which pdf is this?

text: I have a math coursework where we are supposed to find the MLEs of the following distribution: $$ p(x;\alpha) = \frac{1+\alpha x}{2}, x \in [-1,1], \alpha\in[-1,1] $$ I'm just curious as to what it is as I have never seen it before. It is apparently used for modelling something about muon disintegration in physics. Sorry if this should have been posted in the maths or physics sections.


title: Tests for spatial stationarity (homogeneity)?

text: There are many models for spatial point patterns and spatial marked point patterns that assume spatial homogeneity or stationarity. i) Is there a statistical test for determining this, where the null is that of stationarity? ii) Is there a procedural way of looking at the data to gauge/have a say at this, by say using variograms etc?


title: Lognormal Standard Error

text: What is the Standard Error of the Lognormal distribution? I am particularly interested in comparing two probabilities from the distribution. I have the two proportions based on experiments, in general these proportions follow a lognormal fit. The standard error of differences from a normal distribution is straightforward (I've seen one other version of this which I think can see on wikipedia): $SE = \sqrt{\frac{p (1-p)}{n} + \frac{q (1-q)}{m}}$ where $p, q$ are probabilities and $n, m$ sample sizes, expected to be from a normally distributed sample. How would this be extended to the lognormal? Edit: The basic idea is to learn if the difference between two probabilities would be significant or not. I don't think I'm expressing my question correctly as I'm not a strict statistician.


title: When can I leave counterbalancing factors (between participants) in an ANOVA out of a write up?

text: The following is a problem that I keep encountering in different forms. For now, I give one example, but I would also be interested in the answers to the more general question. I have a design with a reaction time based dependent variable $Y$, two hypothesis-relevant between participants factors ($A$ and $B$) and three counterbalancing factor ($C$ [stimuli used in different conditions], $D$ [the hands assigned to different responses], and $E$ [the order of different parts of the reaction time measure]). Except for the DV, none of the variables is continuous. For my research question it is interesting:   1. whether the DV x is overall above zero    2. whether factor $A$ has a significant influence on the DV    3. whether $A$ and $B$ interact in influencing the DV  I perform the analyses with ANOVA. Writing up the results of the whole analysis with all five factors takes up a lot of journal space, so I wonder in which case I can leave out the counterbalancing factors from the analysis:   * Only if there is no effect or interaction whatsoever involving the counterbalancing factors?    * Only if there is no interaction between a counterbalancing factor and one of the hypothesis-relevant effects (e.g., no interaction with $A$, no three-way interaction with $A$ and $B$?    * Does it matter whether an interaction between a counterbalancing factor and a hypothesis-relevant effect is ordinal or disordinal (i.e., whether the effect is "turned around" in one condition)?    * And if I encounter a relevant effect or interaction with a counterbalancing factor, how should I continue: reporting the ANOVA including this factor, or should I use a covariate analysis?    * Would the approach be the same if I had a within-participants design?


title: Can we use IPS test with Kao test?

text: I have a question about using panel unit root test and panel co integration test. I used IPS unit root test then I realized my variables are I(1) can I now use kao co-integration test for checking long run relationship? I think IPS is for heterogeneous panels & kao is for homogeneous panels so we can not use both of them? Is my vision true? if anybody can help me I really appreciate it.


title: With simple random sampling, how to approximate variance of R=avg(Y)/avg(X)?

text: Recenly I am reading "Mathematical statistics and data analysis" written by Rice myself. At page 207, theorem A said: With simple random sampling, approimxate variance of $R=\frac{\overline{Y}}{\overline{X}}$ is $Var(R)\approx\dfrac{r^2\sigma^2_{\overline{X}}+\sigma^2_{\overline{Y}}-2r\sigma^2_{\overline{XY}}}{\mu^2_{x}}$ where $r=\dfrac{\sum^N_{i=1}{y_i}}{\sum^N_{i=1}{x_i}}$. Would anyone give me some hint about the procedure of the proof? Thanks a lot!


title: Learning from small number of sources

text: I have data from two sources; the two data sets were obtained from different geographical locations and at different historical times. I do not want to pool them into one data set sample and lose information about differing values of parameters that generated those collections, but I also cannot use Bayesian hierarchical model, since there would be just two unobservables from which I could not learn anything about hyperparameters. Is there any possible models that could be used in such a problems? I know that the question is a bit vague, but I am open to various possibilities.   EDIT   The data are counts, distributed as Poisson with intensities $\lambda _{1}$ and $\lambda _{2}$. Suppose I have two data samples that comes from similar objects and represents number of failures over some period of time. I want to estimate both rates. I don't want to pool them because, for example, the operating conditions (environmental conditions) were different for those two data sources so that intensities were not identical. But I want to be able to borrow some information from each data sample like in hierarchical Bayesian model. But, as I already said, for hierarchical model there is just not enough data sources.


title: Using the Fisher's sign test to select the better of two classifiers

text: I have two classifiers that classify some data points in the test set into either class 0 or 1. They each might have the same or different predictions for each test item. I want to compare them using Fisher's sign test in order to tell which is performing better. I am not sure how to do that and I am not able to understand how Fisher's sign test works.


title: Confidence interval proportion estimate

text: You are interested in estimating the proportion of female students at the University. Determine the sample size to be 95% Confident that the estimated value will be within +/- 0.098 of the real value. There is no prior information about the population proportion The answer is 100. How is that correct?


title: Normal distribution calculation problem

text: Stats question I need help with... Any help would be great thanks > Scores for a civil service exam are normally distributed, with a standard > deviation of 6.9. To be eligible for civil service employment, you must > score in the top 5%. If the lowest score you can earn and still be eligible > for employment is an 84, what is the mean score for the exam? Round your > answer to the nearest tenth.